OpenSourceAWE · 1-Bart-1 · Sep 17, 2026 · Sep 16, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,13 @@
 
 ## Unreleased
 
+### Added
+
+- Spanwise-flow viscous drag correction (Gaunaa et al. 2024,
+  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.
+
 ### Changed
 
 - Requires Julia 1.12 or 1.13; 1.10 and 1.11 keep resolving v5.1.1.

diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
@@ -52,6 +52,7 @@ effective_alpha
 panel_inflow
 dynamic_pressure
 flow_curvature_cm
+spanwise_flow_drag
 panel_force_directions
 panel_moment
 panel_couple_force

diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
@@ -748,12 +748,15 @@ end
                       core_radius_fraction, mu, alpha_dist, v_a_dist, chord_dist,
                       x_airf_dist, z_airf_dist, va_dist, va_norm_dist, va_unit_dist,
                       panels::Vector{<:Panel}, is_only_f_and_gamma_output::Bool;
-                      correct_aoa=false, flow_curvature=false)
+                      correct_aoa=false, flow_curvature=false,
+                      is_with_viscous_drag_correction=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`.
+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`.
 
 Returns:
     Dict: Results including forces, coefficients and distributions
@@ -777,6 +780,8 @@ function calculate_results(
     is_only_f_and_gamma_output::Bool;
     correct_aoa::Bool=false,
     flow_curvature::Bool=false,
+    is_with_viscous_drag_correction::Bool=false,
+    v_span_dist=nothing,
 )
 
     n_panels = length(panels)
