From 030ee892c8f72016cb45a5ab0cf9ff0837cd14f6 Mon Sep 17 00:00:00 2001
From: Bart <bart@vandelint.net>
Date: Mon, 17 Aug 2026 10:48:44 +0200
Subject: [PATCH 1/2] Optional flow-curvature pitch-rate moment

A section rotating about its own spanwise axis sees an incidence that varies
linearly along the chord. In thin airfoil theory that is the same boundary
condition as parabolic camber, giving A1 = q*c/(2V) and A2 = 0, hence
cm_c/4 = (pi/4)(A2 - A1) = -(pi/4) q_hat.

The lift half of the same effect was already exact: A0 + A1/2 collapses to the
effective incidence at three-quarter chord, which is where the inflow is
sampled, so no correction belongs there. Only the moment was missing, and one
control point per section cannot represent it because it is an A1 term rather
than an A0 shift.

Off by default. Enabled through the `flow_curvature` solver setting, applied in
both force paths (`calc_forces!` behind `solve!`, and `calculate_results`
behind `solve`). Independent of the pivot, so it needs the rate only and not a
rotation centre.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 CHANGELOG.md                       | 14 +++++++
 docs/src/private_functions.md      |  1 +
 src/body_aerodynamics.jl           | 26 +++++++++++-
 src/settings.jl                    |  3 ++
 src/solver.jl                      | 12 +++++-
 test/runtests.jl                   |  1 +
 test/solver/test_flow_curvature.jl | 65 ++++++++++++++++++++++++++++++
 7 files changed, 120 insertions(+), 2 deletions(-)
 create mode 100644 test/solver/test_flow_curvature.jl

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 50f48fdf..72deecf1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,19 @@
 # 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. Read from `body_aero.omega`, so it is only meaningful
+  after `set_va!(body_aero, va, omega)`; the matrix form of `set_va!` does not
+  set `omega`.
+
 ## VortexStepMethod v4.0.0 2026-08-03
 
 ### Added
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index bec6cb73..d678cded 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -17,6 +17,7 @@ calculate_cl
 calculate_cd
 calculate_cm
 calculate_cd_cm
