Per-panel pitch rates for flow curvature

CHANGELOG.md

@@ -1,5 +1,28 @@
 # Changelog
 
+## Unreleased
+
+### Added
+- `flow_curvature` solver setting (default `false`). When enabled, each section
+  gets the thin-airfoil pitch-rate moment increment `Δcm = -(π/4) q̂` with
+  `q̂ = q c / (2 v_rel)` and `q = ω ⋅ y_airf`. A section rotating about its own
+  spanwise axis sees an incidence that varies linearly along the chord, which is
+  equivalent to parabolic camber and produces a quarter-chord moment that a
+  single control point cannot represent. The lift response to `q` was already
+  exact because the inflow is sampled at the three-quarter-chord point, so only
+  the moment was missing.
+- `pitch_rate_dist` field on `BodyAerodynamics`: each panel's rotation rate about
+  its own spanwise axis, which is what `flow_curvature` reads. `set_va!(body_aero,
+  va, omega)` fills it by projecting `omega` onto every panel's `y_airf`, so
+  panels at different dihedral see different rates from one body rate. The
+  distributed `set_va!(body_aero, va_distribution; pitch_rate_dist)` takes it
+  directly, so twist and flapping rates of a deforming wing — which no single
+  body rate can express — reach the moment. Omitting the keyword zeroes it rather
+  than reusing a stale `omega`.
+- `section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord)` builds
+  one entry of that distribution from a section's edge velocities, and reduces to
+  `ω ⋅ y_airf` for rigid motion.
+
 ## VortexStepMethod v4.0.0 2026-08-03
 
 ### Added

docs/src/functions.md

@@ -82,6 +82,7 @@ CurrentModule = VortexStepMethod
 ## Setting the inflow conditions and solving
 ```@docs
 set_va!
+section_pitch_rate
 solve
 solve!
 reinit!(body_aero::BodyAerodynamics{P, W, T}) where {P, W, T}

docs/src/private_functions.md

@@ -17,6 +17,8 @@ calculate_cl
 calculate_cd
 calculate_cm
 calculate_cd_cm
+flow_curvature_cm
+set_pitch_rate_dist!
 calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!

src/VortexStepMethod.jl

@@ -30,7 +30,7 @@ export ObjWing, Section, Wing, refine!, reinit!
 export BodyAerodynamics
 export Solver, VSMSolution, linearize, solve, solve!, solve_base!, calc_forces!
 export calculate_results
-export add_section!, set_va!
+export add_section!, set_va!, section_pitch_rate
 export calculate_projected_area, calculate_span
 export MVec3
 

src/body_aerodynamics.jl

@@ -15,6 +15,9 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru
 - `stall_angle_list`=zeros(Float64, P):  stall angle per panel
 - `alpha_dist::MVector{P, Float64}` = zeros(Float64, P)
 - `v_a_dist::MVector{P, Float64}` = zeros(Float64, P)
+- `pitch_rate_dist::MVector{P, Float64}` = zeros(Float64, P): rotation rate of each
+    panel about its own spanwise axis, positive nose-up [rad/s]; set by
+    [set_va!](@ref) and read when the solver has `flow_curvature` enabled
 - `work_vectors`::NTuple{10, MVec3} = ntuple(_ -> zeros(MVec3), 10)
 - `AIC::Array{Float64, 3}` = zeros(3, P, P)
 - `projected_area::Float64` = 1.0: The area projected onto the xy-plane of the kite body reference frame [m²]
@@ -34,6 +37,7 @@ Main structure for calculating aerodynamic properties of bodies. Use the constru
     stall_angle_list::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})
+    pitch_rate_dist::MVector{P, T} = zeros(MVector{P, T})
     work_vectors::NTuple{10, MVector{3, T}} = ntuple(_ -> zeros(MVector{3, T}), 10)
     AIC::Array{T, 3} = zeros(T, 3, P, P)
     projected_area::T = one(T)
@@ -714,7 +718,53 @@ function compute_panel_center_of_pressures(
 end
 
 """
-    calculate_results(body_aero::BodyAerodynamics, gamma_new, 
+    flow_curvature_cm(pitch_rate, chord, v_rel)
+
+Quarter-chord moment increment of a section rotating about its own spanwise axis,
+from thin airfoil theory. The rotation makes the local incidence vary linearly
+along the chord, which is equivalent to parabolic camber and yields
+`Δcm = -(π/4) q̂` with `q̂ = q c / (2 v_rel)` and `q` positive nose-up.
+Independent of the pivot location; the lift response to `q` needs no correction
+because the inflow is already sampled at the three-quarter-chord control point.
+"""
+@inline function flow_curvature_cm(pitch_rate, chord, v_rel)
+    v_rel > 0 || return zero(chord)
+    return -0.25π * pitch_rate * chord / (2v_rel)
+end
+
+"""
+    section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord)
+
+Rate at which a section rotates about its own spanwise axis, from the velocities
+of its leading and trailing edge. Positive nose-up, matching
+[`flow_curvature_cm`](@ref). Use this to build a `pitch_rate_dist` for
+[`set_va!`](@ref) from a deforming structure, where twist and flapping rates
+differ per section and no single body rate describes them.
+"""
+@inline function section_pitch_rate(velocity_leading, velocity_trailing,
+                                    z_airf, chord)
+    chord > 0 || return zero(chord)
+    normal_rate = dot3(velocity_trailing, z_airf) -
+                  dot3(velocity_leading, z_airf)
+    return -normal_rate / chord
+end
+
+"""
+    set_pitch_rate_dist!(body_aero, omega)
+
+Fill `body_aero.pitch_rate_dist` from a rigid-body turn rate by projecting it
+onto each panel's own spanwise axis. Panels with different dihedral see
+different rates from the same `omega`.
+"""
+function set_pitch_rate_dist!(body_aero::BodyAerodynamics, omega)
+    for (i, panel) in enumerate(body_aero.panels)
+        body_aero.pitch_rate_dist[i] = dot3(omega, panel.y_airf)
+    end
+    return nothing
+end
+
+"""
+    calculate_results(body_aero::BodyAerodynamics, gamma_new,
                      density,
                      core_radius_fraction, mu,
                      alpha_dist, v_a_dist,
@@ -726,6 +776,9 @@ end
 
 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`.
