From 6a18dbfbe2c400254bfe3363306ede0a9408d858 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 23 Sep 2026 11:18:50 +0200 Subject: [PATCH 1/3] Add opt-in force on the attached trailed vortex segments (Gaunaa 2026) Kutta-Joukowski on each panel's two chordwise bound-to-trailing-edge filaments, with the total velocity at their three-quarter-chord point, behind is_with_attached_trailed_force (default false). solve! and calculate_results now assemble each panel's force and moment through one panel_force_moment. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 5 + docs/src/private_functions.md | 3 + src/body_aerodynamics.jl | 89 +++++++++++++++-- src/settings.jl | 3 + src/solver.jl | 14 ++- test/runtests.jl | 1 + test/solver/test_attached_trailed_force.jl | 107 +++++++++++++++++++++ 7 files changed, 208 insertions(+), 14 deletions(-) create mode 100644 test/solver/test_attached_trailed_force.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 4578e6ad..d7f1afd0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,11 @@ doi:10.1088/1742-6596/2767/2/022068): each section gets a drag increment and a force along its span from the flow across it, in `solve!`, `solve` and `linearize`. Opt-in via `is_with_viscous_drag_correction` (default `false`) on the solver settings. +- Force on the chordwise trailed vortex segments bound to each panel (Gaunaa, Li & + Pirrung 2026, doi:10.1088/1742-6596/3224/4/042030): Kutta–Joukowski with the total + velocity at their three-quarter-chord point, in `solve!`, `solve` and `linearize`. It + lowers `CL` by about 2 % on a swept, anhedral wing whose sections are perpendicular to + the quarter-chord line. Opt-in via `is_with_attached_trailed_force` (default `false`). - `plot_section_polars(body_aero; panels, alphas, delta)` draws cl, cd and cm against α per panel through `calculate_cl`/`calculate_cd`/`calculate_cm`, for every aero model and at flap deflection `delta`, in one figure instead of one coefficient per call. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index f7482533..e5784e62 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -34,6 +34,9 @@ bisect_sign_change calculate_relative_alpha_and_velocity calculate_relative_alpha_and_relative_velocity update_effective_angle_of_attack! +induced_velocity +attached_trailed_loads +panel_force_moment calculate_stall_angle_list wing_span_flip spanwise_extent diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 2ce94f20..7558fbab 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -579,6 +579,77 @@ function update_effective_angle_of_attack!(alpha_corrected, nothing end +""" + induced_velocity(body_aero::BodyAerodynamics, point, gamma, core_radius_fraction) + +Velocity [m/s] induced at `point` by the horseshoe vortices of all panels, panel `j` +carrying circulation `gamma[j]` [m²/s] and its frozen wake. +""" +function induced_velocity(body_aero::BodyAerodynamics{P, W, T}, point, gamma, + core_radius_fraction) where {P, W, T} + velocity_ring = body_aero.work_vectors[8] + velocity_filament = body_aero.work_vectors[9] + velocity = zero(SVector{3, T}) + for (j, panel) in enumerate(body_aero.panels) + wake = panel.filaments[4] + calculate_velocity_induced_single_ring_semiinfinite!(velocity_ring, + velocity_filament, panel.filaments, point, false, wake.va, wake.direction, + gamma[j], core_radius_fraction, body_aero.work_vectors) + velocity += SVector{3, T}(velocity_ring) + end + return velocity +end + +""" + attached_trailed_loads(body_aero::BodyAerodynamics, i, gamma, density, + core_radius_fraction, reference_point) + +Kutta–Joukowski force [N] on the two chordwise trailed vortex segments of panel `i`, +from its quarter-chord bound points to its trailing edge, and its moment [N·m] about +`reference_point`, as `(; force, moment)`. Each segment carries `gamma[i]` and sees the +panel inflow plus the [`induced_velocity`](@ref) at its three-quarter-chord point. +""" +function attached_trailed_loads(body_aero::BodyAerodynamics{P, W, T}, i, gamma, density, + core_radius_fraction, reference_point) where {P, W, T} + panel = body_aero.panels[i] + force = zero(SVector{3, T}) + moment = zero(SVector{3, T}) + segments = ((panel.bound_point_1, panel.TE_point_1, 1), + (panel.bound_point_2, panel.TE_point_2, -1)) + for (bound_point, te_point, orientation) in segments + chordwise = SVector{3, T}(te_point) - SVector{3, T}(bound_point) + point = SVector{3, T}(bound_point) + (2 / 3) * chordwise + velocity = SVector{3, T}(panel.va_vec) + + induced_velocity(body_aero, point, gamma, core_radius_fraction) + segment_force = (density * orientation * gamma[i]) * cross(velocity, chordwise) + force += segment_force + moment += cross(point - SVector{3, T}(reference_point), segment_force) + end + return (; force, moment) +end + +""" + panel_force_moment(body_aero::BodyAerodynamics, i, loads, y_airf, gamma, density, + core_radius_fraction, reference_point, + is_with_attached_trailed_force) + +Force [N] and moment [N·m] about `reference_point` of panel `i`: its section `loads` +from [`panel_loads`](@ref) acting at the aerodynamic centre, plus +[`attached_trailed_loads`](@ref) when `is_with_attached_trailed_force`. +Returns `(; force, moment)`. +""" +function panel_force_moment(body_aero::BodyAerodynamics{P, W, T}, i, loads, y_airf, gamma, + density, core_radius_fraction, reference_point, + is_with_attached_trailed_force) where {P, W, T} + arm = SVector{3, T}(body_aero.panels[i].aero_center) - SVector{3, T}(reference_point) + force = loads.force + moment = loads.pitching_moment .* y_airf .+ cross(arm, force) + is_with_attached_trailed_force || return (; force, moment) + attached = attached_trailed_loads(body_aero, i, gamma, density, core_radius_fraction, + reference_point) + return (; force=force + attached.force, moment=moment + attached.moment) +end + @inline function intersect_line_with_plane( x_cp, f_unit, plane_point, plane_normal; tol=1e-6 ) @@ -796,14 +867,15 @@ end 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, v_span_dist=nothing) + 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`. +`v_span_dist`. `is_with_attached_trailed_force` adds [`attached_trailed_loads`](@ref). Returns: Dict: Results including forces, coefficients and distributions @@ -828,6 +900,7 @@ function calculate_results( 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, ) @@ -926,11 +999,7 @@ function calculate_results( va_ref = norm3(va_ref_vec) va_ref > 0.0 || throw(ArgumentError( "Reference freestream magnitude must be positive.")) - va_ref_unit = body_aero.work_vectors[1] - inv_va_ref = 1.0 / va_ref - @inbounds for k in 1:3 - va_ref_unit[k] = va_ref_vec[k] * inv_va_ref - end + 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 " * @@ -956,7 +1025,9 @@ function calculate_results( 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 + (; 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( @@ -981,8 +1052,6 @@ function calculate_results( 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] diff --git a/src/settings.jl b/src/settings.jl index 2e73deb7..354d76c9 100644 --- a/src/settings.jl +++ b/src/settings.jl @@ -179,6 +179,8 @@ Solver configuration, used within [`VSMSettings`](@ref). section (default `false`) - `is_with_viscous_drag_correction`: Add the spanwise-flow viscous drag and side force to each section (default `false`) +- `is_with_attached_trailed_force`: Add the force on the chordwise trailed vortex + segments bound to each panel (default `false`) """ @with_kw mutable struct SolverSettings n_panels::Int64 = 40 @@ -199,6 +201,7 @@ Solver configuration, used within [`VSMSettings`](@ref). 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) + is_with_attached_trailed_force::Bool=false # chordwise bound-vortex force (Gaunaa 2026) end """ diff --git a/src/solver.jl b/src/solver.jl index 7cc8810c..9357dfde 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -137,6 +137,8 @@ Main solver structure for the Vortex Step Method.See also: [`solve`](@ref) increment `-(π/4) q̂` to each section, see: [`flow_curvature_cm`](@ref) - `is_with_viscous_drag_correction`::Bool = false: Add the spanwise-flow viscous drag and side force to each section, see: [`spanwise_flow_drag`](@ref) +- `is_with_attached_trailed_force`::Bool = false: Add the force on the chordwise trailed + vortex segments bound to each panel, see: [`attached_trailed_loads`](@ref) - `reference_point`::MVec3 = [0.0, 0.0, 0.0]: Moment reference point in body frame ## Solution @@ -173,6 +175,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref) correct_aoa::Bool = false flow_curvature::Bool = false is_with_viscous_drag_correction::Bool = false + is_with_attached_trailed_force::Bool = false reference_point::MVector{3, T} = zeros(MVector{3, T}) # Intermediate results @@ -260,6 +263,7 @@ function solver_kwargs(solver_settings::SolverSettings) correct_aoa=solver_settings.correct_aoa, flow_curvature=solver_settings.flow_curvature, is_with_viscous_drag_correction=solver_settings.is_with_viscous_drag_correction, + is_with_attached_trailed_force=solver_settings.is_with_attached_trailed_force, ) end @@ -456,16 +460,16 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics; 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) + (; force, moment) = panel_force_moment(body_aero, i, loads, axes.y_airf, + gamma_new, density, solver.core_radius_fraction, reference_point, + solver.is_with_attached_trailed_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 + moment_dist[i] = dot(loads.force, axes.z_airf) * arm + loads.pitching_moment end end @@ -634,6 +638,7 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n 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 @@ -1237,6 +1242,7 @@ function make_dual_shadow(solver::Solver{P, U, Float64}, 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, reference_point = MVector{3, TD}(solver.reference_point), ) return body_aero_d, solver_d diff --git a/test/runtests.jl b/test/runtests.jl index 76981ba9..d119cdc0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -54,6 +54,7 @@ function include_selected_tests() should_run_test("solver/test_solver.jl") && include("solver/test_solver.jl") should_run_test("solver/test_flow_curvature.jl") && include("solver/test_flow_curvature.jl") should_run_test("solver/test_viscous_drag_correction.jl") && include("solver/test_viscous_drag_correction.jl") + should_run_test("solver/test_attached_trailed_force.jl") && include("solver/test_attached_trailed_force.jl") 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") diff --git a/test/solver/test_attached_trailed_force.jl b/test/solver/test_attached_trailed_force.jl new file mode 100644 index 00000000..318e1374 --- /dev/null +++ b/test/solver/test_attached_trailed_force.jl @@ -0,0 +1,107 @@ +using VortexStepMethod +using VortexStepMethod: attached_trailed_loads, calc_forces!, solver_kwargs +using DifferentiationInterface +using LinearAlgebra +using Test + +""" + curved_elliptic_wing(; sweep, anhedral, aspect_ratio=4.0, n_sections=61) + +Elliptic wing of unit semi-span whose quarter-chord line runs back by `sweep` and down +by `anhedral` at the tips, parabolic in span, with each section perpendicular to that +line. Returns the wing and its aspect ratio. +""" +function curved_elliptic_wing(; sweep, anhedral, aspect_ratio=4.0, n_sections=61) + max_chord = 8 / (π * aspect_ratio) + wing = Wing(n_sections - 1; spanwise_distribution=UNCHANGED) + for theta in range(0, π, length=n_sections) + eta = cos(theta) * (1 - 1e-5) + quarter_chord = [sweep * eta^2, eta, -anhedral * eta^2] + tangent = normalize([2sweep * eta, 1.0, -2anhedral * eta]) + chordwise = normalize([1.0, 0.0, 0.0] .- tangent[1] .* tangent) + chord = max_chord * sqrt(1 - eta^2) + add_section!(wing, quarter_chord .- 0.25chord .* chordwise, + quarter_chord .+ 0.75chord .* chordwise, INVISCID) + end + refine!(wing) + return wing, aspect_ratio +end + +""" + lift_and_oswald(wing, aspect_ratio, model, is_with_attached_trailed_force) + +Wing `CL` and blade-force Oswald factor `CL²/(π AR CD)` at α = 10°. +""" +function lift_and_oswald(wing, aspect_ratio, model, is_with_attached_trailed_force) + body_aero = BodyAerodynamics([wing]) + 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"]) +end + +@testset "Attached trailed vortex force" begin + @testset "defaults to off and follows the solver settings" begin + @test VortexStepMethod.SolverSettings().is_with_attached_trailed_force == false + @test Solver(4, 2).is_with_attached_trailed_force == false + settings = VortexStepMethod.SolverSettings(is_with_attached_trailed_force=true) + @test solver_kwargs(settings).is_with_attached_trailed_force + end + + # Gaunaa, Li & Pirrung, TORQUE 2026, §5.1-5.2: omitting the force overpredicts CL + # by about 2 % on a swept, anhedral wing and barely changes a straight one. + @testset "lowers lift by about 2 % on a swept, anhedral wing ($model)" for + model in (VSM, LLT) + straight, aspect_ratio = curved_elliptic_wing(; sweep=0.0, anhedral=0.0) + CL_off, _ = lift_and_oswald(straight, aspect_ratio, model, false) + CL_on, _ = lift_and_oswald(straight, aspect_ratio, model, true) + @test abs(CL_on / CL_off - 1) < 0.005 + + curved, aspect_ratio = curved_elliptic_wing(; sweep=0.1, anhedral=0.25) + CL_off, oswald_off = lift_and_oswald(curved, aspect_ratio, model, false) + CL_on, oswald_on = lift_and_oswald(curved, aspect_ratio, model, true) + @test -0.03 < CL_on / CL_off - 1 < -0.015 + @test oswald_on > oswald_off + end + + wing, _ = curved_elliptic_wing(; sweep=0.1, anhedral=0.25, n_sections=21) + va_vec = 20.0 .* [cosd(10), 0.0, sind(10)] + + @testset "adds attached_trailed_loads to every panel in solve! and solve" begin + body_aero = BodyAerodynamics([wing]) + set_va!(body_aero, va_vec) + reference_point = [0.1, 0.0, -0.2] + solver_off = Solver(wing.n_panels, wing.n_unrefined_sections; use_gamma_prev=false, + reference_point) + solver_on = Solver(wing.n_panels, wing.n_unrefined_sections; use_gamma_prev=false, + reference_point, is_with_attached_trailed_force=true) + force_off = copy(solve!(solver_off, body_aero).f_body_3D) + moment_off = copy(solver_off.sol.m_body_3D) + solve!(solver_on, body_aero) + for i in eachindex(body_aero.panels) + attached = attached_trailed_loads(body_aero, i, solver_on.lr.gamma_new, + solver_on.density, solver_on.core_radius_fraction, reference_point) + @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 + + @testset "ForwardDiff matches FiniteDiff with the force on" begin + body_aero = BodyAerodynamics([wing]) + set_va!(body_aero, va_vec) + solver = Solver(wing.n_panels, wing.n_unrefined_sections; use_gamma_prev=false, + type_initial_gamma_distribution=ELLIPTIC, is_with_attached_trailed_force=true) + jac_fwd, _ = VortexStepMethod.linearize(solver, body_aero, va_vec; + theta_idxs=nothing, va_vec_idxs=1:3, aero_coeffs=true, + backend=AutoForwardDiff()) + jac_fd, _ = VortexStepMethod.linearize(solver, body_aero, va_vec; + theta_idxs=nothing, va_vec_idxs=1:3, aero_coeffs=true, + backend=AutoFiniteDiff(absstep=1e-5, relstep=1e-5)) + @test maximum(abs.(jac_fwd .- jac_fd)) / maximum(abs, jac_fwd) < 1e-3 + end +end From d589a7deed21e0cf20acb19bea47504a257b9f48 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 23 Sep 2026 11:22:28 +0200 Subject: [PATCH 2/3] Shorten the paper reference in the attached trailed force test Co-Authored-By: Claude Opus 5.5 --- test/solver/test_attached_trailed_force.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/solver/test_attached_trailed_force.jl b/test/solver/test_attached_trailed_force.jl index 318e1374..610262b7 100644 --- a/test/solver/test_attached_trailed_force.jl +++ b/test/solver/test_attached_trailed_force.jl @@ -49,8 +49,7 @@ end @test solver_kwargs(settings).is_with_attached_trailed_force end - # Gaunaa, Li & Pirrung, TORQUE 2026, §5.1-5.2: omitting the force overpredicts CL - # by about 2 % on a swept, anhedral wing and barely changes a straight one. + # Gaunaa, Li & Pirrung 2026 §5: without the force CL is ~2 % high on the curved wing. @testset "lowers lift by about 2 % on a swept, anhedral wing ($model)" for model in (VSM, LLT) straight, aspect_ratio = curved_elliptic_wing(; sweep=0.0, anhedral=0.0) From 675f0e7b928973da7ef7657503297b50a64fd3f7 Mon Sep 17 00:00:00 2001 From: 1-Bort-1 <323661610+1-Bort-1@users.noreply.github.com> Date: Wed, 23 Sep 2026 11:34:00 +0200 Subject: [PATCH 3/3] Turn the attached trailed inflow by omega at its sample point Name the three-quarter-chord fraction, rename induced_velocity to induced_velocity_at, and drop the benchmark number from the changelog. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 4 ++-- docs/src/private_functions.md | 2 +- src/body_aerodynamics.jl | 18 +++++++++++------- 3 files changed, 14 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d7f1afd0..50b26ddd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,8 +24,8 @@ - Force on the chordwise trailed vortex segments bound to each panel (Gaunaa, Li & Pirrung 2026, doi:10.1088/1742-6596/3224/4/042030): Kutta–Joukowski with the total velocity at their three-quarter-chord point, in `solve!`, `solve` and `linearize`. It - lowers `CL` by about 2 % on a swept, anhedral wing whose sections are perpendicular to - the quarter-chord line. Opt-in via `is_with_attached_trailed_force` (default `false`). + lowers `CL` on a swept, anhedral wing whose sections are perpendicular to the + quarter-chord line. Opt-in via `is_with_attached_trailed_force` (default `false`). - `plot_section_polars(body_aero; panels, alphas, delta)` draws cl, cd and cm against α per panel through `calculate_cl`/`calculate_cd`/`calculate_cm`, for every aero model and at flap deflection `delta`, in one figure instead of one coefficient per call. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index e5784e62..ea497f91 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -34,7 +34,7 @@ bisect_sign_change calculate_relative_alpha_and_velocity calculate_relative_alpha_and_relative_velocity update_effective_angle_of_attack! -induced_velocity +induced_velocity_at attached_trailed_loads panel_force_moment calculate_stall_angle_list diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index 7558fbab..6cfdceec 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -580,12 +580,12 @@ function update_effective_angle_of_attack!(alpha_corrected, end """ - induced_velocity(body_aero::BodyAerodynamics, point, gamma, core_radius_fraction) + induced_velocity_at(body_aero::BodyAerodynamics, point, gamma, core_radius_fraction) Velocity [m/s] induced at `point` by the horseshoe vortices of all panels, panel `j` carrying circulation `gamma[j]` [m²/s] and its frozen wake. """ -function induced_velocity(body_aero::BodyAerodynamics{P, W, T}, point, gamma, +function induced_velocity_at(body_aero::BodyAerodynamics{P, W, T}, point, gamma, core_radius_fraction) where {P, W, T} velocity_ring = body_aero.work_vectors[8] velocity_filament = body_aero.work_vectors[9] @@ -606,21 +606,25 @@ end Kutta–Joukowski force [N] on the two chordwise trailed vortex segments of panel `i`, from its quarter-chord bound points to its trailing edge, and its moment [N·m] about -`reference_point`, as `(; force, moment)`. Each segment carries `gamma[i]` and sees the -panel inflow plus the [`induced_velocity`](@ref) at its three-quarter-chord point. +`reference_point`, as `(; force, moment)`. Each segment carries `gamma[i]` and sees, at +its three-quarter-chord point, the inflow turned by `body_aero.omega` plus the +[`induced_velocity_at`](@ref) that point. """ function attached_trailed_loads(body_aero::BodyAerodynamics{P, W, T}, i, gamma, density, core_radius_fraction, reference_point) where {P, W, T} panel = body_aero.panels[i] + three_quarter_chord = (0.75 - 0.25) / (1 - 0.25) force = zero(SVector{3, T}) moment = zero(SVector{3, T}) segments = ((panel.bound_point_1, panel.TE_point_1, 1), (panel.bound_point_2, panel.TE_point_2, -1)) for (bound_point, te_point, orientation) in segments chordwise = SVector{3, T}(te_point) - SVector{3, T}(bound_point) - point = SVector{3, T}(bound_point) + (2 / 3) * chordwise - velocity = SVector{3, T}(panel.va_vec) + - induced_velocity(body_aero, point, gamma, core_radius_fraction) + point = SVector{3, T}(bound_point) + three_quarter_chord * chordwise + arm = point - SVector{3, T}(panel.control_point) + inflow = SVector{3, T}(panel.va_vec) - cross(SVector{3, T}(body_aero.omega), arm) + velocity = inflow + + induced_velocity_at(body_aero, point, gamma, core_radius_fraction) segment_force = (density * orientation * gamma[i]) * cross(velocity, chordwise) force += segment_force moment += cross(point - SVector{3, T}(reference_point), segment_force)