From c3c1fa7f12fb35268ca70d65f9b15bf22023f2b0 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Thu, 24 Sep 2026 10:54:27 +0200 Subject: [PATCH 1/4] Fill VSMSolution with everything solve() returned, and remove solve() solve! now fills lift/drag/side, cl/cd/cs and their inflow-projected distributions, alpha_uncorrected, va_ref_vec, q_ref, rey, the areas, span, projected aspect ratio and both centers of pressure, allocation-free. Plotting, generate_polar_data, examples, docs and tests read the struct; solve() and calculate_results are gone. Co-Authored-By: Claude Opus 5.5 --- docs/src/examples.md | 14 +- docs/src/functions.md | 2 - docs/src/private_functions.md | 2 + docs/src/tips_and_tricks.md | 2 +- examples/V3_kite.jl | 2 +- examples/bench.jl | 3 - examples/billowing.jl | 12 +- examples/obj_to_yaml_kite.jl | 2 +- examples/pyramid_model.jl | 2 +- examples/readme_figure.jl | 2 +- examples/rectangular_wing.jl | 18 +- examples/stall_model.jl | 6 +- ext/VortexStepMethodMakieExt.jl | 42 +-- src/VortexStepMethod.jl | 5 +- src/body_aerodynamics.jl | 328 ++---------------- src/plotting_helpers.jl | 26 +- src/solver.jl | 149 ++++---- src/wing_geometry.jl | 13 +- test/bench.jl | 55 +-- .../test_body_aerodynamics.jl | 39 +-- test/plotting/test_plotting.jl | 4 +- test/solver/test_attached_trailed_force.jl | 7 +- test/solver/test_moment_units.jl | 14 - test/solver/test_solver.jl | 40 ++- test/solver/test_viscous_drag_correction.jl | 8 - test/solver/test_wing_directions.jl | 25 +- test/verification/test_verification.jl | 6 +- 27 files changed, 262 insertions(+), 566 deletions(-) diff --git a/docs/src/examples.md b/docs/src/examples.md index 03cc167c..6b9cddab 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -91,20 +91,20 @@ julia> vsm_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic #### Step 6: Solve using both methods ```julia -julia> results_llt = solve(llt_solver, body_aero) -julia> results_vsm = solve(vsm_solver, body_aero) +julia> results_llt = solve!(llt_solver, body_aero) +julia> results_vsm = solve!(vsm_solver, body_aero) ``` ##### Print results comparison ```julia julia> println("\nLifting Line Theory Results:") -julia> println("CL = $(round(results_llt["cl"], digits=4))") -julia> println("CD = $(round(results_llt["cd"], digits=4))") +julia> println("CL = $(round(results_llt.cl, digits=4))") +julia> println("CD = $(round(results_llt.cd, digits=4))") julia> println("\nVortex Step Method Results:") -julia> println("CL = $(round(results_vsm["cl"], digits=4))") -julia> println("CD = $(round(results_vsm["cd"], digits=4))") -julia> println("Projected area = $(round(results_vsm["projected_area"], digits=4)) m²") +julia> println("CL = $(round(results_vsm.cl, digits=4))") +julia> println("CD = $(round(results_vsm.cd, digits=4))") +julia> println("Projected area = $(round(results_vsm.projected_area, digits=4)) m²") ``` #### Step 7: Plot combined analysis diff --git a/docs/src/functions.md b/docs/src/functions.md index 8f934a66..8d363a23 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -98,14 +98,12 @@ CurrentModule = VortexStepMethod set_va! apparent_wind section_pitch_rate -solve solve! solve_base! reinit!(body_aero::BodyAerodynamics{P, W, T}) where {P, W, T} linearize stability_derivatives trim_angle -calculate_results ``` ## Main Plotting Functions diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index ea497f91..a7912813 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -64,6 +64,8 @@ flow_curvature_cm spanwise_flow_drag panel_force_directions prescribed_va_directions +inflow_loads! +find_center_of_pressure! panel_moment panel_couple_force panel_loads diff --git a/docs/src/tips_and_tricks.md b/docs/src/tips_and_tricks.md index 386f7ea2..3a49b033 100644 --- a/docs/src/tips_and_tricks.md +++ b/docs/src/tips_and_tricks.md @@ -37,7 +37,7 @@ This approach is useful for: If running the example `ram_air_kite.jl` fails, try to run the `cleanup.jl` script and then try again. Background: this example caches the calculated polars. Reading cached polars can fail after an update. ## Output formats -Currently, the `solve!()` function returns the results as [`VSMSolution`](@ref) struct. The function solve() returns a dictionary with the results. The `solve!()` function is faster, and the `solve()` contains many more entries, therefore the first function is good for integration in dynamic models and the second one better suited for aerodynamic analysis. +The `solve!()` function returns the results as a [`VSMSolution`](@ref) struct, which it overwrites on every call. Keep a result from one solve past the next by solving it with a separate `Solver`. ## Performance Calling `reinit!(body_aero; init_aero=false)` is very fast. After calling `unrefined_deform!(wing, theta_angles, delta_angles)`, you have to run `reinit!(body_aero; init_aero=false)` to apply the deformed wing to the body aerodynamics. This is in turn necessary for the linearization from deformation to aerodynamic coefficients for RAM-air kites. diff --git a/examples/V3_kite.jl b/examples/V3_kite.jl index 46b74245..af01538e 100644 --- a/examples/V3_kite.jl +++ b/examples/V3_kite.jl @@ -76,7 +76,7 @@ sideslip_deg = settings.condition.beta yaw_rate = settings.condition.yaw_rate set_va!(body_aero, settings) -results = VortexStepMethod.solve(solver, body_aero; log=true) +results = solve!(solver, body_aero; log=true) PLOT && plot_polars( solvers, diff --git a/examples/bench.jl b/examples/bench.jl index 2f17b216..ca738597 100644 --- a/examples/bench.jl +++ b/examples/bench.jl @@ -45,7 +45,6 @@ llt_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_ vsm_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=VSM) # Step 5: Solve using both methods -results_vsm = solve(vsm_solver, body_aero, nothing) sol = solve!(vsm_solver, body_aero, nothing) results_vsm_base = solve_base!(vsm_solver, body_aero, nothing) println("Rectangular wing, solve_base!:") @@ -54,8 +53,6 @@ println("Rectangular wing, solve_base!:") # time Julia: 0.35 ms Ryzen 7950x println("Rectangular wing, solve!:") @time sol = solve!(vsm_solver, body_aero, nothing) -println("Rectangular wing, solve:") -@time solve(vsm_solver, body_aero, nothing) # Create wing geometry (convert-then-load: obj -> per-section NeuralFoil POLAR_MATRICES) ram_yaml = obj_to_yaml( diff --git a/examples/billowing.jl b/examples/billowing.jl index 74a8b781..aa06fda9 100644 --- a/examples/billowing.jl +++ b/examples/billowing.jl @@ -94,15 +94,15 @@ set_va!(body_aero_flat, va_vec) set_va!(body_aero_bill, va_vec) # --- Solve and compare --- -results_flat = VortexStepMethod.solve( +results_flat = solve!( solver_flat, body_aero_flat; log=true) -results_bill = VortexStepMethod.solve( +results_bill = solve!( solver_bill, body_aero_bill; log=true) -println("\nFlat wing: CL=$(round(results_flat["cl"]; digits=4)), " * - "CD=$(round(results_flat["cd"]; digits=4))") -println("Billowed: CL=$(round(results_bill["cl"]; digits=4)), " * - "CD=$(round(results_bill["cd"]; digits=4))") +println("\nFlat wing: CL=$(round(results_flat.cl; digits=4)), " * + "CD=$(round(results_flat.cd; digits=4))") +println("Billowed: CL=$(round(results_bill.cl; digits=4)), " * + "CD=$(round(results_bill.cd; digits=4))") if PLOT # Plot geometry (flat wing) diff --git a/examples/obj_to_yaml_kite.jl b/examples/obj_to_yaml_kite.jl index 13f2db61..546ed536 100644 --- a/examples/obj_to_yaml_kite.jl +++ b/examples/obj_to_yaml_kite.jl @@ -75,7 +75,7 @@ VortexStepMethod.reinit!(body_aero) solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=VSM, rtol=1e-5, solver_type=LOOP) set_va!(body_aero, apparent_wind(deg2rad(8), 0.0, va)) -results = VortexStepMethod.solve(solver, body_aero; log=true) +results = solve!(solver, body_aero; log=true) if PLOT plot_geometry(body_aero, "Ram air kite (converted from .obj)"; is_show=true, diff --git a/examples/pyramid_model.jl b/examples/pyramid_model.jl index cd19b268..96f8dd5f 100644 --- a/examples/pyramid_model.jl +++ b/examples/pyramid_model.jl @@ -26,7 +26,7 @@ sideslip_deg = vsm_settings.condition.beta yaw_rate = vsm_settings.condition.yaw_rate # Run the solver -results = VortexStepMethod.solve(solver, body_aero; log=true) +results = solve!(solver, body_aero; log=true) # Using plotting modules, to create more comprehensive plots PLOT = true diff --git a/examples/readme_figure.jl b/examples/readme_figure.jl index e270a50b..074ac00b 100644 --- a/examples/readme_figure.jl +++ b/examples/readme_figure.jl @@ -13,7 +13,7 @@ labels = ["VSM Julia", "CFD Re=5e5", "CFD Re=10e5", "VSM Python Re=5e5", "Wind tunnel Re=5e5"] set_va!(body_aero, settings) -results = VortexStepMethod.solve(solver, body_aero) +results = solve!(solver, body_aero) plot_combined_analysis(solver, body_aero, results; labels, diff --git a/examples/rectangular_wing.jl b/examples/rectangular_wing.jl index 2e45aa99..bb64e625 100644 --- a/examples/rectangular_wing.jl +++ b/examples/rectangular_wing.jl @@ -49,19 +49,19 @@ llt_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_ vsm_solver = Solver(wing.n_panels, wing.n_unrefined_sections; aerodynamic_model_type=VSM) # Step 5: Solve using both methods -results_llt = solve(llt_solver, body_aero) -results_vsm = solve(vsm_solver, body_aero) -@time solve(llt_solver, body_aero) -@time solve(vsm_solver, body_aero) +results_llt = solve!(llt_solver, body_aero) +results_vsm = solve!(vsm_solver, body_aero) +@time solve!(llt_solver, body_aero) +@time solve!(vsm_solver, body_aero) # Print results comparison println("\nLifting Line Theory Results:") -println("CL = $(round(results_llt["cl"], digits=4))") -println("CD = $(round(results_llt["cd"], digits=4))") +println("CL = $(round(results_llt.cl, digits=4))") +println("CD = $(round(results_llt.cd, digits=4))") println("\nVortex Step Method Results:") -println("CL = $(round(results_vsm["cl"], digits=4))") -println("CD = $(round(results_vsm["cd"], digits=4))") -println("Projected area = $(round(results_vsm["projected_area"], digits=4)) m²") +println("CL = $(round(results_vsm.cl, digits=4))") +println("CD = $(round(results_vsm.cd, digits=4))") +println("Projected area = $(round(results_vsm.projected_area, digits=4)) m²") # Step 6: Plot geometry PLOT && plot_geometry( diff --git a/examples/stall_model.jl b/examples/stall_model.jl index 3b7fb226..b9796c77 100644 --- a/examples/stall_model.jl +++ b/examples/stall_model.jl @@ -77,9 +77,9 @@ PLOT && plot_geometry( ) # Solving and plotting distributions -results = solve(vsm_solver, body_aero) -@time results_with_stall = solve(VSM_with_stall_correction, body_aero) -@time results_with_stall = solve(VSM_with_stall_correction, body_aero) +results = solve!(vsm_solver, body_aero) +@time results_with_stall = solve!(VSM_with_stall_correction, body_aero) +@time results_with_stall = solve!(VSM_with_stall_correction, body_aero) CAD_y_coordinates = [panel.aero_center[2] for panel in body_aero.panels] diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index 111e0fcb..a1c6306c 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/ext/VortexStepMethodMakieExt.jl @@ -764,7 +764,7 @@ Makie implementation of [`plot_distribution`](@ref). # Arguments - `y_coordinates_list`: List of spanwise coordinates -- `results_list`: List of result dictionaries +- `results_list`: List of [`VSMSolution`](@ref)s - `label_list`: List of labels for different results # Keyword arguments @@ -811,39 +811,39 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la # Plot CL for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - value = round(results["cl"], digits=2) - lines!(ax_cl, Vector(y_coords), Vector(results["cl_distribution"]), + value = round(results.cl, digits=2) + lines!(ax_cl, Vector(y_coords), Vector(results.cl_distribution), label="$label CL: $value") end # Plot CD for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - value = round(results["cd"], digits=2) - lines!(ax_cd, Vector(y_coords), Vector(results["cd_distribution"]), + value = round(results.cd, digits=2) + lines!(ax_cd, Vector(y_coords), Vector(results.cd_distribution), label="$label CD: $value") end # Plot Gamma for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - lines!(ax_gamma, Vector(y_coords), Vector(results["gamma_distribution"]), + lines!(ax_gamma, Vector(y_coords), Vector(results.gamma_distribution), label=label) end # Plot alpha geometric for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - lines!(ax_alpha_geo, Vector(y_coords), rad2deg.(Vector(results["alpha_geometric"])), + lines!(ax_alpha_geo, Vector(y_coords), rad2deg.(Vector(results.alpha_geometric_dist)), label=label) end # Plot alpha at ac for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - lines!(ax_alpha_ac, Vector(y_coords), rad2deg.(Vector(results["alpha_at_ac"])), + lines!(ax_alpha_ac, Vector(y_coords), rad2deg.(Vector(results.alpha_dist)), label=label) end # Plot alpha uncorrected for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - lines!(ax_alpha_unc, Vector(y_coords), rad2deg.(Vector(results["alpha_uncorrected"])), + lines!(ax_alpha_unc, Vector(y_coords), rad2deg.(Vector(results.alpha_uncorrected)), label=label) end @@ -852,12 +852,12 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la components = ["x", "y", "z"] for (idx, (ax, comp)) in enumerate(zip(force_axes, components)) for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - forces = results["F_distribution"][idx, :] + forces = results.f_body_3D[idx, :] if length(y_coords) != length(forces) @warn "Dimension mismatch" length(y_coords) length(forces) comp continue end - total_force = round(results["F$comp"], digits=2) + total_force = round(results.force[idx], digits=2) lines!(ax, Vector(y_coords), Vector(forces), label="$label ΣF$comp: $total_force N") end @@ -1128,7 +1128,8 @@ Makie implementation of [`plot_combined_analysis`](@ref). # Arguments - `solver`: Aerodynamic solver - `body_aero`: BodyAerodynamics object -- `results`: Solution dictionary from solve() +- `results`: [`VSMSolution`](@ref) of each solver, used for the spanwise plots when + `angle_of_attack_for_spanwise_distribution` is `nothing` # Keyword arguments - `labels`: Optional label string or label vector. If a vector with length @@ -1232,8 +1233,7 @@ function VortexStepMethod.plot_combined_analysis( omega_old = copy(ba.omega) set_va!(ba, [cos(α_span) * cos(β_span), sin(β_span), sin(α_span)] * va) - results_spanwise_list[i] = solve(s, ba, - s.sol.gamma_distribution) + results_spanwise_list[i] = solve!(s, ba) set_va!(ba, va_vec_old, omega_old) end end @@ -1363,18 +1363,18 @@ function VortexStepMethod.plot_combined_analysis( color = colors[mod1(si, length(colors))] y_si = [panel.aero_center[2] for panel in body_aeros[si].panels] - lines!(ax_cl, Vector(y_si), Vector(rs["cl_distribution"]); color) - lines!(ax_cd, Vector(y_si), Vector(rs["cd_distribution"]); color) + lines!(ax_cl, Vector(y_si), Vector(rs.cl_distribution); color) + lines!(ax_cd, Vector(y_si), Vector(rs.cd_distribution); color) lines!(ax_gamma, Vector(y_si), - Vector(rs["gamma_distribution"]); label=lbl, color) + Vector(rs.gamma_distribution); label=lbl, color) lines!(ax_alpha_geo, Vector(y_si), - rad2deg.(Vector(rs["alpha_geometric"])); color) + rad2deg.(Vector(rs.alpha_geometric_dist)); color) lines!(ax_alpha_ac, Vector(y_si), - rad2deg.(Vector(rs["alpha_at_ac"])); color) + rad2deg.(Vector(rs.alpha_dist)); color) lines!(ax_alpha_unc, Vector(y_si), - rad2deg.(Vector(rs["alpha_uncorrected"])); color) + rad2deg.(Vector(rs.alpha_uncorrected)); color) for (idx, ax) in enumerate((ax_fx, ax_fy, ax_fz)) - lines!(ax, Vector(y_si), Vector(rs["F_distribution"][idx, :]); color) + lines!(ax, Vector(y_si), Vector(rs.f_body_3D[idx, :]); color) end end diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 6f473f1c..f0b085ac 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -12,7 +12,7 @@ using Measures using LaTeXStrings using NonlinearSolve using SciMLBase -import NonlinearSolve: solve, solve! +import NonlinearSolve: solve! using Interpolations using Parameters using Printf: @sprintf @@ -32,10 +32,9 @@ export airfoil_solver, alpha_range, delta_range, reynolds, rotation_matrix export slice_args, preview_args export ObjWing, Section, Wing, refine!, reinit! export BodyAerodynamics -export Solver, VSMSolution, linearize, solve, solve!, solve_base!, calc_forces! +export Solver, VSMSolution, linearize, solve!, solve_base!, calc_forces! export stability_derivatives, trim_angle export SolveFailure -export calculate_results export add_section!, apparent_wind, set_va!, section_pitch_rate export calculate_projected_area, calculate_span export MVec3 diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 6cfdceec..0e73c49a 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -663,9 +663,9 @@ end denominator = dot3(plane_normal, f_unit) abs(denominator) < tol && return nothing λ = numerator / denominator - return MVec3(x_cp[1] + λ*f_unit[1], - x_cp[2] + λ*f_unit[2], - x_cp[3] + λ*f_unit[3]) + return SVector{3}(x_cp[1] + λ*f_unit[1], + x_cp[2] + λ*f_unit[2], + x_cp[3] + λ*f_unit[3]) end @inline function point_in_triangle(pt, v0, v1, v2; tol=1e-8) @@ -717,7 +717,15 @@ end return (s >= -tol) && (t >= -tol) && (s + t <= 1 + tol) end -function find_center_of_pressure( +""" + find_center_of_pressure!(center_of_pressure, body_aero::BodyAerodynamics, force, + moment, reference_point; force_tol=1e-12) + +Set `center_of_pressure` to the first point where the line of action of `force` and +`moment` about `reference_point` crosses a panel, or to `NaN` where it crosses none. +""" +function find_center_of_pressure!( + center_of_pressure, body_aero::BodyAerodynamics, force, moment, @@ -728,9 +736,10 @@ function find_center_of_pressure( M0 = moment r0 = reference_point F_norm_sq = dot3(F, F) + center_of_pressure .= NaN # Treat near-zero forces as "CoP undefined" if !(isfinite(F_norm_sq)) || F_norm_sq ≤ force_tol^2 - return nothing + return center_of_pressure end wv = body_aero.work_vectors @@ -771,25 +780,22 @@ function find_center_of_pressure( r0_moment, f_unit, cross_tmp, normal) if !isnothing(intersection) && point_in_quad(intersection, corners) - return MVec3(intersection) + center_of_pressure .= intersection + return center_of_pressure end end end - - @warn "No intersection found with any panel " * - "in center-of-pressure calculation." - return nothing + return center_of_pressure end -function compute_panel_center_of_pressures( +function compute_panel_center_of_pressures!( + panel_cp_locations, body_aero::BodyAerodynamics, f_distribution::AbstractMatrix, m_distribution::AbstractMatrix, reference_point ) - n = length(body_aero.panels) - panel_cp_locations = Vector{MVec3}(undef, n) - for i in 1:n + for i in eachindex(body_aero.panels) panel = body_aero.panels[i] @views F = f_distribution[:, i] @views M_ref = m_distribution[:, i] @@ -800,7 +806,7 @@ function compute_panel_center_of_pressures( # Guard against non-finite forces and near-zero forces if !all(isfinite, F) || dot3(F, F) ≤ 1e-24 - panel_cp_locations[i] = MVec3(ac) + panel_cp_locations[i] .= ac continue end @@ -825,13 +831,13 @@ function compute_panel_center_of_pressures( cz*span_dir[3] if abs(F_perp_mag) < 1e-12 - panel_cp_locations[i] = MVec3(ac) + panel_cp_locations[i] .= ac continue end lever = clamp(m_pitch / F_perp_mag, -0.25 * c, 0.75 * c) - panel_cp_locations[i] = MVec3( + panel_cp_locations[i] .= ( ac[1] + lever*chord_dir[1], ac[2] + lever*chord_dir[2], ac[3] + lever*chord_dir[3]) @@ -865,294 +871,6 @@ Lift and side unit vectors, `(; dir_lift, dir_side)`, of inflow `va` on a wing a return (; dir_lift, dir_side=cross(dir_lift, va) / norm(va)) end -""" - calculate_results(body_aero::BodyAerodynamics, gamma_new, reference_point, density, - core_radius_fraction, mu, alpha_dist, v_rel_dist, chord_dist, - x_airf_dist, z_airf_dist, va_vec_dist, va_dist, va_unit_dist, - panels::Vector{<:Panel}, is_only_f_and_gamma_output::Bool; - correct_aoa=false, flow_curvature=false, - is_with_viscous_drag_correction=false, - is_with_attached_trailed_force=false, v_span_dist=nothing) - -Calculate final aerodynamic results. Reference point is in the kite body (KB) frame. - -`flow_curvature` adds [`flow_curvature_cm`](@ref) to every section moment, read -from `body_aero.omega`. `is_with_viscous_drag_correction` adds -[`spanwise_flow_drag`](@ref) to every section, from the velocity along `y_airf` in -`v_span_dist`. `is_with_attached_trailed_force` adds [`attached_trailed_loads`](@ref). - -Returns: - Dict: Results including forces, coefficients and distributions -""" -function calculate_results( - body_aero::BodyAerodynamics, - gamma_new, - reference_point, - density, - core_radius_fraction, - mu, - alpha_dist, - v_rel_dist, - chord_dist, - x_airf_dist, - z_airf_dist, - va_vec_dist, - va_dist, - va_unit_dist, - panels::Vector{<:Panel}, - is_only_f_and_gamma_output::Bool; - correct_aoa::Bool=false, - flow_curvature::Bool=false, - is_with_viscous_drag_correction::Bool=false, - is_with_attached_trailed_force::Bool=false, - v_span_dist=nothing, -) - - n_panels = length(panels) - if length(body_aero.cache) < 15 - append!(body_aero.cache, [LazyBufferCache() for _ in 1:(15 - length(body_aero.cache))]) - end - - cl_dist = body_aero.cache[5][alpha_dist] - cd_dist = body_aero.cache[6][alpha_dist] - cm_dist = body_aero.cache[7][alpha_dist] - panel_width_dist = body_aero.cache[8][alpha_dist] - alpha_corrected = body_aero.cache[9][alpha_dist] - cl_prescribed_va = body_aero.cache[10][alpha_dist] - cd_prescribed_va = body_aero.cache[11][alpha_dist] - cs_prescribed_va = body_aero.cache[12][alpha_dist] - f_body_3D = body_aero.cache[13][alpha_dist, (3, length(alpha_dist))] - m_body_3D = body_aero.cache[14][alpha_dist, (3, length(alpha_dist))] - alpha_geometric = body_aero.cache[15][alpha_dist] - - fill!(f_body_3D, 0.0) - fill!(m_body_3D, 0.0) - - # Calculate coefficients and geometric AoA for each panel - for (i, panel) in enumerate(panels) - cl_dist[i] = calculate_cl(panel, alpha_dist[i]) - cd_dist[i], cm_dist[i] = calculate_cd_cm( - panel, alpha_dist[i]) - if flow_curvature - cm_dist[i] += flow_curvature_cm( - body_aero.pitch_rate_dist[i], chord_dist[i], v_rel_dist[i]) - end - panel_width_dist[i] = panel.width - va = va_dist[i] - x_norm = norm3(panel.x_airf) - z_norm = norm3(panel.z_airf) - if va == 0.0 || x_norm == 0.0 || z_norm == 0.0 - alpha_geometric[i] = NaN - else - inv_va = 1.0 / va - v_tangential = -dot3(panel.x_airf, panel.va_vec) * - inv_va / x_norm - v_normal = -dot3(panel.z_airf, panel.va_vec) * - inv_va / z_norm - alpha_geometric[i] = atan(-v_normal, -v_tangential) - end - end - - # Calculate alpha corrections based on model type - if correct_aoa - update_effective_angle_of_attack!( - alpha_corrected, - body_aero, - gamma_new, - core_radius_fraction, - z_airf_dist, - x_airf_dist, - va_vec_dist, - va_dist, - va_unit_dist - ) - else - alpha_corrected .= alpha_dist - end - - area_all_panels = 0.0 - lift_wing_3D_sum = 0.0 - drag_wing_3D_sum = 0.0 - side_wing_3D_sum = 0.0 - - reference_spanwise = SVector{3}(body_aero.wings[1].spanwise_direction) - va_ref_vec = MVec3(0.0, 0.0, 0.0) - weighted_speed_sq = 0.0 - total_area = 0.0 - @inbounds for i in 1:n_panels - area_i = chord_dist[i] * panel_width_dist[i] - total_area += area_i - speed_i = va_dist[i] - weighted_speed_sq += area_i * speed_i^2 - va_ref_vec[1] += area_i * va_vec_dist[i, 1] - va_ref_vec[2] += area_i * va_vec_dist[i, 2] - va_ref_vec[3] += area_i * va_vec_dist[i, 3] - end - total_area > 0.0 || throw(ArgumentError( - "Total panel area must be positive.")) - reference_speed = sqrt(weighted_speed_sq / total_area) - direction_norm = norm3(va_ref_vec) - if direction_norm <= 0.0 - va_ref_vec .= (1.0, 0.0, 0.0) - direction_norm = 1.0 - end - @inbounds for k in 1:3 - va_ref_vec[k] = va_ref_vec[k] / direction_norm * - reference_speed - end - va_ref = norm3(va_ref_vec) - va_ref > 0.0 || throw(ArgumentError( - "Reference freestream magnitude must be positive.")) - va_ref_unit = SVector{3}(va_ref_vec) / va_ref - reference_dirs = prescribed_va_directions(SVector{3}(va_ref_vec), reference_spanwise) - all(isfinite, reference_dirs.dir_lift) || throw(ArgumentError( - "Reference lift direction is undefined because " * - "reference flow is parallel to spanwise direction.")) - q_ref = 0.5 * density * va_ref^2 - - for (wing_idx, wing) in enumerate(body_aero.wings) - spanwise_unit = SVector{3}(wing.spanwise_direction) - for i in panel_range(body_aero, wing_idx) - panel = panels[i] - panel_area = panel.chord * panel.width - area_all_panels += panel_area - - axes = panel_axes(panel) - dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit) - c_span = 0.0 - if is_with_viscous_drag_correction - viscous = spanwise_flow_drag(v_rel_dist[i], v_span_dist[i], panel.chord, - density, mu) - cd_dist[i] += viscous.delta_cd - c_span = viscous.c_span - end - loads = panel_loads(axes, dirs, - dynamic_pressure(density, density, v_rel_dist[i]), - cl_dist[i], cd_dist[i], cm_dist[i]; c_span) - (; force, moment) = panel_force_moment(body_aero, i, loads, axes.y_airf, - gamma_new, density, core_radius_fraction, reference_point, - is_with_attached_trailed_force) - - va_panel = va_dist[i] - va_panel > 0.0 || throw(ArgumentError( - "Panel $i has non-positive apparent " * - "velocity magnitude.")) - q_panel = 0.5 * density * va_panel^2 - panel_va = SVector{3}(panel.va_vec) - inv_va_panel = 1.0 / va_panel - drag_prescribed_va = dot(force, panel_va) * inv_va_panel - wing_dirs = prescribed_va_directions(panel_va, spanwise_unit) - body_dirs = prescribed_va_directions(panel_va, reference_spanwise) - - lift_wing_3D_sum += dot(force, body_dirs.dir_lift) * - dot(body_dirs.dir_lift, reference_dirs.dir_lift) - drag_wing_3D_sum += drag_prescribed_va * - (dot(panel_va, va_ref_unit) * inv_va_panel) - side_wing_3D_sum += dot(force, body_dirs.dir_side) * - dot(body_dirs.dir_side, reference_dirs.dir_side) - - inv_q_area = 1.0 / (q_panel * panel_area) - cl_prescribed_va[i] = dot(force, wing_dirs.dir_lift) * inv_q_area - cd_prescribed_va[i] = drag_prescribed_va * inv_q_area - cs_prescribed_va[i] = dot(force, wing_dirs.dir_side) * inv_q_area - - @inbounds for k in 1:3 - f_body_3D[k, i] = force[k] - m_body_3D[k, i] = moment[k] - end - end - end - - if is_only_f_and_gamma_output - return Dict{String,Any}( - "F_distribution" => copy(f_body_3D), - "gamma_distribution" => gamma_new - ) - end - - # Calculate wing geometry properties - projected_area = body_aero.projected_area - wing_span = calculate_span(body_aero.wings, reference_spanwise) - aspect_ratio_projected = wing_span^2 / projected_area - - # Calculate Reynolds number - c_ref = body_aero.c_ref - reynolds_number = density * va_ref * c_ref / mu - - force_total = body_aero.work_vectors[9] - moment_total = body_aero.work_vectors[10] - force_total .= 0.0 - moment_total .= 0.0 - @inbounds for i in 1:n_panels - force_total[1] += f_body_3D[1, i] - force_total[2] += f_body_3D[2, i] - force_total[3] += f_body_3D[3, i] - moment_total[1] += m_body_3D[1, i] - moment_total[2] += m_body_3D[2, i] - moment_total[3] += m_body_3D[3, i] - end - center_of_pressure = try - find_center_of_pressure(body_aero, force_total, moment_total, reference_point) - catch err - @warn "Center-of-pressure calculation failed: $(err)" - nothing - end - panel_cp_locations = compute_panel_center_of_pressures( - body_aero, - f_body_3D, - m_body_3D, - reference_point - ) - - # Create results dictionary - results = Dict{String,Any}( - "Fx" => force_total[1], - "Fy" => force_total[2], - "Fz" => force_total[3], - "Mx" => moment_total[1], - "My" => moment_total[2], - "Mz" => moment_total[3], - "lift" => lift_wing_3D_sum, - "drag" => drag_wing_3D_sum, - "side" => side_wing_3D_sum, - "cl" => lift_wing_3D_sum / (q_ref * projected_area), - "cd" => drag_wing_3D_sum / (q_ref * projected_area), - "cs" => side_wing_3D_sum / (q_ref * projected_area), - "cmx" => moment_total[1] / (q_ref * projected_area * c_ref), - "cmy" => moment_total[2] / (q_ref * projected_area * c_ref), - "cmz" => moment_total[3] / (q_ref * projected_area * c_ref), - "cl_distribution" => copy(cl_prescribed_va), - "cd_distribution" => copy(cd_prescribed_va), - "cs_distribution" => copy(cs_prescribed_va), - "F_distribution" => copy(f_body_3D), - "M_distribution" => copy(m_body_3D), - "cfx" => (force_total[1] / (q_ref * projected_area)), - "cfy" => (force_total[2] / (q_ref * projected_area)), - "cfz" => (force_total[3] / (q_ref * projected_area)), - "alpha_at_ac" => copy(alpha_corrected), - "alpha_uncorrected" => alpha_dist, - "alpha_geometric" => copy(alpha_geometric), - "gamma_distribution" => gamma_new, - "area_all_panels" => area_all_panels, - "projected_area" => projected_area, - "wing_span" => wing_span, - "aspect_ratio_projected" => aspect_ratio_projected, - "Rey" => reynolds_number, - "q_ref" => q_ref, - "va_ref_vec" => va_ref_vec, - "center_of_pressure" => center_of_pressure, - "panel_cp_locations" => panel_cp_locations - ) - - @debug "Results summary:" cl=results["cl"] cd=results["cd"] cs=results["cs"] - @debug "Forces:" lift=lift_wing_3D_sum drag=drag_wing_3D_sum side=side_wing_3D_sum - @debug "Areas:" total=area_all_panels projected=projected_area - @debug "Aspect ratio:" ar=aspect_ratio_projected - - return results -end - - """ set_va!(body_aero::BodyAerodynamics, va_vec::VelVector, omega=zeros(MVec3); reference_point=body_aero.reference_point) diff --git a/src/plotting_helpers.jl b/src/plotting_helpers.jl index bbc61a72..1f9fdf7a 100644 --- a/src/plotting_helpers.jl +++ b/src/plotting_helpers.jl @@ -123,21 +123,17 @@ function generate_polar_data( set_va!(body_aero, [cos(α) * cos(β), sin(β), sin(α)] * va) - results = solve(solver, body_aero, - gamma_distribution[i, :]) - - cl[i] = results["cl"] - cd[i] = results["cd"] - cs[i] = results["cs"] - cmx[i] = get(results, "cmx", NaN) - cmy[i] = get(results, "cmy", NaN) - cmz[i] = get(results, "cmz", NaN) - gamma_distribution[i, :] = - results["gamma_distribution"] - cl_distribution[i, :] = results["cl_distribution"] - cd_distribution[i, :] = results["cd_distribution"] - cs_distribution[i, :] = results["cs_distribution"] - reynolds_number[i] = results["Rey"] + sol = solve!(solver, body_aero, gamma_distribution[i, :]) + + cl[i] = sol.cl + cd[i] = sol.cd + cs[i] = sol.cs + cmx[i], cmy[i], cmz[i] = sol.moment_coeffs + gamma_distribution[i, :] = sol.gamma_distribution + cl_distribution[i, :] = sol.cl_distribution + cd_distribution[i, :] = sol.cd_distribution + cs_distribution[i, :] = sol.cs_distribution + reynolds_number[i] = sol.rey end polar_data = [ diff --git a/src/solver.jl b/src/solver.jl index 9357dfde..7779e616 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -24,6 +24,24 @@ Struct for storing the solution of the [`solve!`](@ref) function. Must contain a - moment::MVec3: Aerodynamic moments [Mx, My, Mz] around the reference point [Nm] - force_coeffs::MVec3: Aerodynamic force coefficients [CFx, CFy, CFz] [-] - `moment_coeffs`::MVec3: Aerodynamic moment coefficients [CMx, CMy, CMz] [-] +- `lift`, `drag`, `side`: Total force along the lift, drag and side directions of the + reference inflow `va_ref_vec` [N] +- `cl`, `cd`, `cs`: `lift`, `drag` and `side` divided by `q_ref * projected_area` [-] +- `cl_distribution`, `cd_distribution`, `cs_distribution`::Vector{Float64}: Panel force along + the lift, drag and side directions of the panel's own inflow, divided by its dynamic + pressure and area [-]; unlike `cl_dist` and `cd_dist`, these include every force term +- `alpha_uncorrected`::Vector{Float64}: Angle of attack of each panel at its evaluation + point, before the aerodynamic-center correction [rad] +- `va_ref_vec`::MVec3: Area-weighted reference inflow velocity [m/s] +- `q_ref`: Dynamic pressure of `va_ref_vec` [Pa] +- `rey`: Reynolds number of `va_ref_vec` on the reference chord [-] +- `area_all_panels`: Sum of the panel areas [m²] +- `projected_area`: Projected area of the body [m²] +- `wing_span`: Span along the first wing's spanwise direction [m] +- `aspect_ratio_projected`: `wing_span^2 / projected_area` [-] +- `center_of_pressure`::MVec3: Point where the line of action of `force` crosses a + panel, `NaN` where it crosses none [m] +- `panel_cp_locations`::Vector{MVec3}: Center of pressure of each panel [m] - `moment_dist`::Vector{Float64}: Pitching moments around the spanwise vector of each panel. [Nm] - `moment_coeff_dist`::Vector{Float64}: Pitching moment coefficient around the spanwise vector of each panel. [-] - `moment_unrefined_dist`::MVector{U, Float64}: Averaged moments for unrefined sections [Nm] @@ -59,8 +77,25 @@ Struct for storing the solution of the [`solve!`](@ref) function. Must contain a moment::MVector{3, T} = zeros(MVector{3, T}) force_coeffs::MVector{3, T} = zeros(MVector{3, T}) moment_coeffs::MVector{3, T} = zeros(MVector{3, T}) - center_of_pressure::Union{Nothing, MVector{3, T}} = nothing - panel_cp_locations::Vector{MVector{3, T}} = MVector{3, T}[] + lift::T = zero(T) + drag::T = zero(T) + side::T = zero(T) + cl::T = zero(T) + cd::T = zero(T) + cs::T = zero(T) + cl_distribution::Vector{T} = zeros(T, P) + cd_distribution::Vector{T} = zeros(T, P) + cs_distribution::Vector{T} = zeros(T, P) + alpha_uncorrected::Vector{T} = zeros(T, P) + va_ref_vec::MVector{3, T} = zeros(MVector{3, T}) + q_ref::T = zero(T) + rey::T = zero(T) + area_all_panels::T = zero(T) + projected_area::T = zero(T) + wing_span::T = zero(T) + aspect_ratio_projected::T = zero(T) + center_of_pressure::MVector{3, T} = zeros(MVector{3, T}) + panel_cp_locations::Vector{MVector{3, T}} = [zeros(MVector{3, T}) for _ in 1:P] moment_dist::MVector{P, T} = zeros(MVector{P, T}) moment_coeff_dist::MVector{P, T} = zeros(MVector{P, T}) moment_unrefined_dist::MVector{U, T} = zeros(MVector{U, T}) @@ -108,7 +143,7 @@ end """ Solver -Main solver structure for the Vortex Step Method.See also: [`solve`](@ref) +Main solver structure for the Vortex Step Method. See also: [`solve!`](@ref) # Attributes @@ -295,8 +330,7 @@ finite_full(x::ForwardDiff.Dual) = throw_on_fail=false) Main solving routine for the aerodynamic model. Reference point is in the kite body (KB) frame. -This version is modifying the `solver.sol` struct and is faster than the `solve` function which returns -a dictionary. +Fills and returns `solver.sol`. # Arguments: - solver::Solver: The solver to use, could be a VSM or LLT solver. See: [`Solver`](@ref) @@ -474,7 +508,8 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; end # Python parity: normalize with area-weighted reference velocity for distributed inflow. - va_ref_vec = _compute_reference_velocity_from_distribution( + va_ref_vec = solver.sol.va_ref_vec + va_ref_vec .= _compute_reference_velocity_from_distribution( solver.sol.va_vec_dist, length(panels), panel_areas @@ -482,6 +517,8 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; va_ref = norm(va_ref_vec) va_ref > 0.0 || throw(ArgumentError("Reference freestream magnitude must be positive.")) q_ref = 0.5 * density * va_ref^2 + solver.sol.q_ref = q_ref + solver.sol.rey = density * va_ref * c_ref / solver.mu moment_coeff_dist .= moment_dist ./ (q_ref * projected_area * c_ref) # Only compute unrefined arrays if there are unrefined sections @@ -579,9 +616,17 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; end solver.sol.force_coeffs .= solver.sol.force ./ (q_ref * projected_area) solver.sol.moment_coeffs .= solver.sol.moment ./ (q_ref * projected_area * c_ref) - # Keep solve! fast: center-of-pressure is only computed in solve() dictionary path. - solver.sol.center_of_pressure = nothing - empty!(solver.sol.panel_cp_locations) + solver.sol.alpha_uncorrected .= alpha_dist + solver.sol.area_all_panels = area_all_panels + solver.sol.projected_area = projected_area + reference_spanwise = body_aero.wings[1].spanwise_direction + solver.sol.wing_span = calculate_span(body_aero.wings, reference_spanwise) + solver.sol.aspect_ratio_projected = solver.sol.wing_span^2 / projected_area + inflow_loads!(solver.sol, body_aero, density) + find_center_of_pressure!(solver.sol.center_of_pressure, body_aero, solver.sol.force, + solver.sol.moment, reference_point) + compute_panel_center_of_pressures!(solver.sol.panel_cp_locations, body_aero, + solver.sol.f_body_3D, solver.sol.m_body_3D, reference_point) if converged # TODO: Check if the result if feasible if converged solver.sol.solver_status = FEASIBLE @@ -593,59 +638,45 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; end """ - solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=nothing; - log=false, reference_point=solver.reference_point) + inflow_loads!(sol::VSMSolution, body_aero::BodyAerodynamics, density) -Main solving routine for the aerodynamic model. Reference point is in the kite body (KB) frame. -See also: [`solve!`](@ref) - -# Arguments: -- solver::Solver: The solver to use, could be a VSM or LLT solver. See: [`Solver`](@ref) -- body_aero::BodyAerodynamics: The aerodynamic body. See: [`BodyAerodynamics`](@ref) -- gamma_distribution: Initial circulation vector or nothing; Length: Number of segments. [m²/s] - -# Keyword Arguments: -- log=false: If true, print the number of iterations and other info. -- reference_point=solver.reference_point - -# Returns -A dictionary with the results. +Fill `lift`, `drag`, `side`, `cl`, `cd`, `cs` and `cl_distribution`, `cd_distribution`, +`cs_distribution` of `sol` by projecting the panel forces `sol.f_body_3D` on the lift, drag +and side directions of the reference inflow `sol.va_ref_vec` and of each panel's own inflow. """ -function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=nothing; - log=false, reference_point=solver.reference_point) - reference_point_checked = check_reference_point(reference_point) - # calculate intermediate result - solve_base!(solver, body_aero, gamma_distribution; log) - - # Calculate final results as dictionary - results = calculate_results( - body_aero, - solver.lr.gamma_new, - reference_point_checked, - solver.density, - solver.core_radius_fraction, - solver.mu, - solver.lr.alpha_dist, - solver.lr.v_rel_dist, - solver.sol._chord_dist, - solver.sol._x_airf_dist, - solver.sol._z_airf_dist, - solver.sol.va_vec_dist, - solver.br.va_dist, - solver.br.va_unit_dist, - body_aero.panels, - solver.is_only_f_and_gamma_output; - correct_aoa=solver.correct_aoa, - flow_curvature=solver.flow_curvature, - is_with_viscous_drag_correction=solver.is_with_viscous_drag_correction, - is_with_attached_trailed_force=solver.is_with_attached_trailed_force, - v_span_dist=solver.lr.v_span_dist, - ) - # Attach geometric AoA (already computed in calculate_results) to solver.sol - if haskey(results, "alpha_geometric") - solver.sol.alpha_geometric_dist .= results["alpha_geometric"] +function inflow_loads!(sol::VSMSolution, body_aero::BodyAerodynamics, density) + reference_spanwise = SVector{3}(body_aero.wings[1].spanwise_direction) + va_ref = SVector{3}(sol.va_ref_vec) + va_ref_unit = va_ref / norm(va_ref) + reference_dirs = prescribed_va_directions(va_ref, reference_spanwise) + lift = drag = side = zero(sol.q_ref) + for (wing_idx, wing) in enumerate(body_aero.wings) + spanwise_unit = SVector{3}(wing.spanwise_direction) + for i in panel_range(body_aero, wing_idx) + force = SVector{3}(sol.f_body_3D[1, i], sol.f_body_3D[2, i], sol.f_body_3D[3, i]) + panel_va = SVector{3}(sol.va_vec_dist[i, 1], sol.va_vec_dist[i, 2], + sol.va_vec_dist[i, 3]) + va_panel = norm(panel_va) + panel_drag = dot(force, panel_va) / va_panel + wing_dirs = prescribed_va_directions(panel_va, spanwise_unit) + body_dirs = prescribed_va_directions(panel_va, reference_spanwise) + + lift += dot(force, body_dirs.dir_lift) * + dot(body_dirs.dir_lift, reference_dirs.dir_lift) + drag += panel_drag * dot(panel_va, va_ref_unit) / va_panel + side += dot(force, body_dirs.dir_side) * + dot(body_dirs.dir_side, reference_dirs.dir_side) + + q_area = dynamic_pressure(density, density, va_panel) * sol.panel_area_dist[i] + sol.cl_distribution[i] = dot(force, wing_dirs.dir_lift) / q_area + sol.cd_distribution[i] = panel_drag / q_area + sol.cs_distribution[i] = dot(force, wing_dirs.dir_side) / q_area + end end - return results + q_area_ref = sol.q_ref * sol.projected_area + sol.lift, sol.drag, sol.side = lift, drag, side + sol.cl, sol.cd, sol.cs = lift / q_area_ref, drag / q_area_ref, side / q_area_ref + return nothing end @inline @inbounds function calc_norm_dist!(va_dist, va_vec_dist) diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index c849c001..57214a32 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -1692,10 +1692,15 @@ Lowest and highest projection of the unrefined sections' LE and TE points on `wi `spanwise_direction`, or of all `wings` on `spanwise_direction`, as `(lo, hi)` [m]. """ function spanwise_extent(wings, spanwise_direction) - axis = normalize(spanwise_direction) - return extrema(dot(point, axis) for wing in wings - for section in wing.unrefined_sections - for point in (section.LE_point, section.TE_point)) + axis = SVector{3}(spanwise_direction) / norm(spanwise_direction) + lo, hi = Inf, -Inf + for wing in wings, section in wing.unrefined_sections + for point in (section.LE_point, section.TE_point) + distance = dot(point, axis) + lo, hi = min(lo, distance), max(hi, distance) + end + end + return lo, hi end spanwise_extent(wing::AbstractWing) = spanwise_extent((wing,), wing.spanwise_direction) diff --git a/test/bench.jl b/test/bench.jl index 93347449..c0f12f99 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -6,7 +6,7 @@ end using BenchmarkTools using StaticArrays using VortexStepMethod -using VortexStepMethod: calculate_AIC_matrices!, gamma_loop!, calculate_results, +using VortexStepMethod: calculate_AIC_matrices!, gamma_loop!, update_effective_angle_of_attack!, calculate_projected_area, calculate_cl, calculate_cd_cm, calculate_velocity_induced_single_ring_semiinfinite!, @@ -163,58 +163,7 @@ using LinearAlgebra end end - @testset "Results Calculation" begin - # Pre-allocate arrays - alpha_dist = zeros(n_panels) - v_rel_dist = zeros(n_panels) - chord_dist = zeros(n_panels) - x_airf_dist = zeros(n_panels, 3) - y_airf_dist = zeros(n_panels, 3) - z_airf_dist = zeros(n_panels, 3) - va_vec_dist = zeros(n_panels, 3) - va_dist = zeros(n_panels) - va_unit_dist = zeros(n_panels, 3) - reference_point = zeros(3) - - - set_va!(body_aero, va_vec) - # Fill arrays with panel data to satisfy calculate_results preconditions. - for (i, panel) in enumerate(body_aero.panels) - chord_dist[i] = panel.chord - x_airf_dist[i, :] .= panel.x_airf - y_airf_dist[i, :] .= panel.y_airf - z_airf_dist[i, :] .= panel.z_airf - va_vec_dist[i, :] .= panel.va_vec - va_dist[i] = norm(panel.va_vec) - va_unit_dist[i, :] .= - va_dist[i] > 0.0 ? panel.va_vec ./ va_dist[i] : [1.0, 0.0, 0.0] - v_rel_dist[i] = va_dist[i] - end - results = @MVector zeros(3) - - result = @benchmark calculate_results( - $body_aero, - $gamma, - $reference_point, - $density, - 1e-20, - 0.0, - $alpha_dist, - $v_rel_dist, - $chord_dist, - $x_airf_dist, - $z_airf_dist, - $va_vec_dist, - $va_dist, - $va_unit_dist, - $body_aero.panels, - false - ) samples=1 evals=1 - @info "Calculate Results Allocations: $(result.allocs) Memory: $(result.memory)" - @test result.allocs ≤ 700 - end - - @testset "Allocation Tests for solve() and solve!()" begin + @testset "Allocation Tests for solve_base!() and solve!()" begin result = @benchmark solve_base!($solver, $body_aero, nothing) samples=1 evals=1 @test result.allocs <= 55 # time Python: 32.0 ms Ryzen 7950x diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl index 893853c4..3183011e 100644 --- a/test/body_aerodynamics/test_body_aerodynamics.jl +++ b/test/body_aerodynamics/test_body_aerodynamics.jl @@ -343,10 +343,7 @@ end atol=1e-8, rtol=1e-8 ) - results_NEW = solve(loop_solver, body_aero; reference_point=[0,1,0]) - # println(results_NEW) - - @test results_NEW isa Dict + results_NEW = solve!(loop_solver, body_aero; reference_point=[0,1,0]) @testset "Loop and nonlin solve!" begin loop_sol = solve!(loop_solver, body_aero; reference_point=[0,1,0]) @@ -375,7 +372,7 @@ end end # Calculate forces using uncorrected alpha - alpha = results_NEW["alpha_uncorrected"] + alpha = results_NEW.alpha_uncorrected dyn_visc = 0.5 * density * norm(va_vec)^2 n_panels = length(body_aero.panels) lift = zeros(n_panels) @@ -392,7 +389,7 @@ end Fmag = hcat(lift, drag, moment) # Calculate coefficients using corrected alpha - alpha = results_NEW["alpha_at_ac"] + alpha = results_NEW.alpha_dist aero_coeffs = hcat( [alpha[i] for (i, panel) in enumerate(body_aero.panels)], [calculate_cl(panel, alpha[i]) for (i, panel) in enumerate(body_aero.panels)], @@ -410,25 +407,25 @@ end # Compare results @info "Comparing results" - @info "cl_calculated: $(results_NEW["cl"]), CL_ref: $CL_ref" - @info "cd_calculated: $(results_NEW["cd"]), CD_ref: $CD_ref" - @info "cs_calculated: $(results_NEW["cs"]), CS_ref: $CS_ref" - @info "L_calculated: $(results_NEW["lift"]), Ltot_ref: $Ltot_ref" - @info "D_calculated: $(results_NEW["drag"]), Dtot_ref: $Dtot_ref" + @info "cl_calculated: $(results_NEW.cl), CL_ref: $CL_ref" + @info "cd_calculated: $(results_NEW.cd), CD_ref: $CD_ref" + @info "cs_calculated: $(results_NEW.cs), CS_ref: $CS_ref" + @info "L_calculated: $(results_NEW.lift), Ltot_ref: $Ltot_ref" + @info "D_calculated: $(results_NEW.drag), Dtot_ref: $Dtot_ref" # Assert results - @test isapprox(results_NEW["cl"], CL_ref, rtol=1e-4) - @test isapprox(results_NEW["cd"], CD_ref, rtol=1e-4) - @test isapprox(results_NEW["cs"], CS_ref, rtol=1e-4) - @test isapprox(results_NEW["lift"], Ltot_ref, rtol=1e-4) - @test isapprox(results_NEW["drag"], Dtot_ref, rtol=1e-4) - @test isapprox(results_NEW["Fx"], results_NEW["Mz"], rtol=1e-4) # 1 meter arm - @test isapprox(results_NEW["My"], 0.0, atol=1e-3) - @test isapprox(results_NEW["Fz"], -results_NEW["Mx"], rtol=1e-4) # 1 meter arm + @test isapprox(results_NEW.cl, CL_ref, rtol=1e-4) + @test isapprox(results_NEW.cd, CD_ref, rtol=1e-4) + @test isapprox(results_NEW.cs, CS_ref, rtol=1e-4) + @test isapprox(results_NEW.lift, Ltot_ref, rtol=1e-4) + @test isapprox(results_NEW.drag, Dtot_ref, rtol=1e-4) + @test isapprox(results_NEW.force[1], results_NEW.moment[3], rtol=1e-4) # 1 meter arm + @test isapprox(results_NEW.moment[2], 0.0, atol=1e-3) + @test isapprox(results_NEW.force[3], -results_NEW.moment[1], rtol=1e-4) # 1 meter arm # Check array shapes - @test length(results_NEW["cl_distribution"]) == length(body_aero.panels) - @test length(results_NEW["cd_distribution"]) == length(body_aero.panels) + @test length(results_NEW.cl_distribution) == length(body_aero.panels) + @test length(results_NEW.cd_distribution) == length(body_aero.panels) end @testset "set_va! with VSMSettings applies the yaw rate about body z" begin diff --git a/test/plotting/test_plotting.jl b/test/plotting/test_plotting.jl index 6c780c60..0eb8f131 100644 --- a/test/plotting/test_plotting.jl +++ b/test/plotting/test_plotting.jl @@ -85,8 +85,8 @@ end aerodynamic_model_type=LLT) # Solve the VSM and LLT - results_vsm = solve(vsm_solver, body_aero) - results_llt = solve(llt_solver, body_aero) + results_vsm = solve!(vsm_solver, body_aero) + results_llt = solve!(llt_solver, body_aero) # Plot spanwise distributions y_coordinates = [panel.aero_center[2] diff --git a/test/solver/test_attached_trailed_force.jl b/test/solver/test_attached_trailed_force.jl index 610262b7..77c820a8 100644 --- a/test/solver/test_attached_trailed_force.jl +++ b/test/solver/test_attached_trailed_force.jl @@ -37,8 +37,8 @@ function lift_and_oswald(wing, aspect_ratio, model, is_with_attached_trailed_for solver = Solver(wing.n_panels, wing.n_unrefined_sections; use_gamma_prev=false, aerodynamic_model_type=model, is_with_attached_trailed_force) set_va!(body_aero, 20.0 .* [cosd(10), 0.0, sind(10)]) - results = solve(solver, body_aero) - return results["cl"], results["cl"]^2 / (π * aspect_ratio * results["cd"]) + sol = solve!(solver, body_aero) + return sol.cl, sol.cl^2 / (π * aspect_ratio * sol.cd) end @testset "Attached trailed vortex force" begin @@ -84,9 +84,6 @@ end @test solver_on.sol.f_body_3D[:, i] ≈ force_off[:, i] .+ attached.force @test solver_on.sol.m_body_3D[:, i] ≈ moment_off[:, i] .+ attached.moment end - results = solve(solver_on, body_aero) - @test results["F_distribution"] ≈ solver_on.sol.f_body_3D - @test results["M_distribution"] ≈ solver_on.sol.m_body_3D @test (@allocated calc_forces!(solver_on, body_aero)) == 0 end diff --git a/test/solver/test_moment_units.jl b/test/solver/test_moment_units.jl index 62bf2f9b..d03001ed 100644 --- a/test/solver/test_moment_units.jl +++ b/test/solver/test_moment_units.jl @@ -48,18 +48,4 @@ end @test large.moment_coeff_dist ≈ small.moment_coeff_dist rtol = 1e-6 end - @testset "solve" begin - small, large = map((1.0, k)) do scale - body_aero = scaled_wing_aero(scale) - wing = only(body_aero.wings) - solver = Solver(wing.n_panels, wing.n_unrefined_sections) - solve(solver, body_aero; reference_point=reference_point(scale)) - end - for key in ("Mx", "My", "Mz", "M_distribution") - @test large[key] ≈ k^3 .* small[key] rtol = 1e-6 - end - for key in ("cmx", "cmy", "cmz") - @test large[key] ≈ small[key] rtol = 1e-6 - end - end end diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index 089b3353..d0781491 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -60,7 +60,6 @@ end Solver(wing.n_panels, n_sections + 1)) @test_throws DimensionMismatch solve!(other, body_aero) @test_throws "Solver built for" solve!(other, body_aero) - @test_throws DimensionMismatch solve(other, body_aero) end end finally @@ -198,6 +197,45 @@ calc_forces_allocs(solver, body_aero) = end end +""" + flat_plate_aero(y_ranges; chord=1.0) + +`BodyAerodynamics` of one flat inviscid rectangle per `(y_start, y_end)` in `y_ranges`, +with its leading edge on `x = 0`, in a 5° inflow of 12 m/s. +""" +function flat_plate_aero(y_ranges; chord=1.0) + wings = map(y_ranges) do (y_start, y_end) + wing = Wing(10) + add_section!(wing, [0.0, y_start, 0.0], [chord, y_start, 0.0], INVISCID) + add_section!(wing, [0.0, y_end, 0.0], [chord, y_end, 0.0], INVISCID) + refine!(wing) + wing + end + body_aero = BodyAerodynamics(collect(wings)) + set_va!(body_aero, 12.0 .* [cosd(5), 0.0, sind(5)]) + return body_aero +end + +@testset "solve! fills the reference inflow and the centers of pressure" begin + body_aero = flat_plate_aero([(4.0, -4.0)]) + solver = Solver(length(body_aero.panels), 2) + sol = solve!(solver, body_aero) + @test sol.va_ref_vec ≈ 12.0 .* [cosd(5), 0.0, sind(5)] + @test sol.q_ref ≈ 0.5 * solver.density * 12.0^2 + @test sol.rey ≈ solver.density * 12.0 * body_aero.c_ref / solver.mu + @test sol.alpha_uncorrected == solver.lr.alpha_dist + # a flat plate carries no section moment, so each load acts at its quarter chord + @test sol.center_of_pressure ≈ [0.25, 0.0, 0.0] atol = 1e-6 + for (location, panel) in zip(sol.panel_cp_locations, body_aero.panels) + @test location ≈ panel.aero_center + end + + # two plates with a gap between them: the line of action runs through the gap + gapped = flat_plate_aero([(6.0, 2.0), (-2.0, -6.0)]) + solver = Solver(length(gapped.panels), 4) + @test all(isnan, solve!(solver, gapped).center_of_pressure) +end + @testset "Spanwise Laplacian tip closures" begin # Interior three-point stencil plus the Eq. 15 tip closures. n = 5 diff --git a/test/solver/test_viscous_drag_correction.jl b/test/solver/test_viscous_drag_correction.jl index 628a87c4..06f1e560 100644 --- a/test/solver/test_viscous_drag_correction.jl +++ b/test/solver/test_viscous_drag_correction.jl @@ -71,14 +71,6 @@ end end end - @testset "solve reports the corrected forces" begin - set_va!(body_aero, va_vec_sideslip) - solve!(solver_on, body_aero) - results = solve(solver_on, body_aero) - @test [results["Fx"], results["Fy"], results["Fz"]] ≈ solver_on.sol.force - @test results["F_distribution"] ≈ solver_on.sol.f_body_3D - end - @testset "linearize reports the corrected forces" begin y = [va_vec_sideslip; zeros(3)] results_for(solver) = VortexStepMethod.linearize(solver, body_aero, y; diff --git a/test/solver/test_wing_directions.jl b/test/solver/test_wing_directions.jl index e1ba4e61..0926de00 100644 --- a/test/solver/test_wing_directions.jl +++ b/test/solver/test_wing_directions.jl @@ -44,22 +44,13 @@ end @test pair.f_body_3D[:, rotated] ≈ rotation * solo.f_body_3D rtol = 1e-5 @test pair.lift_dist[rotated] ≈ solo.lift_dist rtol = 1e-5 @test pair.drag_dist[rotated] ≈ solo.drag_dist rtol = 1e-5 - end - - @testset "solve" begin - solo, pair = map((solo_aero, pair_aero)) do body_aero - solver = Solver(length(body_aero.panels), 2length(body_aero.wings)) - solve(solver, body_aero) - end - rotated_forces = pair["F_distribution"][:, rotated] - @test rotated_forces ≈ rotation * solo["F_distribution"] rtol = 1e-5 - for key in ("cl_distribution", "cd_distribution", "cs_distribution") - @test pair[key][rotated] ≈ solo[key] rtol = 1e-5 atol = 1e-8 + for field in (:cl_distribution, :cd_distribution, :cs_distribution) + @test getfield(pair, field)[rotated] ≈ getfield(solo, field) rtol = 1e-5 atol = 1e-8 end # along x the inflow makes z the body's lift and y its side direction - @test pair["lift"] ≈ pair["Fz"] rtol = 1e-10 - @test pair["side"] ≈ pair["Fy"] rtol = 1e-10 - @test pair["drag"] ≈ pair["Fx"] rtol = 1e-10 + @test pair.lift ≈ pair.force[3] rtol = 1e-10 + @test pair.side ≈ pair.force[2] rtol = 1e-10 + @test pair.drag ≈ pair.force[1] rtol = 1e-10 end end @@ -68,7 +59,7 @@ end rectangular_wing(I(3), [0.0, 6.0, 0.0])]) set_va!(body_aero, [20.0, 0.0, 0.0]) solver = Solver(length(body_aero.panels), 2length(body_aero.wings)) - results = solve(solver, body_aero) - @test results["wing_span"] ≈ 12.0 - @test results["aspect_ratio_projected"] ≈ 12.0^2 / 18.0 + sol = solve!(solver, body_aero) + @test sol.wing_span ≈ 12.0 + @test sol.aspect_ratio_projected ≈ 12.0^2 / 18.0 end diff --git a/test/verification/test_verification.jl b/test/verification/test_verification.jl index 8dcff408..cdf39467 100644 --- a/test/verification/test_verification.jl +++ b/test/verification/test_verification.jl @@ -62,9 +62,9 @@ function lift_drag_polar(wing, model, alphas; va, relaxation_factor) CD = zeros(length(alphas)) for (i, alpha) in enumerate(alphas) set_va!(body_aero, va .* [cosd(alpha), 0.0, sind(alpha)]) - results = solve(solver, body_aero) - CL[i] = results["cl"] - CD[i] = results["cd"] + sol = solve!(solver, body_aero, nothing) + CL[i] = sol.cl + CD[i] = sol.cd end return CL, CD end From 5e56e6f5c57ef37794369ce65b955affcb96b223 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Thu, 24 Sep 2026 11:26:19 +0200 Subject: [PATCH 2/4] Let calc_only_f_and_gamma skip the analysis fields of VSMSolution The setting's only consumer was solve()'s short Dict; it now keeps solve! from filling the projections, span and centers of pressure. Adds the changelog entries and wraps the lines over 92 characters. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 28 ++++++++++++++ docs/src/private_functions.md | 2 +- ext/VortexStepMethodMakieExt.jl | 3 +- src/VortexStepMethod.jl | 4 +- src/settings.jl | 5 ++- src/solver.jl | 59 +++++++++++++++++------------ test/solver/test_solver.jl | 7 ++++ test/solver/test_wing_directions.jl | 3 +- 8 files changed, 79 insertions(+), 32 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 50b26ddd..e07c5df6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,33 @@ # Changelog +## Unreleased + +### Added + +- `VSMSolution` holds everything the removed `solve` dictionary held, filled by every + `solve!`: `lift`, `drag`, `side`, `cl`, `cd`, `cs`, `cl_distribution`, + `cd_distribution`, `cs_distribution`, `alpha_uncorrected`, `va_ref_vec`, `q_ref`, `rey`, + `area_all_panels`, `projected_area`, `wing_span`, `aspect_ratio_projected`, + `center_of_pressure` and `panel_cp_locations`. `calc_only_f_and_gamma` skips the + projections, span and centers of pressure. + +### Changed + +- BREAKING: `solve` and `calculate_results` are removed; use `solve!`, which returns the + solver's `VSMSolution`, and `solve!(solver, body_aero, nothing)` to start from a fresh + circulation as `solve` did. Keys that became other fields: `Fx`/`Fy`/`Fz` → `force`, + `Mx`/`My`/`Mz` → `moment`, `cfx`… → `force_coeffs`, `cmx`… → `moment_coeffs`, + `F_distribution` → `f_body_3D`, `M_distribution` → `m_body_3D`, `alpha_at_ac` → + `alpha_dist`, `alpha_geometric` → `alpha_geometric_dist`, `Rey` → `rey`. +- BREAKING: `plot_distribution` and `plot_combined_analysis` take `VSMSolution`s instead of + dictionaries. `plot_combined_analysis`, `plot_polars` and `generate_polar_data` solve + through `solve!`, so they leave `solver.sol` at their last angle. +- An LLT solver with `correct_aoa` no longer corrects the angle of attack in any result, + as `solve!` already did; `solve` did, and gave `cl` up to 0.3 % and `cmx` up to 4.5 % + apart from `solve!`. +- `center_of_pressure` is `NaN` where the line of action crosses no panel, instead of + `nothing` with a warning. + ## VortexStepMethod v6.0.0 2026-09-23 ### Added diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index a7912813..b9da3b72 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -64,7 +64,7 @@ flow_curvature_cm spanwise_flow_drag panel_force_directions prescribed_va_directions -inflow_loads! +analysis_fields! find_center_of_pressure! panel_moment panel_couple_force diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index a1c6306c..e69ac82b 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/ext/VortexStepMethodMakieExt.jl @@ -831,7 +831,8 @@ function VortexStepMethod.plot_distribution(y_coordinates_list, results_list, la # Plot alpha geometric for (y_coords, results, label) in zip(y_coordinates_list, results_list, label_list) - lines!(ax_alpha_geo, Vector(y_coords), rad2deg.(Vector(results.alpha_geometric_dist)), + lines!(ax_alpha_geo, Vector(y_coords), + rad2deg.(Vector(results.alpha_geometric_dist)), label=label) end diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index f0b085ac..b29870d5 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -86,7 +86,7 @@ Plot spanwise distributions of aerodynamic properties. # Arguments - `y_coordinates_list`: list of spanwise coordinate arrays -- `results_list`: list of result dictionaries from [`solve!`](@ref) +- `results_list`: list of [`VSMSolution`](@ref)s from [`solve!`](@ref) - `label_list`: list of labels for each result # Keyword arguments @@ -151,7 +151,7 @@ in sequence. # Arguments - `solver`: solver or vector of solvers - `body_aero`: [`BodyAerodynamics`](@ref) object or vector thereof -- `results`: results dictionary (or vector) from [`solve!`](@ref) +- `results`: [`VSMSolution`](@ref) (or vector of them) from [`solve!`](@ref) # Keyword arguments - `solver_label`: label string for the solver (backward-compatible alias for `labels`) diff --git a/src/settings.jl b/src/settings.jl index 354d76c9..4a97a54c 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -171,7 +171,8 @@ Solver configuration, used within [`VSMSettings`](@ref). - `core_radius_fraction`: Bound vortex core cut-off, as a fraction of the filament length, following Damiani et al. (2019) (default `0.05`) - `mu`: Dynamic viscosity (N*s/m^2) (default `1.81e-5`) -- `calc_only_f_and_gamma`: Only output forces and circulation +- `calc_only_f_and_gamma`: Leave the analysis fields of [`VSMSolution`](@ref) (`lift`, + `cl`, their distributions, span and centers of pressure) at their last values (default `false`) - `correct_aoa`: Perform angle of attack correction (default `false`) @@ -197,7 +198,7 @@ Solver configuration, used within [`VSMSettings`](@ref). use_gamma_prev::Bool = true # if false, always reinitialize gamma from type_initial_gamma_distribution core_radius_fraction::Float64 = 0.05 mu::Float64 = 1.81e-5 # dynamic viscosity [N·s/m²] - calc_only_f_and_gamma::Bool=false # whether to only output f and gamma + calc_only_f_and_gamma::Bool=false # skip the analysis fields of VSMSolution correct_aoa::Bool=false # perform aoa correction flow_curvature::Bool=false # thin-airfoil pitch-rate moment increment is_with_viscous_drag_correction::Bool=false # spanwise-flow viscous force (Gaunaa 2024) diff --git a/src/solver.jl b/src/solver.jl index 7779e616..6df4627d 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -24,12 +24,13 @@ Struct for storing the solution of the [`solve!`](@ref) function. Must contain a - moment::MVec3: Aerodynamic moments [Mx, My, Mz] around the reference point [Nm] - force_coeffs::MVec3: Aerodynamic force coefficients [CFx, CFy, CFz] [-] - `moment_coeffs`::MVec3: Aerodynamic moment coefficients [CMx, CMy, CMz] [-] -- `lift`, `drag`, `side`: Total force along the lift, drag and side directions of the +- `lift`, `drag`, `side`¹: Total force along the lift, drag and side directions of the reference inflow `va_ref_vec` [N] -- `cl`, `cd`, `cs`: `lift`, `drag` and `side` divided by `q_ref * projected_area` [-] -- `cl_distribution`, `cd_distribution`, `cs_distribution`::Vector{Float64}: Panel force along - the lift, drag and side directions of the panel's own inflow, divided by its dynamic - pressure and area [-]; unlike `cl_dist` and `cd_dist`, these include every force term +- `cl`, `cd`, `cs`¹: `lift`, `drag` and `side` divided by `q_ref * projected_area` [-] +- `cl_distribution`, `cd_distribution`, `cs_distribution`¹::Vector{Float64}: Panel force + along the lift, drag and side directions of the panel's own inflow, divided by its + dynamic pressure and area [-]; unlike `cl_dist` and `cd_dist`, these include every + force term - `alpha_uncorrected`::Vector{Float64}: Angle of attack of each panel at its evaluation point, before the aerodynamic-center correction [rad] - `va_ref_vec`::MVec3: Area-weighted reference inflow velocity [m/s] @@ -37,11 +38,11 @@ Struct for storing the solution of the [`solve!`](@ref) function. Must contain a - `rey`: Reynolds number of `va_ref_vec` on the reference chord [-] - `area_all_panels`: Sum of the panel areas [m²] - `projected_area`: Projected area of the body [m²] -- `wing_span`: Span along the first wing's spanwise direction [m] -- `aspect_ratio_projected`: `wing_span^2 / projected_area` [-] -- `center_of_pressure`::MVec3: Point where the line of action of `force` crosses a +- `wing_span`¹: Span along the first wing's spanwise direction [m] +- `aspect_ratio_projected`¹: `wing_span^2 / projected_area` [-] +- `center_of_pressure`¹::MVec3: Point where the line of action of `force` crosses a panel, `NaN` where it crosses none [m] -- `panel_cp_locations`::Vector{MVec3}: Center of pressure of each panel [m] +- `panel_cp_locations`¹::Vector{MVec3}: Center of pressure of each panel [m] - `moment_dist`::Vector{Float64}: Pitching moments around the spanwise vector of each panel. [Nm] - `moment_coeff_dist`::Vector{Float64}: Pitching moment coefficient around the spanwise vector of each panel. [-] - `moment_unrefined_dist`::MVector{U, Float64}: Averaged moments for unrefined sections [Nm] @@ -51,6 +52,8 @@ Struct for storing the solution of the [`solve!`](@ref) function. Must contain a - `moment_coeff_unrefined_dist`::MVector{U, Float64}: Summed `moment_frac`-referenced pitching-moment coefficient per unrefined section [-] - `alpha_unrefined_dist`::MVector{U, Float64}: Averaged angles of attack for unrefined sections [rad] - `solver_status`::SolverStatus: enum, see [`SolverStatus`](@ref) + +¹ Kept at its last value when the solver's `calc_only_f_and_gamma` is set. """ @with_kw mutable struct VSMSolution{P, U, T} ### private vectors of solve_base! @@ -167,7 +170,8 @@ Main solver structure for the Vortex Step Method. See also: [`solve!`](@ref) - `core_radius_fraction`::Float64 = 0.05: Bound vortex core cut-off, as a fraction of the filament length, following Damiani et al. (2019) - mu::Float64 = 1.81e-5: Dynamic viscosity [N·s/m²] -- `is_only_f_and_gamma_output`::Bool = false: Whether to only output f and gamma +- `is_only_f_and_gamma_output`::Bool = false: Whether `solve!` skips the fields of + [`VSMSolution`](@ref) that only analysis reads, see [`SolverSettings`](@ref) - `flow_curvature`::Bool = false: Add the thin-airfoil pitch-rate moment increment `-(π/4) q̂` to each section, see: [`flow_curvature_cm`](@ref) - `is_with_viscous_drag_correction`::Bool = false: Add the spanwise-flow viscous drag @@ -619,14 +623,8 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; solver.sol.alpha_uncorrected .= alpha_dist solver.sol.area_all_panels = area_all_panels solver.sol.projected_area = projected_area - reference_spanwise = body_aero.wings[1].spanwise_direction - solver.sol.wing_span = calculate_span(body_aero.wings, reference_spanwise) - solver.sol.aspect_ratio_projected = solver.sol.wing_span^2 / projected_area - inflow_loads!(solver.sol, body_aero, density) - find_center_of_pressure!(solver.sol.center_of_pressure, body_aero, solver.sol.force, - solver.sol.moment, reference_point) - compute_panel_center_of_pressures!(solver.sol.panel_cp_locations, body_aero, - solver.sol.f_body_3D, solver.sol.m_body_3D, reference_point) + solver.is_only_f_and_gamma_output || + analysis_fields!(solver.sol, body_aero, density, reference_point) if converged # TODO: Check if the result if feasible if converged solver.sol.solver_status = FEASIBLE @@ -638,14 +636,24 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; end """ - inflow_loads!(sol::VSMSolution, body_aero::BodyAerodynamics, density) - -Fill `lift`, `drag`, `side`, `cl`, `cd`, `cs` and `cl_distribution`, `cd_distribution`, -`cs_distribution` of `sol` by projecting the panel forces `sol.f_body_3D` on the lift, drag -and side directions of the reference inflow `sol.va_ref_vec` and of each panel's own inflow. + analysis_fields!(sol::VSMSolution, body_aero::BodyAerodynamics, density, + reference_point) + +Fill the fields of `sol` that `calc_only_f_and_gamma` skips: `wing_span`, +`aspect_ratio_projected`, both centers of pressure, and `lift`, `drag`, `side`, `cl`, `cd`, +`cs` and their `_distribution`s, which project the panel forces `sol.f_body_3D` on the +lift, drag and side directions of the reference inflow `sol.va_ref_vec` and of each +panel's own inflow. """ -function inflow_loads!(sol::VSMSolution, body_aero::BodyAerodynamics, density) +function analysis_fields!(sol::VSMSolution, body_aero::BodyAerodynamics, density, + reference_point) reference_spanwise = SVector{3}(body_aero.wings[1].spanwise_direction) + sol.wing_span = calculate_span(body_aero.wings, reference_spanwise) + sol.aspect_ratio_projected = sol.wing_span^2 / sol.projected_area + find_center_of_pressure!(sol.center_of_pressure, body_aero, sol.force, sol.moment, + reference_point) + compute_panel_center_of_pressures!(sol.panel_cp_locations, body_aero, sol.f_body_3D, + sol.m_body_3D, reference_point) va_ref = SVector{3}(sol.va_ref_vec) va_ref_unit = va_ref / norm(va_ref) reference_dirs = prescribed_va_directions(va_ref, reference_spanwise) @@ -653,7 +661,8 @@ function inflow_loads!(sol::VSMSolution, body_aero::BodyAerodynamics, density) for (wing_idx, wing) in enumerate(body_aero.wings) spanwise_unit = SVector{3}(wing.spanwise_direction) for i in panel_range(body_aero, wing_idx) - force = SVector{3}(sol.f_body_3D[1, i], sol.f_body_3D[2, i], sol.f_body_3D[3, i]) + force = SVector{3}(sol.f_body_3D[1, i], sol.f_body_3D[2, i], + sol.f_body_3D[3, i]) panel_va = SVector{3}(sol.va_vec_dist[i, 1], sol.va_vec_dist[i, 2], sol.va_vec_dist[i, 3]) va_panel = norm(panel_va) diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index d0781491..55ec6db8 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -230,6 +230,13 @@ end @test location ≈ panel.aero_center end + # calc_only_f_and_gamma leaves the analysis fields alone + skipping = Solver(length(body_aero.panels), 2; is_only_f_and_gamma_output=true) + skipped = solve!(skipping, body_aero) + @test skipped.force ≈ sol.force + @test skipped.cl == 0.0 + @test all(iszero, skipped.cl_distribution) + # two plates with a gap between them: the line of action runs through the gap gapped = flat_plate_aero([(6.0, 2.0), (-2.0, -6.0)]) solver = Solver(length(gapped.panels), 4) diff --git a/test/solver/test_wing_directions.jl b/test/solver/test_wing_directions.jl index 0926de00..47314192 100644 --- a/test/solver/test_wing_directions.jl +++ b/test/solver/test_wing_directions.jl @@ -45,7 +45,8 @@ end @test pair.lift_dist[rotated] ≈ solo.lift_dist rtol = 1e-5 @test pair.drag_dist[rotated] ≈ solo.drag_dist rtol = 1e-5 for field in (:cl_distribution, :cd_distribution, :cs_distribution) - @test getfield(pair, field)[rotated] ≈ getfield(solo, field) rtol = 1e-5 atol = 1e-8 + rotated_values = getfield(pair, field)[rotated] + @test rotated_values ≈ getfield(solo, field) rtol = 1e-5 atol = 1e-8 end # along x the inflow makes z the body's lift and y its side direction @test pair.lift ≈ pair.force[3] rtol = 1e-10 From 8a35a0741deee63cde56dcfa4a3cb382ef481f78 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Thu, 24 Sep 2026 11:52:06 +0200 Subject: [PATCH 3/4] Pin the calc_only_f_and_gamma keep-last-value contract and tidy after review The skip test now flips the flag on a solved solver and asserts the analysis fields keep their values at a new inflow. Adds the compute_panel_center_of_pressures! docstring, says where cl_distribution is NaN and that correct_aoa is VSM-only, restates the CHANGELOG line, and solves the loop case in test_body_aerodynamics once. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 5 +-- docs/src/private_functions.md | 1 + src/body_aerodynamics.jl | 9 +++++ src/settings.jl | 2 +- src/solver.jl | 4 +- .../test_body_aerodynamics.jl | 37 +++++++++---------- test/solver/test_moment_units.jl | 1 - test/solver/test_solver.jl | 23 +++++++----- 8 files changed, 46 insertions(+), 36 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e07c5df6..ee2ef710 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,9 +22,8 @@ - BREAKING: `plot_distribution` and `plot_combined_analysis` take `VSMSolution`s instead of dictionaries. `plot_combined_analysis`, `plot_polars` and `generate_polar_data` solve through `solve!`, so they leave `solver.sol` at their last angle. -- An LLT solver with `correct_aoa` no longer corrects the angle of attack in any result, - as `solve!` already did; `solve` did, and gave `cl` up to 0.3 % and `cmx` up to 4.5 % - apart from `solve!`. +- `correct_aoa` applies to the VSM model only. With LLT, results that came from `solve` + move by up to 0.34 % in `cl`, 0.43 % in `cd` and 4.5 % in `cmx`. - `center_of_pressure` is `NaN` where the line of action crosses no panel, instead of `nothing` with a warning. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index b9da3b72..8905a7ac 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -66,6 +66,7 @@ panel_force_directions prescribed_va_directions analysis_fields! find_center_of_pressure! +compute_panel_center_of_pressures! panel_moment panel_couple_force panel_loads diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 0e73c49a..7f8dea82 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -788,6 +788,15 @@ function find_center_of_pressure!( return center_of_pressure end +""" + compute_panel_center_of_pressures!(panel_cp_locations, body_aero::BodyAerodynamics, + f_distribution, m_distribution, reference_point) + +Set each entry of `panel_cp_locations` to the point on its panel's chord, clamped between +leading and trailing edge, where the panel's column of `f_distribution` gives its column of +`m_distribution` about `reference_point`; to the aerodynamic center where that force has +no finite component normal to the chord. +""" function compute_panel_center_of_pressures!( panel_cp_locations, body_aero::BodyAerodynamics, diff --git a/src/settings.jl b/src/settings.jl index 4a97a54c..4062d56d 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -174,7 +174,7 @@ Solver configuration, used within [`VSMSettings`](@ref). - `calc_only_f_and_gamma`: Leave the analysis fields of [`VSMSolution`](@ref) (`lift`, `cl`, their distributions, span and centers of pressure) at their last values (default `false`) -- `correct_aoa`: Perform angle of attack correction +- `correct_aoa`: Perform angle of attack correction (VSM model only) (default `false`) - `flow_curvature`: Add the thin-airfoil pitch-rate moment increment to each section (default `false`) diff --git a/src/solver.jl b/src/solver.jl index 6df4627d..6af0da23 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -29,8 +29,8 @@ Struct for storing the solution of the [`solve!`](@ref) function. Must contain a - `cl`, `cd`, `cs`¹: `lift`, `drag` and `side` divided by `q_ref * projected_area` [-] - `cl_distribution`, `cd_distribution`, `cs_distribution`¹::Vector{Float64}: Panel force along the lift, drag and side directions of the panel's own inflow, divided by its - dynamic pressure and area [-]; unlike `cl_dist` and `cd_dist`, these include every - force term + dynamic pressure and area [-], `NaN` on a panel without inflow; unlike `cl_dist` and + `cd_dist`, these include every force term - `alpha_uncorrected`::Vector{Float64}: Angle of attack of each panel at its evaluation point, before the aerodynamic-center correction [rad] - `va_ref_vec`::MVec3: Area-weighted reference inflow velocity [m/s] diff --git a/test/body_aerodynamics/test_body_aerodynamics.jl b/test/body_aerodynamics/test_body_aerodynamics.jl index 3183011e..a47650a6 100644 --- a/test/body_aerodynamics/test_body_aerodynamics.jl +++ b/test/body_aerodynamics/test_body_aerodynamics.jl @@ -343,10 +343,9 @@ end atol=1e-8, rtol=1e-8 ) - results_NEW = solve!(loop_solver, body_aero; reference_point=[0,1,0]) + loop_sol = solve!(loop_solver, body_aero; reference_point=[0,1,0]) @testset "Loop and nonlin solve!" begin - loop_sol = solve!(loop_solver, body_aero; reference_point=[0,1,0]) nonlin_sol = solve!(nonlin_solver, body_aero; reference_point=[0,1,0]) @test all(isapprox.(nonlin_sol.gamma_distribution, loop_sol.gamma_distribution; atol=1e-4)) @@ -372,7 +371,7 @@ end end # Calculate forces using uncorrected alpha - alpha = results_NEW.alpha_uncorrected + alpha = loop_sol.alpha_uncorrected dyn_visc = 0.5 * density * norm(va_vec)^2 n_panels = length(body_aero.panels) lift = zeros(n_panels) @@ -389,7 +388,7 @@ end Fmag = hcat(lift, drag, moment) # Calculate coefficients using corrected alpha - alpha = results_NEW.alpha_dist + alpha = loop_sol.alpha_dist aero_coeffs = hcat( [alpha[i] for (i, panel) in enumerate(body_aero.panels)], [calculate_cl(panel, alpha[i]) for (i, panel) in enumerate(body_aero.panels)], @@ -407,25 +406,25 @@ end # Compare results @info "Comparing results" - @info "cl_calculated: $(results_NEW.cl), CL_ref: $CL_ref" - @info "cd_calculated: $(results_NEW.cd), CD_ref: $CD_ref" - @info "cs_calculated: $(results_NEW.cs), CS_ref: $CS_ref" - @info "L_calculated: $(results_NEW.lift), Ltot_ref: $Ltot_ref" - @info "D_calculated: $(results_NEW.drag), Dtot_ref: $Dtot_ref" + @info "cl_calculated: $(loop_sol.cl), CL_ref: $CL_ref" + @info "cd_calculated: $(loop_sol.cd), CD_ref: $CD_ref" + @info "cs_calculated: $(loop_sol.cs), CS_ref: $CS_ref" + @info "L_calculated: $(loop_sol.lift), Ltot_ref: $Ltot_ref" + @info "D_calculated: $(loop_sol.drag), Dtot_ref: $Dtot_ref" # Assert results - @test isapprox(results_NEW.cl, CL_ref, rtol=1e-4) - @test isapprox(results_NEW.cd, CD_ref, rtol=1e-4) - @test isapprox(results_NEW.cs, CS_ref, rtol=1e-4) - @test isapprox(results_NEW.lift, Ltot_ref, rtol=1e-4) - @test isapprox(results_NEW.drag, Dtot_ref, rtol=1e-4) - @test isapprox(results_NEW.force[1], results_NEW.moment[3], rtol=1e-4) # 1 meter arm - @test isapprox(results_NEW.moment[2], 0.0, atol=1e-3) - @test isapprox(results_NEW.force[3], -results_NEW.moment[1], rtol=1e-4) # 1 meter arm + @test isapprox(loop_sol.cl, CL_ref, rtol=1e-4) + @test isapprox(loop_sol.cd, CD_ref, rtol=1e-4) + @test isapprox(loop_sol.cs, CS_ref, rtol=1e-4) + @test isapprox(loop_sol.lift, Ltot_ref, rtol=1e-4) + @test isapprox(loop_sol.drag, Dtot_ref, rtol=1e-4) + @test isapprox(loop_sol.force[1], loop_sol.moment[3], rtol=1e-4) # 1 meter arm + @test isapprox(loop_sol.moment[2], 0.0, atol=1e-3) + @test isapprox(loop_sol.force[3], -loop_sol.moment[1], rtol=1e-4) # 1 meter arm # Check array shapes - @test length(results_NEW.cl_distribution) == length(body_aero.panels) - @test length(results_NEW.cd_distribution) == length(body_aero.panels) + @test length(loop_sol.cl_distribution) == length(body_aero.panels) + @test length(loop_sol.cd_distribution) == length(body_aero.panels) end @testset "set_va! with VSMSettings applies the yaw rate about body z" begin diff --git a/test/solver/test_moment_units.jl b/test/solver/test_moment_units.jl index d03001ed..a6f69af2 100644 --- a/test/solver/test_moment_units.jl +++ b/test/solver/test_moment_units.jl @@ -47,5 +47,4 @@ end @test large.moment_coeffs ≈ small.moment_coeffs rtol = 1e-6 @test large.moment_coeff_dist ≈ small.moment_coeff_dist rtol = 1e-6 end - end diff --git a/test/solver/test_solver.jl b/test/solver/test_solver.jl index 55ec6db8..9eedcd37 100644 --- a/test/solver/test_solver.jl +++ b/test/solver/test_solver.jl @@ -198,12 +198,12 @@ calc_forces_allocs(solver, body_aero) = end """ - flat_plate_aero(y_ranges; chord=1.0) + inviscid_plates_aero(y_ranges; chord=1.0) `BodyAerodynamics` of one flat inviscid rectangle per `(y_start, y_end)` in `y_ranges`, with its leading edge on `x = 0`, in a 5° inflow of 12 m/s. """ -function flat_plate_aero(y_ranges; chord=1.0) +function inviscid_plates_aero(y_ranges; chord=1.0) wings = map(y_ranges) do (y_start, y_end) wing = Wing(10) add_section!(wing, [0.0, y_start, 0.0], [chord, y_start, 0.0], INVISCID) @@ -217,7 +217,7 @@ function flat_plate_aero(y_ranges; chord=1.0) end @testset "solve! fills the reference inflow and the centers of pressure" begin - body_aero = flat_plate_aero([(4.0, -4.0)]) + body_aero = inviscid_plates_aero([(4.0, -4.0)]) solver = Solver(length(body_aero.panels), 2) sol = solve!(solver, body_aero) @test sol.va_ref_vec ≈ 12.0 .* [cosd(5), 0.0, sind(5)] @@ -230,15 +230,18 @@ end @test location ≈ panel.aero_center end - # calc_only_f_and_gamma leaves the analysis fields alone - skipping = Solver(length(body_aero.panels), 2; is_only_f_and_gamma_output=true) - skipped = solve!(skipping, body_aero) - @test skipped.force ≈ sol.force - @test skipped.cl == 0.0 - @test all(iszero, skipped.cl_distribution) + # calc_only_f_and_gamma keeps the analysis fields at their last value + cl, cl_distribution, lift = sol.cl, copy(sol.cl_distribution), sol.lift + solver.is_only_f_and_gamma_output = true + set_va!(body_aero, 12.0 .* [cosd(10), 0.0, sind(10)]) + skipped = solve!(solver, body_aero) + @test skipped.force[3] > 1.5 * lift + @test skipped.cl == cl + @test skipped.lift == lift + @test skipped.cl_distribution == cl_distribution # two plates with a gap between them: the line of action runs through the gap - gapped = flat_plate_aero([(6.0, 2.0), (-2.0, -6.0)]) + gapped = inviscid_plates_aero([(6.0, 2.0), (-2.0, -6.0)]) solver = Solver(length(gapped.panels), 4) @test all(isnan, solve!(solver, gapped).center_of_pressure) end From e9bd10c5f50df0d5549ad3f2985f71c3b5c69426 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 28 Sep 2026 18:38:30 +0200 Subject: [PATCH 4/4] Draw the spanwise plot before the polar sweep in V3_kite and pyramid_model MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit solve! returns the solver's own VSMSolution, and plot_polars re-solves the same solver over the whole angle range, so the kept `results` held the 25° solution by the time plot_distribution drew it under the configured angle. Co-Authored-By: Claude Opus 5.5 --- examples/V3_kite.jl | 38 +++++++++++++++++++------------------- examples/pyramid_model.jl | 34 +++++++++++++++++----------------- 2 files changed, 36 insertions(+), 36 deletions(-) diff --git a/examples/V3_kite.jl b/examples/V3_kite.jl index af01538e..46ed8772 100644 --- a/examples/V3_kite.jl +++ b/examples/V3_kite.jl @@ -78,25 +78,6 @@ yaw_rate = settings.condition.yaw_rate set_va!(body_aero, settings) results = solve!(solver, body_aero; log=true) -PLOT && plot_polars( - solvers, - bodies, - labels, - literature_path_list=literature_paths, - angle_range=range(-5, 25, length=31), - angle_type="angle_of_attack", - angle_of_attack=angle_of_attack_deg, - side_slip=sideslip_deg, - va=va, - title="$(wing.n_panels)_panels_$(wing.spanwise_distribution)_from_yaml_settings", - save_path=OUTPUT_DIR, - is_save=false || SAVE_ALL, - is_show=true, - use_tex=USE_TEX, - show_moments=false, - cl_over_cd=true -) - # Plotting geometry PLOT && plot_geometry( body_aero, @@ -124,6 +105,25 @@ PLOT && plot_distribution( use_tex=USE_TEX ) +PLOT && plot_polars( + solvers, + bodies, + labels, + literature_path_list=literature_paths, + angle_range=range(-5, 25, length=31), + angle_type="angle_of_attack", + angle_of_attack=angle_of_attack_deg, + side_slip=sideslip_deg, + va=va, + title="$(wing.n_panels)_panels_$(wing.spanwise_distribution)_from_yaml_settings", + save_path=OUTPUT_DIR, + is_save=false || SAVE_ALL, + is_show=true, + use_tex=USE_TEX, + show_moments=false, + cl_over_cd=true +) + # --- Beta sweep --- PLOT && plot_polars( solvers, diff --git a/examples/pyramid_model.jl b/examples/pyramid_model.jl index 96f8dd5f..e9389426 100644 --- a/examples/pyramid_model.jl +++ b/examples/pyramid_model.jl @@ -34,23 +34,6 @@ SAVE_ALL = false USE_TEX = false OUTPUT_DIR = joinpath(dirname(@__DIR__), "output") -# Plotting polars -PLOT && plot_polars( - [solver], - [body_aero], - ["VSM Pyramid Model"], - angle_range=range(-5, 25, length=30), - angle_type="angle_of_attack", - angle_of_attack=angle_of_attack_deg, - side_slip=sideslip_deg, - va=va, - title="$(wing.n_panels)_panels_$(wing.spanwise_distribution)_pyramid_model", - save_path=OUTPUT_DIR, - is_save=false || SAVE_ALL, - is_show=true, - use_tex=USE_TEX -) - # Plotting geometry PLOT && plot_geometry( body_aero, @@ -77,4 +60,21 @@ PLOT && plot_distribution( use_tex=USE_TEX ) +# Plotting polars +PLOT && plot_polars( + [solver], + [body_aero], + ["VSM Pyramid Model"], + angle_range=range(-5, 25, length=30), + angle_type="angle_of_attack", + angle_of_attack=angle_of_attack_deg, + side_slip=sideslip_deg, + va=va, + title="$(wing.n_panels)_panels_$(wing.spanwise_distribution)_pyramid_model", + save_path=OUTPUT_DIR, + is_save=false || SAVE_ALL, + is_show=true, + use_tex=USE_TEX +) + nothing \ No newline at end of file