+
 Returns:
     Dict: Results including forces, coefficients and distributions
 """
@@ -747,6 +800,7 @@ function calculate_results(
     panels::Vector{<:Panel},
     is_only_f_and_gamma_output::Bool;
     correct_aoa::Bool=false,
+    flow_curvature::Bool=false,
 )
 
     n_panels = length(panels)
@@ -775,6 +829,10 @@ function calculate_results(
         cl_array[i] = calculate_cl(panel, alpha_dist[i])
         cd_array[i], cm_array[i] = calculate_cd_cm(
             panel, alpha_dist[i])
+        if flow_curvature
+            cm_array[i] += flow_curvature_cm(
+                body_aero.pitch_rate_dist[i], chord_array[i], v_a_dist[i])
+        end
         panel_width_array[i] = panel.width
         va_norm = va_norm_array[i]
         x_norm = norm3(panel.x_airf)
@@ -1058,11 +1116,15 @@ Set velocity array and update wake filaments.
 - body_aero::BodyAerodynamics: The [BodyAerodynamics](@ref) struct to modify
 - `va::VelVector`: Velocity vector of the apparent wind speed           [m/s]
 - `omega::VelVector`: Turn rate vector around x y and z axis            [rad/s]
+
+`omega` is also projected onto each panel's spanwise axis into
+`pitch_rate_dist`, which the solver reads when `flow_curvature` is enabled.
 """
 function set_va!(body_aero::BodyAerodynamics{P, W, T}, va::AbstractVector, omega=zeros(MVector{3, T})) where {P, W, T}
     n_panels = length(body_aero.panels)
     va_distribution = zeros(T, n_panels, 3)
     body_aero.omega .= omega
+    set_pitch_rate_dist!(body_aero, omega)
 
     if all(iszero, omega)
         va_distribution .= reshape(va, 1, 3)