+flow_curvature_cm
 calculate_relative_alpha_and_velocity
 calculate_relative_alpha_and_relative_velocity
 update_effective_angle_of_attack!
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index 044d4df7..360bf6cf 100644
--- a/src/body_aerodynamics.jl
+++ b/src/body_aerodynamics.jl
@@ -714,7 +714,23 @@ function compute_panel_center_of_pressures(
 end
 
 """
-    calculate_results(body_aero::BodyAerodynamics, gamma_new, 
+    flow_curvature_cm(omega, panel, 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 = ω ⋅ y_airf` 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(omega, panel, chord, v_rel)
+    v_rel > 0 || return zero(chord)
+    return -0.25π * dot3(omega, panel.y_airf) * chord / (2v_rel)
+end
+
+"""
+    calculate_results(body_aero::BodyAerodynamics, gamma_new,
                      density,
                      core_radius_fraction, mu,
                      alpha_dist, v_a_dist,
@@ -726,6 +742,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 +766,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 +795,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.omega, panel, 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)
diff --git a/src/settings.jl b/src/settings.jl
index 83b05c90..ebc65ed5 100644
--- a/src/settings.jl
+++ b/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
 
 """
diff --git a/src/solver.jl b/src/solver.jl
index b3e36d03..672c2243 100644
--- a/src/solver.jl
+++ b/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,10 @@ 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.omega, panel, 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 +607,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 +1227,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
diff --git a/test/runtests.jl b/test/runtests.jl
index 6635b965..d613c2b5 100644
--- a/test/runtests.jl
+++ b/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")
diff --git a/test/solver/test_flow_curvature.jl b/test/solver/test_flow_curvature.jl
new file mode 100644
index 00000000..bffd1318
--- /dev/null
+++ b/test/solver/test_flow_curvature.jl
@@ -0,0 +1,65 @@
+using VortexStepMethod
+using VortexStepMethod: flow_curvature_cm
+using LinearAlgebra
+using Test
+
+@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 flow_curvature_cm(body_aero.omega, panel, chord, V) ≈
+              -0.25π * q_local * chord / (2V)
+        @test flow_curvature_cm(zeros(3), panel, chord, V) == 0.0
+        @test flow_curvature_cm(body_aero.omega, panel, 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
+end

From 215ba54b0c104a5af96de2ab0c1e4d231426b364 Mon Sep 17 00:00:00 2001
From: Bart <bart@vandelint.net>
Date: Mon, 17 Aug 2026 11:02:07 +0200
Subject: [PATCH 2/2] Per-panel pitch rates for flow curvature

The moment increment responds to the chordwise gradient of normal velocity
across a panel. For rigid motion that gradient is omega . y_airf, but a
deforming wing has a twist rate that varies along the span and no single body
rate expresses it, so flapping and torsion modes could not reach the term.

BodyAerodynamics now carries pitch_rate_dist, one rate per panel. The rigid
set_va! fills it by projecting omega onto each y_airf, which already differs per
panel under dihedral. The distributed set_va! takes it as a keyword, so a
structural model can pass rates it computed itself.

section_pitch_rate builds one entry from a section's leading and trailing edge
velocities. It reduces to omega . y_airf for rigid motion, so one expression
covers both cases and the sign convention is pinned in one place rather than
rediscovered by every caller.

Omitting the keyword zeroes the rates. The distributed set_va! never set omega,
so reading it there would have fed the moment a stale rate from whenever the
rigid form was last called.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
---
 CHANGELOG.md                       | 15 ++++--
 docs/src/functions.md              |  1 +
 docs/src/private_functions.md      |  1 +
 src/VortexStepMethod.jl            |  2 +-
 src/body_aerodynamics.jl           | 75 ++++++++++++++++++++++++++----
 src/solver.jl                      |  3 +-
 test/solver/test_flow_curvature.jl | 49 +++++++++++++++++--
 7 files changed, 129 insertions(+), 17 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 72deecf1..d593e4a7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,9 +10,18 @@
   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. Read from `body_aero.omega`, so it is only meaningful
-  after `set_va!(body_aero, va, omega)`; the matrix form of `set_va!` does not
-  set `omega`.
+  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
 
diff --git a/docs/src/functions.md b/docs/src/functions.md
index 72de8257..b5cc0dd4 100644
--- a/docs/src/functions.md
+++ b/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}
diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md
index d678cded..c6ba50f9 100644
--- a/docs/src/private_functions.md
+++ b/docs/src/private_functions.md
@@ -18,6 +18,7 @@ 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!
diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl
index c8b8de59..7f707116 100644
--- a/src/VortexStepMethod.jl
+++ b/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
 
diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl
index 360bf6cf..96e4c60b 100644
--- a/src/body_aerodynamics.jl
+++ b/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,19 +718,49 @@ function compute_panel_center_of_pressures(
 end
 
 """
-    flow_curvature_cm(omega, panel, chord, v_rel)
+    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 = ω ⋅ y_airf` 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.
+`Δ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(omega, panel, chord, v_rel)
+@inline function flow_curvature_cm(pitch_rate, chord, v_rel)
     v_rel > 0 || return zero(chord)
-    return -0.25π * dot3(omega, panel.y_airf) * chord / (2v_rel)
+    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
 
 """
@@ -797,7 +831,7 @@ function calculate_results(
             panel, alpha_dist[i])
         if flow_curvature
             cm_array[i] += flow_curvature_cm(
-                body_aero.omega, panel, chord_array[i], v_a_dist[i])
+                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]
@@ -1082,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)
@@ -1116,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, :]
diff --git a/src/solver.jl b/src/solver.jl
index 672c2243..ba31e1fd 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -294,7 +294,8 @@ function calc_forces!(solver::Solver{P, U, T}, body_aero::BodyAerodynamics;
         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.omega, panel, solver.sol._chord_dist[i], v_a_dist[i])
+                body_aero.pitch_rate_dist[i], solver.sol._chord_dist[i],
+                v_a_dist[i])
         end
         width_dist[i] = panel.width
 
diff --git a/test/solver/test_flow_curvature.jl b/test/solver/test_flow_curvature.jl
index bffd1318..7a5c6b08 100644
--- a/test/solver/test_flow_curvature.jl
+++ b/test/solver/test_flow_curvature.jl
@@ -3,6 +3,20 @@ 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
 
@@ -18,10 +32,11 @@ using Test
     q_local = dot(body_aero.omega, panel.y_airf)
 
     @testset "increment matches thin airfoil theory" begin
-        @test flow_curvature_cm(body_aero.omega, panel, chord, V) ≈
+        @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(zeros(3), panel, chord, V) == 0.0
-        @test flow_curvature_cm(body_aero.omega, panel, chord, 0.0) == 0.0
+        @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
@@ -62,4 +77,32 @@ using Test
         @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