Adds an opt-in flow_curvature solver setting (default false) giving each
section the thin-airfoil pitch-rate moment increment, driven by a per-panel rate
so deformation modes reach it.

What is missing today

A section rotating about its own spanwise axis at rate q sees an incidence that
varies linearly along the chord, α_eff(x) = α + q(x − x_p)/V. Thin airfoil theory
cannot tell that apart from parabolic camber, so the Fourier coefficients pick up

A₀ = α + (q/V)(c/2 − x_p),   A₁ = (q/V)(c/2),   A₂ = 0

The lift half is already exact in VSM: Cl = 2π(A₀ + A₁/2) = 2π[α + (q/V)(3c/4 − x_p)]
is just the effective incidence at three-quarter chord — precisely where set_va!
samples the inflow (control_points = LE·0.25 + TE·0.75). Nothing to fix there.

The moment is not, because it lives in A₁ rather than A₀:

Cm_c/4 = (π/4)(A₂ − A₁) = −(π/4) q̂,    q̂ = q c / (2 v_rel)

One control point per section can represent an incidence but not an incidence
gradient. A chordwise-refined VLM recovers this from the discretisation; a
lifting line has to add it analytically. It is independent of the pivot (x_p
only enters A₀), so only the rate is needed — the arm effect is already handled
by the existing (aero_center − reference_point) × F sum.

Per-panel rates

What the moment responds to is the chordwise gradient of normal velocity across a
panel. For rigid motion that gradient is ω ⋅ y_airf, but a deforming wing has a
twist rate that varies along the span and no single body rate expresses it — which
matters because this moment is largest exactly there (a strip's torsional inertia
is small while the moment scales with c²).

So the input is the gradient, not a body rate:

• pitch_rate_dist on BodyAerodynamics, one rate per panel, positive nose-up.
• set_va!(body_aero, va, omega) fills it by projecting onto each y_airf, so
  panels at different dihedral see different rates from one body rate.
• set_va!(body_aero, va_distribution; pitch_rate_dist) takes it directly.
• section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord),
  exported, builds one entry from a section's edge velocities.

section_pitch_rate reduces exactly to ω ⋅ y_airf under rigid motion
(v_TE − v_LE = ω × c·x_airf, and y × x = −z), so one expression covers both
cases and the sign convention — where an error flips damping into divergence — is
pinned in one place instead of being rediscovered by each caller.

It also closes a trap: the distributed set_va! never set omega, so reading it
there would have used a stale rate. Omitting the keyword now zeroes it.

Verification

• The panel frame is right-handed and q > 0 is nose-up: a station 0.1c aft of
  the aero centre gains a positive z_airf velocity component, i.e. more
  incidence, matching the derivation.
• dM · ω < 0 for both signs of ω — it damps rather than drives.
• The summed increment matches the closed form πρVSc²/16 per rad/s
  (predicted 19.24 N·m at the test condition, measured 19.245).
• Exactly zero at ω = 0, so default-off results are untouched.
• Distributed path: an antisymmetric twist rate no rigid ω can produce, per-panel
  cm_dist increments matching flow_curvature_cm, opposite signs on the two
  half-wings, length mismatch throwing.
• The relation itself is checked against measured section data, not just against
  the code. A rotating section is the A₁ = q c/(2 v_rel), A₂ = 0 case of
  cm_c/4 = (π/4)(A₂ − A₁); applying the same relation to NACA 4-digit mean
  lines gives −0.053 for the 2412 against a measured −0.047, −0.106 for the
  4412 against −0.093, and exactly zero for the symmetric 0012. The zero-lift
  angle is asserted from the same coefficients (−2.08° against a measured
  −2.1°), so a sign or factor error cannot hide in the moment alone.

Applied in both force paths — calc_forces! behind solve! and
calculate_results behind solve.

Caveats

• Quasi-steady. Roughly a third of the full sectional pitch-damping moment; the
  apparent-mass part (−(π/2) q̂) is larger and belongs in the dynamics layer.
• Second commit regenerates both default manifests: Arrow had become a direct
  dependency without reaching either, so seeding a fresh checkout from them left
  it unloadable. Unrelated to the physics change.

References

The relation this rests on, and the coefficient definitions used above:

MIT OpenCourseWare 16.01 Unified Engineering I–IV, Fall 2005 – Spring 2006,
Fluids Lecture 3 Notes: Thin-Airfoil Analysis Problem
(OCW).
A₀ = α − (1/π)∫₀^π (dZ/dx) dθ and Aₙ = (2/π)∫₀^π (dZ/dx) cos nθ dθ on p. 2;
c_ℓ = π(2A₀ + A₁) and c_m,c/4 = (π/4)(A₂ − A₁) on p. 3. The notes give
Anderson, Fundamentals of Aerodynamics §4.8 as the reading.

Only the substitution of the rotation-induced camber slope into those definitions
is done here; the integrals are in the description above and the resulting
relation is exercised against measured section data in the test suite.

For the apparent-mass term named in the caveats, which this PR does not add:

Theodorsen, T., General Theory of Aerodynamic Instability and the Mechanism of
Flutter, NACA Report No. 496, 1935.
Section "Velocity potentials, forces, and moments of the noncirculatory flow",
pp. 292–294; the noncirculatory force and moment are eqs. (I) and (III) on
p. 294, with M_α taken about the rotation axis at ab from midchord and
lengths non-dimensionalised on the semichord b.
Scan: https://ntrs.nasa.gov/api/citations/19800006788/downloads/19800006788.pdf

Eq. (III) carries the π(1/8 + a²)b²α̈ and baπḧ terms, i.e. the added moment of
inertia about the rotation axis and the added-mass force acting at midchord —
the −(π/2) q̂ moment and the added-mass matrix discussed above both follow from
it. Note NTRS catalogues the scan under 1949, the year of the Annual Report
volume it was reprinted in; the report itself is 1935.

🤖 Generated with Claude Code