@@ -893,8 +898,6 @@ function calculate_results(
     cross3!(dir_side_ref, dir_lift_ref, va_ref_unit)
     q_ref = 0.5 * density * va_ref_mag^2
 
-    lift_induced_va = body_aero.work_vectors[7]
-    drag_induced_va = body_aero.work_vectors[8]
     dir_lift_prescribed_va = body_aero.work_vectors[9]
     temp_vec = body_aero.work_vectors[10]
     spanwise_unit = SVector{3}(spanwise_direction)
@@ -906,14 +909,17 @@ function calculate_results(
 
         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_a_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_a_dist[i]),
-            cl_dist[i], cd_dist[i], cm_dist[i])
-        moment_i = loads.moment
-        @inbounds for k in 1:3
-            lift_induced_va[k] = loads.lift * dirs.dir_lift[k]
-            drag_induced_va[k] = loads.drag * dirs.dir_drag[k]
-        end
+            cl_dist[i], cd_dist[i], cm_dist[i]; c_span)
+        force = loads.force
 
         va_panel_mag = va_norm_dist[i]
         va_panel_mag > 0.0 || throw(ArgumentError(
@@ -924,54 +930,32 @@ function calculate_results(
                 panel.va, spanwise_direction)
         normalize3!(dir_lift_prescribed_va)
 
-        lift_prescribed_va =
-            dot3(lift_induced_va, dir_lift_prescribed_va) +
-            dot3(drag_induced_va, dir_lift_prescribed_va)
-        drag_prescribed_va =
-            (dot3(lift_induced_va, panel.va) +
-             dot3(drag_induced_va, panel.va)) / va_panel_mag
         cross3!(temp_vec, dir_lift_prescribed_va, panel.va)
         inv_vpm = 1.0 / va_panel_mag
         @inbounds for k in 1:3
             temp_vec[k] *= inv_vpm
         end
-        side_prescribed_va =
-            dot3(lift_induced_va, temp_vec) +
-            dot3(drag_induced_va, temp_vec)
-
-        width = panel.width
-        @inbounds for k in 1:3
-            f_body_3D[k, i] = (lift_induced_va[k] +
-                               drag_induced_va[k]) * width
-        end
+        lift_prescribed_va = dot(force, dir_lift_prescribed_va)
+        drag_prescribed_va = dot(force, panel.va) * inv_vpm
+        side_prescribed_va = dot(force, temp_vec)
 
-        lift_wing_3D_sum += lift_prescribed_va * width *
+        lift_wing_3D_sum += lift_prescribed_va *
             dot3(dir_lift_prescribed_va, dir_lift_ref)
-        drag_wing_3D_sum += drag_prescribed_va * width *
-            (dot3(panel.va, va_ref_unit) / va_panel_mag)
-        side_wing_3D_sum += side_prescribed_va * width *
+        drag_wing_3D_sum += drag_prescribed_va *
+            (dot3(panel.va, va_ref_unit) * inv_vpm)
+        side_wing_3D_sum += side_prescribed_va *
             dot3(temp_vec, dir_side_ref)
 
-        inv_qc = 1.0 / (q_panel * panel.chord)
-        cl_prescribed_va[i] = lift_prescribed_va * inv_qc
-        cd_prescribed_va[i] = drag_prescribed_va * inv_qc
-        cs_prescribed_va[i] = side_prescribed_va * inv_qc
+        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
 
-        ### Moment ###
-        # r_vector = panel.aero_center - reference_point
-        # M_shift = cross(r_vector, f_body_3D[:,i])
-        # m_body_3D[:,i] = moment_i * panel.y_airf * width + M_shift
-        @inbounds for k in 1:3
-            dir_lift_prescribed_va[k] = panel.aero_center[k] -
-                                        reference_point[k]
-            drag_induced_va[k] = f_body_3D[k, i]
-        end
-        cross3!(temp_vec,
-                dir_lift_prescribed_va, drag_induced_va)
-        local_moment_scale = moment_i * width
+        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
-            m_body_3D[k, i] = local_moment_scale *
-                              panel.y_airf[k] + temp_vec[k]
+            f_body_3D[k, i] = force[k]
+            m_body_3D[k, i] = moment[k]
         end
     end
 

diff --git a/src/panel_aerodynamics.jl b/src/panel_aerodynamics.jl
@@ -167,18 +167,34 @@ the panel normal at its leading and trailing edge.
     scale * width * cm * q_dyn * chord
 
 """
-    panel_loads(axes, dirs, q_dyn, cl, cd, cm, scale=1)
+    spanwise_flow_drag(v_a, v_span, chord, density, mu)
+
+Viscous force increments from spanwise flow (Gaunaa et al. 2024,
+doi:10.1088/1742-6596/2767/2/022068), as `(; delta_cd, c_span)`: an addition to the
+section drag coefficient and a force coefficient along `y_airf`. Both refer to `v_a`,
+the speed normal to the span; `v_span` is the velocity along `y_airf`.
+"""
+@inline function spanwise_flow_drag(v_a, v_span, chord, density, mu)
+    f0 = 0.062 * (density * v_a * chord / mu)^(-1 / 7)
+    skew_factor = (hypot(v_a, v_span) / v_a)^(5 / 7)
+    return (; delta_cd = f0 * (skew_factor - 1), c_span = f0 * v_span / v_a * skew_factor)
+end
+
+"""
+    panel_loads(axes, dirs, q_dyn, cl, cd, cm, scale=1; c_span=0)
 
 Panel load from its polar coefficients, [`panel_axes`](@ref) and
 [`panel_force_directions`](@ref), as `(; lift, drag, moment, force,
 pitching_moment)`. The first three are per unit span; `force` and
-`pitching_moment` are the whole panel's, `scale` included.
+`pitching_moment` are the whole panel's, `scale` included. `c_span` adds a force
+along `y_airf`, as [`spanwise_flow_drag`](@ref) gives it.
 """
-@inline function panel_loads(axes, dirs, q_dyn, cl, cd, cm, scale=1)
-    (; chord, width) = axes
+@inline function panel_loads(axes, dirs, q_dyn, cl, cd, cm, scale=1; c_span=0)
+    (; y_airf, chord, width) = axes
     lift = cl * q_dyn * chord
     drag = cd * q_dyn * chord
     moment = panel_moment(cm, q_dyn, chord)
-    force = (scale * width) .* (lift .* dirs.dir_lift .+ drag .* dirs.dir_drag)
+    force = (scale * width) .* (lift .* dirs.dir_lift .+ drag .* dirs.dir_drag .+
+                                (c_span * q_dyn * chord) .* y_airf)
     return (; lift, drag, moment, force, pitching_moment=scale * width * moment)
 end
diff --git a/src/settings.jl b/src/settings.jl
@@ -180,6 +180,8 @@ Solver configuration, used within [`VSMSettings`](@ref).
     (default `false`)
 - `flow_curvature`: Add the thin-airfoil pitch-rate moment increment to each
     section (default `false`)
+- `is_with_viscous_drag_correction`: Add the spanwise-flow viscous drag and side
+    force to each section (default `false`)
 """
 @with_kw mutable struct SolverSettings
     n_panels::Int64 = 40
@@ -202,6 +204,7 @@ Solver configuration, used within [`VSMSettings`](@ref).
     calc_only_f_and_gamma::Bool=false       # whether to only output f and gamma
     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)
 end
 
 """

diff --git a/src/solver.jl b/src/solver.jl
@@ -85,6 +85,7 @@ end
     gamma_new::MVector{P, T}         = zeros(MVector{P, T})
     alpha_dist::MVector{P, T}        = zeros(MVector{P, T})
     v_a_dist::MVector{P, T}          = zeros(MVector{P, T})
+    v_span_dist::MVector{P, T}       = zeros(MVector{P, T})
 end
 
 @with_kw struct BaseResult{P, T}
@@ -138,6 +139,8 @@ Main solver structure for the Vortex Step Method.See also: [`solve`](@ref)
 - `is_only_f_and_gamma_output`::Bool = false: Whether to only output f and gamma
 - `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
+    and side force to each section, see: [`spanwise_flow_drag`](@ref)
 - `reference_point`::MVec3 = [0.0, 0.0, 0.0]: Moment reference point in body frame
 
 ## Solution
@@ -177,6 +180,7 @@ sol::VSMSolution = VSMSolution(): The result of calling [`solve!`](@ref)
     is_only_f_and_gamma_output::Bool = false
     correct_aoa::Bool = false
     flow_curvature::Bool = false
+    is_with_viscous_drag_correction::Bool = false
     reference_point::MVector{3, T} = zeros(MVector{3, T})
 
     # Intermediate results
@@ -219,6 +223,7 @@ function Solver(body_aero, settings::VSMSettings)
         is_only_f_and_gamma_output=ss.calc_only_f_and_gamma,
         correct_aoa=ss.correct_aoa,
         flow_curvature=ss.flow_curvature,
+        is_with_viscous_drag_correction=ss.is_with_viscous_drag_correction,
         reference_point=reference_point,
     )
 end
@@ -405,9 +410,16 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
 
         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_a_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_a_dist[i]),
-            cl_dist[i], cd_dist[i], cm_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
@@ -593,7 +605,9 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n
         body_aero.panels,
         solver.is_only_f_and_gamma_output;
         correct_aoa=solver.correct_aoa,
-        flow_curvature=solver.flow_curvature
+        flow_curvature=solver.flow_curvature,
+        is_with_viscous_drag_correction=solver.is_with_viscous_drag_correction,
+        v_span_dist=solver.lr.v_span_dist,
     )
     # Attach geometric AoA (already computed in calculate_results) to solver.sol
     if haskey(results, "alpha_geometric")
@@ -755,6 +769,10 @@ end
             x_airf_dist[i,1]*relative_velocity_dist[i,1] +
             x_airf_dist[i,2]*relative_velocity_dist[i,2] +
             x_airf_dist[i,3]*relative_velocity_dist[i,3]
+        solver.lr.v_span_dist[i] =
+            y_airf_dist[i,1]*relative_velocity_dist[i,1] +
+            y_airf_dist[i,2]*relative_velocity_dist[i,2] +
+            y_airf_dist[i,3]*relative_velocity_dist[i,3]
     end
     solver.lr.alpha_dist .= atan.(v_normal_dist, v_tangential_dist)
 
@@ -1240,6 +1258,7 @@ function make_dual_shadow(solver::Solver{P, U, Float64},
         is_only_f_and_gamma_output = 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,
         reference_point = MVector{3, TD}(solver.reference_point),
     )
     return body_aero_d, solver_d

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -53,6 +53,7 @@ function include_selected_tests()
     should_run_test("settings/test_settings.jl") && include("settings/test_settings.jl")
     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_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_viscous_drag_correction.jl b/test/solver/test_viscous_drag_correction.jl
@@ -0,0 +1,100 @@
+using VortexStepMethod
+using VortexStepMethod: spanwise_flow_drag, calc_forces!
+using LinearAlgebra
+using Test
+
+@testset "spanwise_flow_drag is the Gaunaa et al. 2024 model" begin
+    # Gaunaa et al. 2024, J. Phys.: Conf. Ser. 2767 022068, with Re and both
+    # coefficients built on the speed normal to the span:
+    # f0 = 0.062 Re^(-1/7), ΔCd = f0 (cos(β)^(-5/7) - 1), C_par = f0 tan(β) cos(β)^(-5/7)
+    density, mu, chord, speed = 1.225, 1.81e-5, 0.8, 25.0
+    for beta in deg2rad.((-40.0, -5.0, 0.0, 10.0, 35.0))
+        v_normal = speed * cos(beta)
+        f0 = 0.062 * (density * v_normal * chord / mu)^(-1 / 7)
+        correction = spanwise_flow_drag(v_normal, speed * sin(beta), chord, density, mu)
+        @test correction.delta_cd ≈ f0 * (cos(beta)^(-5 / 7) - 1) atol = 1e-14
+        @test correction.c_span ≈ f0 * tan(beta) * cos(beta)^(-5 / 7) atol = 1e-14
+    end
+    @test spanwise_flow_drag(speed, 0.0, chord, density, mu) == (delta_cd=0.0, c_span=0.0)
+end
+
+@testset "Spanwise-flow viscous drag correction in the solver" begin
+    chord, span, V = 1.0, 8.0, 20.0
+    alpha, beta = deg2rad(6.0), deg2rad(12.0)
+
+    wing = Wing(12)
+    for y in range(span / 2, -span / 2, length=5)
+        add_section!(wing, [0.0, y, 0.0], [chord, y, 0.0], INVISCID)
+    end
+    refine!(wing)
+    body_aero = BodyAerodynamics([wing])
+
+    solver_off = Solver(body_aero; use_gamma_prev=false)
+    solver_on = Solver(body_aero; use_gamma_prev=false,
+                       is_with_viscous_drag_correction=true)
+    va_sideslip = V .* [cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)]
+    va_straight = V .* [cos(alpha), 0.0, sin(alpha)]
+
+    function force_dist_at(solver, va)
+        set_va!(body_aero, va)
+        solve!(solver, body_aero)
+        return copy(solver.sol.f_body_3D)
+    end
+
+    @testset "adds the model's drag and spanwise force to each panel" begin
+        density, mu = solver_on.density, solver_on.mu
+        for va in (va_sideslip, va_straight)
+            delta_force = force_dist_at(solver_on, va) .- force_dist_at(solver_off, va)
+            for (i, panel) in enumerate(body_aero.panels)
+                v_normal = solver_on.lr.v_a_dist[i]
+                v_span = solver_on.lr.v_span_dist[i]
+                cos_beta = v_normal / hypot(v_normal, v_span)
+                f0 = 0.062 * (density * v_normal * panel.chord / mu)^(-1 / 7)
+                q_chord_width = 0.5 * density * v_normal^2 * panel.chord * panel.width
+                span_force = dot(delta_force[:, i], panel.y_airf)
+                drag_force = norm(delta_force[:, i] .- span_force .* panel.y_airf)
+                @test span_force ≈ q_chord_width * f0 * (v_span / v_normal) *
+                                   cos_beta^(-5 / 7) atol = 1e-12
+                @test drag_force ≈ q_chord_width * f0 * (cos_beta^(-5 / 7) - 1) atol = 1e-12
+            end
+        end
+    end
+
+    @testset "sideslip drives the spanwise flow and raises the drag" begin
+        delta_force = force_dist_at(solver_on, va_sideslip) .-
+                      force_dist_at(solver_off, va_sideslip)
+        for (i, panel) in enumerate(body_aero.panels)
+            @test solver_on.lr.v_span_dist[i] ≈ dot(va_sideslip, panel.y_airf) rtol = 0.05
+            @test dot(delta_force[:, i], va_sideslip) > 0
+        end
+    end
+
+    @testset "solve reports the corrected forces" begin
+        set_va!(body_aero, va_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_sideslip; zeros(3)]
+        results_for(solver) = VortexStepMethod.linearize(solver, body_aero, y;
+            theta_idxs=nothing, va_idxs=1:3, omega_idxs=4:6)[2]
+        results_on, results_off = results_for(solver_on), results_for(solver_off)
+        @test results_on[1:3] ≈ vec(sum(force_dist_at(solver_on, va_sideslip); dims=2))
+        @test !(results_on[1:3] ≈ results_off[1:3])
+    end
+
+    @testset "calc_forces! stays zero-alloc" begin
+        set_va!(body_aero, va_sideslip)
+        solve!(solver_on, body_aero)
+        calc_forces!(solver_on, body_aero)
+        @test (@allocated calc_forces!(solver_on, body_aero)) == 0
+    end
+
+    @testset "defaults to off" begin
+        @test Solver(body_aero).is_with_viscous_drag_correction == false
+        @test VortexStepMethod.SolverSettings().is_with_viscous_drag_correction == false
+    end
+end