Add the attached-trailed vortex force as an opt-in solver setting

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` 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.

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_at
+attached_trailed_loads
+panel_force_moment
 calculate_stall_angle_list
 wing_span_flip
 spanwise_extent

src/body_aerodynamics.jl

@@ -579,6 +579,81 @@ function update_effective_angle_of_attack!(alpha_corrected,
     nothing
 end
 
+"""
+    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_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]
+    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, 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) + 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)
+    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 +871,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 +904,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 +1003,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 +1029,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 +1056,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]

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
 
 """

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

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")

test/solver/test_attached_trailed_force.jl

@@ -0,0 +1,106 @@
+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 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)
+        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