@@ -1092,9 +1154,28 @@ function set_va!(body_aero::BodyAerodynamics{P, W, T}, va::AbstractVector, omega
     return nothing
 end
 
-function set_va!(body_aero::BodyAerodynamics, va_distribution::AbstractMatrix)
+"""
+    set_va!(body_aero::BodyAerodynamics, va_distribution::AbstractMatrix;
+            pitch_rate_dist=nothing)
+
+Set a per-panel inflow distribution. `pitch_rate_dist` gives each panel's rotation
+rate about its own spanwise axis [rad/s], positive nose-up; build it with
+[`section_pitch_rate`](@ref) when the structure deforms, since twist and flapping
+rates differ per section and no single body rate describes them. It is reset to
+zero when omitted, because this method takes no `omega` and a stale one would
+silently feed the `flow_curvature` moment.
+"""
+function set_va!(body_aero::BodyAerodynamics, va_distribution::AbstractMatrix;
+                 pitch_rate_dist=nothing)
     size(va_distribution, 1) != length(body_aero.panels) &&
         throw(ArgumentError("Number of rows in va distribution should be equal to number of panels."))
+    if isnothing(pitch_rate_dist)
+        body_aero.pitch_rate_dist .= 0
+    else
+        length(pitch_rate_dist) != length(body_aero.panels) &&
+            throw(ArgumentError("Length of pitch rate distribution should be equal to number of panels."))
+        body_aero.pitch_rate_dist .= pitch_rate_dist
+    end
 
     for (i, panel) in enumerate(body_aero.panels)
         panel.va .= va_distribution[i, :]

src/settings.jl

@@ -89,6 +89,8 @@ Solver configuration, used within [`VSMSettings`](@ref).
     (default `false`)
 - `correct_aoa`: Perform angle of attack correction
     (default `false`)
+- `flow_curvature`: Add the thin-airfoil pitch-rate moment increment to each
+    section (default `false`)
 """
 @with_kw mutable struct SolverSettings
     n_panels::Int64 = 40
@@ -110,6 +112,7 @@ Solver configuration, used within [`VSMSettings`](@ref).
     mu::Float64 = 1.81e-5                   # dynamic viscosity [N·s/m²]
     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
 end
 
 """

src/solver.jl

@@ -135,6 +135,8 @@ Main solver structure for the Vortex Step Method.See also: [solve](@ref)
 - `core_radius_fraction`::Float64 = 1e-20: 
 - 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
+- `flow_curvature`::Bool = false: Add the thin-airfoil pitch-rate moment
+    increment `-(π/4) q̂` to each section, see: [flow_curvature_cm](@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)
     mu::T = T(1.81e-5)
     is_only_f_and_gamma_output::Bool = false
     correct_aoa::Bool = false
+    flow_curvature::Bool = false
     reference_point::MVector{3, T} = zeros(MVector{3, T})
 
     # Intermediate results
@@ -214,6 +217,7 @@ function Solver(body_aero, settings::VSMSettings)
         mu=ss.mu,
         is_only_f_and_gamma_output=ss.calc_only_f_and_gamma,
         correct_aoa=ss.correct_aoa,
+        flow_curvature=ss.flow_curvature,
         reference_point=reference_point,
     )
 end
@@ -288,6 +292,11 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
     for (i, panel) in enumerate(panels)                                               # zero bytes
         cl_dist[i] = calculate_cl(panel, alpha_dist[i])
         cd_dist[i], cm_dist[i] = calculate_cd_cm(panel, alpha_dist[i])
+        if solver.flow_curvature
+            cm_dist[i] += flow_curvature_cm(
+                body_aero.pitch_rate_dist[i], solver.sol._chord_dist[i],
+                v_a_dist[i])
+        end
         width_dist[i] = panel.width
 
         # Geometric AoA using panel-local axes and prescribed
@@ -599,7 +608,8 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n
         solver.br.va_unit_dist,
         body_aero.panels,
         solver.is_only_f_and_gamma_output;
-        correct_aoa=solver.correct_aoa
+        correct_aoa=solver.correct_aoa,
+        flow_curvature=solver.flow_curvature
     )
     # Attach geometric AoA (already computed in calculate_results) to solver.sol
     if haskey(results, "alpha_geometric")
@@ -1218,6 +1228,7 @@ function make_dual_shadow(solver::Solver{P, U, Float64},
         mu = TD(solver.mu),
         is_only_f_and_gamma_output = solver.is_only_f_and_gamma_output,
         correct_aoa = solver.correct_aoa,
+        flow_curvature = solver.flow_curvature,
         reference_point = MVector{3, TD}(solver.reference_point),
     )
     return body_aero_d, solver_d

test/runtests.jl

@@ -52,6 +52,7 @@ function include_selected_tests()
     should_run_test("ram_geometry/test_kite_geometry.jl") && include("ram_geometry/test_kite_geometry.jl")
     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_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_unrefined_dist.jl") && include("solver/test_unrefined_dist.jl")

test/solver/test_flow_curvature.jl

