From cadd8bc53ec31b5da5fda883ee244749dcdb9f8a Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 21 Sep 2026 14:09:41 +0200 Subject: [PATCH 1/4] Give each panel its own wing's spanwise direction in calc_forces! and calculate_results Co-Authored-By: Claude Opus 5 --- src/body_aerodynamics.jl | 146 +++++++++++++++------------- src/solver.jl | 78 +++++++-------- test/solver/test_wing_directions.jl | 59 +++++++++++ 3 files changed, 172 insertions(+), 111 deletions(-) create mode 100644 test/solver/test_wing_directions.jl diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index be9f9156..fa492024 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -270,6 +270,30 @@ function unrefined_section_range(body_aero::BodyAerodynamics, wing_idx) return offset .+ (1:body_aero.wings[wing_idx].n_unrefined_sections) end +""" + prescribed_va_directions(va, spanwise) + +Lift and side unit vectors, `(; dir_lift, dir_side)`, of inflow `va` on a wing along +`spanwise`: lift normal to both, side normal to lift and `va`. +""" +@inline function prescribed_va_directions(va, spanwise) + dir_lift = normalize(cross(va, spanwise)) + return (; dir_lift, dir_side=cross(dir_lift, va) / norm(va)) +end + +""" + panel_range(body_aero::BodyAerodynamics, wing_idx) + +Indices of the panels of wing `wing_idx` in `body_aero.panels`. +""" +function panel_range(body_aero::BodyAerodynamics, wing_idx) + offset = 0 + for i in 1:wing_idx-1 + offset += body_aero.wings[i].n_panels + end + return offset .+ (1:body_aero.wings[wing_idx].n_panels) +end + """ unrefined_deform!(body_aero::BodyAerodynamics, theta_angles, delta_angles) @@ -888,8 +912,7 @@ function calculate_results( drag_wing_3D_sum = 0.0 side_wing_3D_sum = 0.0 - # Get wing properties and reference velocity - spanwise_direction = body_aero.wings[1].spanwise_direction + 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 @@ -922,77 +945,62 @@ function calculate_results( @inbounds for k in 1:3 va_ref_unit[k] = va_ref_vec[k] * inv_va_ref end - dir_lift_ref = body_aero.work_vectors[2] - cross3!(dir_lift_ref, va_ref_vec, spanwise_direction) - dir_lift_ref_norm = norm3(dir_lift_ref) - dir_lift_ref_norm > 0.0 || throw(ArgumentError( + norm(cross(va_ref_vec, reference_spanwise)) > 0.0 || throw(ArgumentError( "Reference lift direction is undefined because " * "reference flow is parallel to spanwise direction.")) - @inbounds for k in 1:3 - dir_lift_ref[k] /= dir_lift_ref_norm - end - dir_side_ref = body_aero.work_vectors[3] - cross3!(dir_side_ref, dir_lift_ref, va_ref_unit) + reference_dirs = prescribed_va_directions(SVector{3}(va_ref_vec), reference_spanwise) q_ref = 0.5 * density * va_ref^2 - dir_lift_prescribed_va = body_aero.work_vectors[9] - temp_vec = body_aero.work_vectors[10] - spanwise_unit = SVector{3}(spanwise_direction) - - # Main calculation loop - for (i, panel) in enumerate(panels) - 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 = loads.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 - cross3!(dir_lift_prescribed_va, - panel.va, spanwise_direction) - normalize3!(dir_lift_prescribed_va) - - cross3!(temp_vec, dir_lift_prescribed_va, panel.va) - inv_va_panel = 1.0 / va_panel - @inbounds for k in 1:3 - temp_vec[k] *= inv_va_panel - end - lift_prescribed_va = dot(force, dir_lift_prescribed_va) - drag_prescribed_va = dot(force, panel.va) * inv_va_panel - side_prescribed_va = dot(force, temp_vec) - - lift_wing_3D_sum += lift_prescribed_va * - dot3(dir_lift_prescribed_va, dir_lift_ref) - drag_wing_3D_sum += drag_prescribed_va * - (dot3(panel.va, va_ref_unit) * inv_va_panel) - side_wing_3D_sum += side_prescribed_va * - dot3(temp_vec, dir_side_ref) - - inv_q_area = 1.0 / (q_panel * panel_area) - cl_prescribed_va[i] = lift_prescribed_va * inv_q_area - cd_prescribed_va[i] = drag_prescribed_va * inv_q_area - cs_prescribed_va[i] = side_prescribed_va * inv_q_area - - arm = SVector{3}(panel.aero_center) - SVector{3}(reference_point) - moment = loads.pitching_moment .* axes.y_airf .+ cross(arm, force) - @inbounds for k in 1:3 - f_body_3D[k, i] = force[k] - m_body_3D[k, i] = moment[k] + 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 = loads.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) + 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 * + (dot3(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 + + arm = SVector{3}(panel.aero_center) - SVector{3}(reference_point) + moment = loads.pitching_moment .* axes.y_airf .+ cross(arm, force) + @inbounds for k in 1:3 + f_body_3D[k, i] = force[k] + m_body_3D[k, i] = moment[k] + end end end diff --git a/src/solver.jl b/src/solver.jl index 04bcc7f4..f6e4cd72 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -439,46 +439,44 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; area_all_panels = zero(T) panel_areas = solver.sol.panel_area_dist - # Get wing properties - spanwise_direction = body_aero.wings[1].spanwise_direction - - # Calculate wing geometry properties projected_area = body_aero.projected_area c_ref = body_aero.c_ref - - spanwise_unit = SVector{3, T}(spanwise_direction) - for (i, panel) in enumerate(panels) - panel_area = panel.chord * panel.width - area_all_panels += panel_area - panel_areas[i] = panel_area - - axes = panel_axes(panel) - dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit) - c_span = zero(T) - if solver.is_with_viscous_drag_correction - viscous = spanwise_flow_drag(v_rel_dist[i], solver.lr.v_span_dist[i], - panel.chord, density, solver.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) - lift[i] = loads.lift - drag[i] = loads.drag - panel_moment_dist[i] = loads.moment - - force = loads.force - arm_vec = SVector{3, T}(panel.aero_center) - SVector{3, T}(reference_point) - moment = loads.pitching_moment .* axes.y_airf .+ cross(arm_vec, force) - @inbounds for k in 1:3 - solver.sol.f_body_3D[k, i] = force[k] - solver.sol.m_body_3D[k, i] = moment[k] - end + for (wing_idx, wing) in enumerate(body_aero.wings) + spanwise_unit = SVector{3, T}(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 + panel_areas[i] = panel_area + + axes = panel_axes(panel) + dirs = panel_force_directions(axes, alpha_corrected[i], spanwise_unit) + c_span = zero(T) + if solver.is_with_viscous_drag_correction + viscous = spanwise_flow_drag(v_rel_dist[i], solver.lr.v_span_dist[i], + panel.chord, density, solver.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) + lift[i] = loads.lift + drag[i] = loads.drag + panel_moment_dist[i] = loads.moment + + force = loads.force + arm_vec = SVector{3, T}(panel.aero_center) - SVector{3, T}(reference_point) + moment = loads.pitching_moment .* axes.y_airf .+ cross(arm_vec, force) + @inbounds for k in 1:3 + solver.sol.f_body_3D[k, i] = force[k] + solver.sol.m_body_3D[k, i] = moment[k] + end - arm = (moment_frac - 0.25) * panel.chord - moment_dist[i] = dot(force, axes.z_airf) * arm + loads.pitching_moment + arm = (moment_frac - 0.25) * panel.chord + moment_dist[i] = dot(force, axes.z_airf) * arm + loads.pitching_moment + end end # Python parity: normalize with area-weighted reference velocity for distributed inflow. @@ -525,11 +523,11 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; width_unrefined_dist .= 0.0 fill!(unrefined_count_dist, 0) - panel_idx = 1 for (wing_idx, wing) in enumerate(body_aero.wings) if wing.n_unrefined_sections > 0 section_range = unrefined_section_range(body_aero, wing_idx) - for local_panel_idx in 1:wing.n_panels + wing_panels = panel_range(body_aero, wing_idx) + for (local_panel_idx, panel_idx) in enumerate(wing_panels) panel = body_aero.panels[panel_idx] original_section_idx = wing.refined_panel_mapping[local_panel_idx] target_unrefined_idx = section_range[original_section_idx] @@ -551,7 +549,6 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; width_unrefined_dist[target_unrefined_idx] += panel.width unrefined_count_dist[target_unrefined_idx] += 1 - panel_idx += 1 end # Average coefficients and geometry. width and @@ -573,9 +570,6 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; # sum of panel widths in the unrefined section end end - else - # Skip panels for wings with no unrefined sections - panel_idx += wing.n_panels end end end diff --git a/test/solver/test_wing_directions.jl b/test/solver/test_wing_directions.jl new file mode 100644 index 00000000..937c5cb7 --- /dev/null +++ b/test/solver/test_wing_directions.jl @@ -0,0 +1,59 @@ +using VortexStepMethod +using LinearAlgebra +using Test + +""" + rectangular_wing(rotation, offset) + +A rectangular wing spanning `y`, with every point turned by `rotation` and moved by +`offset`, and its spanwise direction turned with it. +""" +function rectangular_wing(rotation, offset) + chord, span = 1.5, 6.0 + alpha_range = deg2rad.([-10.0, 0.0, 10.0]) + polar = (alpha_range, [-0.6, 0.4, 1.4], fill(0.02, 3), fill(-0.08, 3)) + wing = Wing(10; spanwise_direction=rotation * [0.0, 1.0, 0.0]) + for y in (span / 2, -span / 2) + add_section!(wing, rotation * [0.0, y, 0.0] + offset, + rotation * [chord, y, 0.0] + offset, POLAR_VECTORS, polar) + end + refine!(wing) + return wing +end + +@testset "each wing's panels take their own wing's spanwise direction" begin + # 90° about x leaves the inflow along x unchanged; the second wing is far enough + # away that the two wings do not induce on each other + rotation = [1.0 0.0 0.0; 0.0 0.0 -1.0; 0.0 1.0 0.0] + va = [20.0, 0.0, 0.0] + solo_aero = BodyAerodynamics([rectangular_wing(I(3), zeros(3))]) + pair_aero = BodyAerodynamics([rectangular_wing(I(3), zeros(3)), + rectangular_wing(rotation, [0.0, 1000.0, 0.0])]) + set_va!(solo_aero, va) + set_va!(pair_aero, va) + n_panels = length(solo_aero.panels) + rotated = n_panels .+ (1:n_panels) + + @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 + @test pair.solver_status == FEASIBLE + @test pair.f_body_3D[:, 1:n_panels] ≈ solo.f_body_3D rtol = 1e-5 + @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 + @test pair["F_distribution"][:, rotated] ≈ 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 + end + end +end From c3f4be64c321d38f8911416767e57458ca1a4bff Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 21 Sep 2026 14:10:58 +0200 Subject: [PATCH 2/4] Take the projected span from every wing along the first wing's span direction Co-Authored-By: Claude Opus 5 --- CHANGELOG.md | 5 +++++ docs/src/private_functions.md | 2 ++ src/body_aerodynamics.jl | 24 ++++++++++++------------ src/wing_geometry.jl | 23 +++++++++-------------- test/runtests.jl | 1 + test/solver/test_wing_directions.jl | 17 ++++++++++++++++- 6 files changed, 45 insertions(+), 27 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f65c25f9..ee637142 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -35,6 +35,11 @@ ### Fixed +- On a body whose wings span different directions, such as a wing and a vertical fin, + `solve!`, `solve` and `linearize` take each panel's lift, drag and side directions + from its own wing's `spanwise_direction`, not the first wing's. `solve` computes + `wing_span` and `aspect_ratio_projected` from the extent of all wings along the + first wing's span, through the new `calculate_span(wings, spanwise_direction)`. - `solve!` and `solve` throw a `DimensionMismatch` naming both sizes for a `body_aero` whose panel or unrefined-section count differs from the solver's, where they failed on a broadcast partway through or silently left section results at zero. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 43016830..2a1a9132 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -57,6 +57,7 @@ dynamic_pressure flow_curvature_cm spanwise_flow_drag panel_force_directions +prescribed_va_directions panel_moment panel_couple_force panel_loads @@ -88,6 +89,7 @@ calculate_filaments_for_plotting ```@docs unrefined_deform! unrefined_section_range +panel_range deform! compute_refined_panel_mapping! compute_refined_section_interpolation! diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index fa492024..36326486 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -270,17 +270,6 @@ function unrefined_section_range(body_aero::BodyAerodynamics, wing_idx) return offset .+ (1:body_aero.wings[wing_idx].n_unrefined_sections) end -""" - prescribed_va_directions(va, spanwise) - -Lift and side unit vectors, `(; dir_lift, dir_side)`, of inflow `va` on a wing along -`spanwise`: lift normal to both, side normal to lift and `va`. -""" -@inline function prescribed_va_directions(va, spanwise) - dir_lift = normalize(cross(va, spanwise)) - return (; dir_lift, dir_side=cross(dir_lift, va) / norm(va)) -end - """ panel_range(body_aero::BodyAerodynamics, wing_idx) @@ -804,6 +793,17 @@ function set_pitch_rate_dist!(body_aero::BodyAerodynamics, omega) return nothing end +""" + prescribed_va_directions(va, spanwise) + +Lift and side unit vectors, `(; dir_lift, dir_side)`, of inflow `va` on a wing along +`spanwise`: lift normal to both, side normal to lift and `va`. +""" +@inline function prescribed_va_directions(va, spanwise) + dir_lift = normalize(cross(va, spanwise)) + 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, @@ -1013,7 +1013,7 @@ function calculate_results( # Calculate wing geometry properties projected_area = body_aero.projected_area - wing_span = body_aero.wings[1].span + wing_span = calculate_span(body_aero.wings, reference_spanwise) aspect_ratio_projected = wing_span^2 / projected_area # Calculate Reynolds number diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index 7329dcc2..124409cb 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -1686,24 +1686,19 @@ end """ calculate_span(wing::AbstractWing) + calculate_span(wings, spanwise_direction) -Calculate wing span along spanwise direction. - -Returns: - Float64: Wing span +Extent [m] of the unrefined sections of `wing` along its spanwise direction, or of all +`wings` together along `spanwise_direction`. """ -function calculate_span(wing::AbstractWing) - # Normalize spanwise direction - vector_axis = wing.spanwise_direction ./ norm(wing.spanwise_direction) - - # Get all points - all_points = reduce(vcat, [[section.LE_point, section.TE_point] - for section in wing.unrefined_sections]) - - # Project points and calculate span - projections = [dot(point, vector_axis) for point in all_points] +function calculate_span(wings, spanwise_direction) + axis = spanwise_direction ./ norm(spanwise_direction) + projections = [dot(point, axis) for wing in wings + for section in wing.unrefined_sections + for point in (section.LE_point, section.TE_point)] return maximum(projections) - minimum(projections) end +calculate_span(wing::AbstractWing) = calculate_span((wing,), wing.spanwise_direction) # Project point onto plane @inline function project_onto_plane!(point_proj, point, normal) diff --git a/test/runtests.jl b/test/runtests.jl index b77976d4..9d9122f6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -57,6 +57,7 @@ function include_selected_tests() should_run_test("solver/test_moment_units.jl") && include("solver/test_moment_units.jl") should_run_test("solver/test_forwarddiff.jl") && include("solver/test_forwarddiff.jl") should_run_test("solver/test_backend_comparison.jl") && include("solver/test_backend_comparison.jl") + should_run_test("solver/test_wing_directions.jl") && include("solver/test_wing_directions.jl") should_run_test("solver/test_unrefined_dist.jl") && include("solver/test_unrefined_dist.jl") should_run_test("verification/test_verification.jl") && include("verification/test_verification.jl") should_run_test("VortexStepMethod/test_VortexStepMethod.jl") && include("VortexStepMethod/test_VortexStepMethod.jl") diff --git a/test/solver/test_wing_directions.jl b/test/solver/test_wing_directions.jl index 937c5cb7..9fa13cdc 100644 --- a/test/solver/test_wing_directions.jl +++ b/test/solver/test_wing_directions.jl @@ -51,9 +51,24 @@ end solver = Solver(length(body_aero.panels), 2length(body_aero.wings)) solve(solver, body_aero) end - @test pair["F_distribution"][:, rotated] ≈ rotation * solo["F_distribution"] rtol = 1e-5 + 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 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 end end + +@testset "the projected aspect ratio spans every wing" begin + body_aero = BodyAerodynamics([rectangular_wing(I(3), zeros(3)), + 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), 4) + results = solve(solver, body_aero) + @test results["wing_span"] ≈ 12.0 + @test results["aspect_ratio_projected"] ≈ 12.0^2 / 18.0 +end From 5a7df1d4ca3ebd3dcc9dfe492c76ce14cdcfce46 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 21 Sep 2026 14:16:32 +0200 Subject: [PATCH 3/4] Walk the wings' panel and section offsets with one sum each Co-Authored-By: Claude Opus 5 --- src/body_aerodynamics.jl | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 36326486..7a930108 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -263,11 +263,9 @@ Indices of the unrefined sections of wing `wing_idx` in a distribution that runs unrefined sections of all wings in order, such as `moment_unrefined_dist`. """ function unrefined_section_range(body_aero::BodyAerodynamics, wing_idx) - offset = 0 - for i in 1:wing_idx-1 - offset += body_aero.wings[i].n_unrefined_sections - end - return offset .+ (1:body_aero.wings[wing_idx].n_unrefined_sections) + wings = body_aero.wings + offset = sum(wing.n_unrefined_sections for wing in view(wings, 1:wing_idx-1); init=0) + return offset .+ (1:wings[wing_idx].n_unrefined_sections) end """ @@ -276,11 +274,9 @@ end Indices of the panels of wing `wing_idx` in `body_aero.panels`. """ function panel_range(body_aero::BodyAerodynamics, wing_idx) - offset = 0 - for i in 1:wing_idx-1 - offset += body_aero.wings[i].n_panels - end - return offset .+ (1:body_aero.wings[wing_idx].n_panels) + wings = body_aero.wings + offset = sum(wing.n_panels for wing in view(wings, 1:wing_idx-1); init=0) + return offset .+ (1:wings[wing_idx].n_panels) end """ From f3c2339766503c3c4a90ac38524322abb331b3cb Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Mon, 21 Sep 2026 14:24:50 +0200 Subject: [PATCH 4/4] Check the reference lift direction on the vector it is built from Take the panel's inflow from panel_va throughout the loop, and size the aspect-ratio test's solver from its wing count. Co-Authored-By: Claude Opus 5 --- src/body_aerodynamics.jl | 6 +++--- test/solver/test_wing_directions.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 7a930108..684cde44 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -941,10 +941,10 @@ function calculate_results( @inbounds for k in 1:3 va_ref_unit[k] = va_ref_vec[k] * inv_va_ref end - norm(cross(va_ref_vec, reference_spanwise)) > 0.0 || throw(ArgumentError( + 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.")) - reference_dirs = prescribed_va_directions(SVector{3}(va_ref_vec), reference_spanwise) q_ref = 0.5 * density * va_ref^2 for (wing_idx, wing) in enumerate(body_aero.wings) @@ -982,7 +982,7 @@ function calculate_results( 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 * - (dot3(panel.va, va_ref_unit) * inv_va_panel) + (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) diff --git a/test/solver/test_wing_directions.jl b/test/solver/test_wing_directions.jl index 9fa13cdc..e1ba4e61 100644 --- a/test/solver/test_wing_directions.jl +++ b/test/solver/test_wing_directions.jl @@ -67,7 +67,7 @@ end body_aero = BodyAerodynamics([rectangular_wing(I(3), zeros(3)), 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), 4) + 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