@@ -0,0 +1,108 @@
+using VortexStepMethod
+using VortexStepMethod: flow_curvature_cm
+using LinearAlgebra
+using Test
+
+@testset "section_pitch_rate reduces to the rigid-body rate" begin
+    # a rigid section rotating at q about y_airf must give back exactly q, which
+    # is what lets one expression serve both rigid and deforming wings
+    x_airf, y_airf, z_airf = [1.0, 0, 0], [0, 1.0, 0], [0, 0, 1.0]
+    chord, q = 1.7, 0.8
+    omega = q * y_airf
+    leading = [0.0, 0.0, 0.0]
+    trailing = chord * x_airf
+    @test section_pitch_rate(cross(omega, leading), cross(omega, trailing),
+                             z_airf, chord) ≈ q
+    @test section_pitch_rate([0.0, 0, 0], [0.0, 0, 0], z_airf, chord) == 0.0
+    @test section_pitch_rate([0.0, 0, 0], [0.0, 0, 1.0], z_airf, 0.0) == 0.0
+end
+
+@testset "Flow curvature pitch-rate moment" begin
+    chord, span, V, q = 1.0, 8.0, 20.0, 0.5
+
+    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])
+    set_va!(body_aero, [V, 0.0, 0.0], [0.0, q, 0.0])
+
+    panel = body_aero.panels[6]
+    q_local = dot(body_aero.omega, panel.y_airf)
+
+    @testset "increment matches thin airfoil theory" begin
+        @test body_aero.pitch_rate_dist[6] ≈ q_local
+        @test flow_curvature_cm(q_local, chord, V) ≈
+              -0.25π * q_local * chord / (2V)
+        @test flow_curvature_cm(0.0, chord, V) == 0.0
+        @test flow_curvature_cm(q_local, chord, 0.0) == 0.0
+    end
+
+    @testset "positive rate about y_airf raises aft incidence" begin
+        # nose-up rotation loads the aft chord more, which is what makes the
+        # negative increment a damping rather than a driving moment
+        aft = panel.aero_center + 0.1chord * panel.x_airf
+        dv = cross(body_aero.omega, panel.aero_center) -
+             cross(body_aero.omega, aft)
+        @test sign(dot(dv, panel.z_airf)) == sign(q_local)
+    end
+
+    solver_off = Solver(body_aero; flow_curvature=false)
+    solver_on = Solver(body_aero; flow_curvature=true)
+
+    function moment_at(solver, omega)
+        set_va!(body_aero, [V, 0.0, 0.0], omega)
+        solve!(solver, body_aero)
+        return copy(solver.sol.moment)
+    end
+
+    @testset "opposes the rotation" begin
+        # closed form of the summed increment, pi*rho*V*S*c^2/16 per rad/s,
+        # exact in the limit where the induced velocity is small against V
+        expected = π * 1.225 * V * span * chord * chord^2 / 16 * q
+        for omega in ([0.0, q, 0.0], [0.0, -q, 0.0])
+            dM = moment_at(solver_on, omega) .- moment_at(solver_off, omega)
+            @test dot(dM, omega) < 0
+            @test norm(dM) ≈ expected rtol = 0.05
+        end
+    end
+
+    @testset "no effect without rotation" begin
+        dM = moment_at(solver_on, zeros(3)) .- moment_at(solver_off, zeros(3))
+        @test all(iszero, dM)
+    end
+
+    @testset "defaults to off" begin
+        @test Solver(body_aero).flow_curvature == false
+        @test VortexStepMethod.SolverSettings().flow_curvature == false
+    end
+
+    @testset "distributed rates drive a deformation mode" begin
+        n = length(body_aero.panels)
+        va_dist = repeat([V 0.0 0.0], n)
+
+        set_va!(body_aero, va_dist)
+        @test all(iszero, body_aero.pitch_rate_dist)
+        solve!(solver_on, body_aero)
+        base = copy(solver_on.sol.cm_dist)
+
+        # antisymmetric twist rate: no rigid-body omega can express this
+        rates = [panel.aero_center[2] > 0 ? 1.0 : -1.0 for panel in body_aero.panels]
+        set_va!(body_aero, va_dist; pitch_rate_dist=rates)
+        @test body_aero.pitch_rate_dist ≈ rates
+        solve!(solver_on, body_aero)
+
+        for i in 1:n
+            @test solver_on.sol.cm_dist[i] - base[i] ≈
+                  flow_curvature_cm(rates[i], solver_on.sol._chord_dist[i],
+                                    solver_on.lr.v_a_dist[i])
+        end
+        # the two half-wings must be driven in opposite senses
+        @test sign(solver_on.sol.cm_dist[1] - base[1]) ==
+              -sign(solver_on.sol.cm_dist[n] - base[n])
+
+        @test_throws ArgumentError set_va!(body_aero, va_dist;
+                                           pitch_rate_dist=rates[1:end-1])
+    end
+end