feat(vpm): compute in-plane hub force (H-force) in the free-wake march

AGENTS.md

@@ -166,6 +166,25 @@ into the disk plane by the local flap angle beta(psi) -- and its sign
 convention are documented in the blade-flapping section below, since it
 shares the flap-angle sign convention and the harmonic flap solve.
 
+### VPM in-plane H-force
+
+`VpmRotor` assembles the same in-plane hub force directly in its
+time-domain loads loop (`dynbem_rs/src/vpm/march.rs`), NOT via
+`apply_flap_dynamics`. Per element it accumulates
+`fx += dF_t*sin(psi) - dT*sin(beta_b)*cos(psi)` and
+`fy += dF_t*cos(psi) + dT*sin(beta_b)*sin(psi)`, where `dF_t` is the
+tangential force (same one whose moment is the torque) and `beta_b` is
+the instantaneous per-blade flap angle. The first term is the
+profile/induced-drag H-force (identical pairing to `SweepCtx::run`); the
+second is the flapping-tilt term, resolved instantaneously from the real
+flap DOF (so it needs no b1c/b1s harmonics and is exactly zero for a
+rigid rotor). Note `FlightCondition::v_hub` is the freestream air
+velocity in the hub frame, so a `+X` freestream yields a downwind `+X`
+in-plane force -- matching the level-disk crosswind direction validated
+in `validation_rs/src/checks/h_force.rs`. The VPM direction is validated
+by `validation_rs/src/checks/vpm_hforce_directional.rs`; if you flip any
+of these signs re-run both checks.
+
 ## Pitt-Peters inflow model
 
 Implementation: `dynbem_rs/src/pitt_peters.rs`.

docs/VPM_DESIGN.md

@@ -523,13 +523,26 @@ revolution (AGENTS.md hub-frame convention, thrust along $-Z$):
 ```math
 \begin{aligned}
 dT &= dL\cos\phi - dD\sin\phi \\
+dF_t &= dL\sin\phi + dD\cos\phi \\
 T &= \sum dT \\
-Q &= \sum r\,(dL\sin\phi + dD\cos\phi) \\
+Q &= \sum r\,dF_t \\
 M_x &= \sum r\,dT\,\sin\psi \quad (\text{roll-right positive}) \\
-M_y &= \sum r\,dT\,\cos\psi \quad (\text{pitch-up positive})
+M_y &= \sum r\,dT\,\cos\psi \quad (\text{pitch-up positive}) \\
+F_x &= \sum \big(dF_t\,\sin\psi - dT\,\sin\beta_b\,\cos\psi\big) \\
+F_y &= \sum \big(dF_t\,\cos\psi + dT\,\sin\beta_b\,\sin\psi\big)
 \end{aligned}
 ```
 
+`F_x`, `F_y` are the in-plane hub force (the "H-force"), summed with the same
+convention as the BEM-family models (`bem_common::SweepCtx::run` and
+`apply_flap_dynamics`). Two physical contributions add per element: the
+profile/induced-drag term projects the tangential force `dF_t` onto the fixed
+hub axes, and the flapping-tilt term projects the thrust that leans in-plane
+when the blade is flapped up by `beta_b`. Because the VPM marches the flap DOF
+in the time domain (Section 5.6), the tilt term is resolved instantaneously
+from the real disk geometry -- no separate harmonic flap solve is needed, and
+it is exactly zero for a rigid rotor (`beta_b = 0`).
+
 Blade flapping and servo-flap feathering are modeled as optional per-blade
 time-domain DOFs (Section 5.6). The rigid-blade path remains the default when
 those properties are not supplied.
@@ -753,6 +766,7 @@ vs direct ($<5\%$ of peak at $\theta = 0.5$).
 | `flapping_harmonics` | Bramwell / Seddon flap theory | Coning `a0`, longitudinal flap `a1` vs closed form | `a0`, `a1` within ~14% |
 | `cyclic_sign` | Directional (AGENTS.md) | Collective monotone; cyclic tilt signs | PASS |
 | `flap_directional` | Directional | Flap coning in hover; hub-moment relief under cyclic | PASS |
+| `vpm_hforce_directional` | Directional | In-plane hub force (H-force) downwind in edgewise flow; ~0 in hover | PASS |
 | `servo_flap` | Directional | Kaman servo-flap feathering (zero / collective / cyclic) | PASS |
 | `cyclic_phase_servo` | Directional | Direct-mech pitching `My` vs servo-flap rolling `Mx` | PASS |
 

dynbem_rs/src/vpm/march.rs

@@ -152,6 +152,8 @@ impl<P: Polar> VpmRotor<P> {
         let mut q_acc = 0.0;
         let mut mx_acc = 0.0;
         let mut my_acc = 0.0;
+        let mut fx_acc = 0.0;
+        let mut fy_acc = 0.0;
         let mut avg_count = 0usize;
 
         // Scratch reused each step.
@@ -179,6 +181,10 @@ impl<P: Polar> VpmRotor<P> {
             let mut torque_step = 0.0;
             let mut mx_step = 0.0;
             let mut my_step = 0.0;
+            // In-plane hub force (H-force) accumulators, projected onto the
+            // fixed hub axes each element below.
+            let mut fx_step = 0.0;
+            let mut fy_step = 0.0;
 
             // Aerodynamic flap moment about the hinge, per blade, accumulated
             // in the loads loop below (M = sum r * dF_z). Drives the flap ODE.
@@ -386,11 +392,31 @@ impl<P: Polar> VpmRotor<P> {
                     let dl = q_dyn * c * cl * self.dr[i];
                     let dd = q_dyn * c * cd * self.dr[i];
                     let d_thrust = dl * phi.cos() - dd * phi.sin(); // up (-Z)
+                                                                    // In-plane tangential aero force (the same force whose
+                                                                    // moment gives the torque, opposing blade motion).
+                    let d_ft = dl * phi.sin() + dd * phi.cos();
                     thrust_step += d_thrust;
-                    torque_step += (dl * phi.sin() + dd * phi.cos()) * r;
+                    torque_step += d_ft * r;
                     // Hub moments (AGENTS.md): Mx = r dT sin psi, My = r dT cos psi.
                     mx_step += r * d_thrust * spsi;
                     my_step += r * d_thrust * cpsi;
+                    // In-plane hub force (H-force), two physical contributions
+                    // summed per element (same convention as
+                    // bem_common::SweepCtx::run and the flapping-tilt term in
+                    // apply_flap_dynamics):
+                    //   1. profile/induced-drag: the tangential force d_ft acts
+                    //      along -t_hat = [sin psi, cos psi, 0], so it projects
+                    //      as (d_ft*sin psi, d_ft*cos psi).
+                    //   2. flapping-tilt: when the blade is flapped up by beta_b
+                    //      the normal (thrust) force tilts by beta about t_hat,
+                    //      giving an in-plane component -d_thrust*sin(beta)*r_hat
+                    //      with r_hat = [cos psi, -sin psi, 0]. This term is zero
+                    //      when flap is inactive (beta_b = 0). Unlike the BEM
+                    //      harmonic solve this is resolved instantaneously from
+                    //      the real flap DOF, so it needs no b1c/b1s harmonics.
+                    let tilt = d_thrust * beta_b.sin();
+                    fx_step += d_ft * spsi - tilt * cpsi;
+                    fy_step += d_ft * cpsi + tilt * spsi;
                     // Aero flap moment about the hinge: out-of-plane force
                     // (d_thrust, up) at arm r. Positive -> flaps blade up.
                     m_flap[b] += r * d_thrust;
@@ -518,6 +544,8 @@ impl<P: Polar> VpmRotor<P> {
                 q_acc += torque_step;
                 mx_acc += mx_step;
                 my_acc += my_step;
+                fx_acc += fx_step;
+                fy_acc += fy_step;
                 avg_count += 1;
             }
         }
@@ -529,6 +557,8 @@ impl<P: Polar> VpmRotor<P> {
             torque: q_acc * inv,
             mx_hub: mx_acc * inv,
             my_hub: my_acc * inv,
+            fx_hub: fx_acc * inv,
+            fy_hub: fy_acc * inv,
             n_particles: wake.len(),
             wake_centroid: centroid,
         };

dynbem_rs/src/vpm/mod.rs

@@ -243,6 +243,10 @@ pub struct VpmRotorResult {
     pub mx_hub: f64,
     /// Hub pitching moment (pitch-up positive), N*m.
     pub my_hub: f64,
+    /// In-plane hub force along hub +X (North-ish), N.
+    pub fx_hub: f64,
+    /// In-plane hub force along hub +Y (East-ish), N.
+    pub fy_hub: f64,
     /// Wake particle count at the end of the run.
     pub n_particles: usize,
     /// Mean wake position (for diagnostics, e.g. skew).
@@ -546,16 +550,19 @@ impl<P: Polar> AeroModel for VpmRotor<P> {
         let (fc, kin) = self.flight_condition(inputs);
         // One sub-step per call; dt IS the sub-step duration.
         let (res, out_state) = self.march_window(&fc, Some(state), state.psi, dt, 1, 1);
-        // TODO: VPM's free-wake march doesn't yet accumulate the in-plane
-        // hub force (H-force) the BEM-family models compute in
-        // bem_common::SweepCtx::run -- pass 0.0, 0.0 until that's added here.
+        // The free-wake march accumulates the in-plane hub force (H-force)
+        // directly from the blade-element loads (profile/induced drag plus,
+        // when the flap DOF is active, the instantaneous thrust-vector tilt) --
+        // see the loads loop in `march_window`. No separate harmonic flap solve
+        // is needed here (unlike the BEM-family models) because VPM marches the
+        // flap DOF in the time domain.
         let result = assemble_result(
             res.thrust,
             res.torque,
             res.mx_hub,
             res.my_hub,
-            0.0,
-            0.0,
+            res.fx_hub,
+            res.fy_hub,
             kin.hub_axis,
             &inputs.R_hub,
         );

validation_rs/src/checks/cyclic_phase_servo.rs

@@ -150,6 +150,8 @@ fn march_nan_debug(
         torque: results.iter().map(|r| r.torque).sum::<f64>() / n,
         mx_hub: results.iter().map(|r| r.mx_hub).sum::<f64>() / n,
         my_hub: results.iter().map(|r| r.my_hub).sum::<f64>() / n,
+        fx_hub: results.iter().map(|r| r.fx_hub).sum::<f64>() / n,
+        fy_hub: results.iter().map(|r| r.fy_hub).sum::<f64>() / n,
         n_particles: results.last().map(|r| r.n_particles).unwrap_or(0),
         wake_centroid: results.last().map(|r| r.wake_centroid).unwrap_or([0.0; 3]),
     }

validation_rs/src/checks/mod.rs

@@ -20,6 +20,7 @@ mod servo_flap;
 mod vpm_aging_survey;
 mod vpm_empirical;
 mod vpm_engine_compare;
+mod vpm_hforce_directional;
 mod wake_skew;
 
 pub use autorotation::check_autorotation;
@@ -42,11 +43,12 @@ pub use servo_flap::check_servo_flap;
 pub use vpm_aging_survey::check_vpm_aging_survey;
 pub use vpm_empirical::check_vpm_empirical;
 pub use vpm_engine_compare::check_vpm_engine_compare;
+pub use vpm_hforce_directional::check_vpm_hforce_directional;
 pub use wake_skew::check_wake_skew;
 
 use crate::report::Report;
 
-/// Run all 18 validation checks in order.
+/// Run all 19 validation checks in order.
 pub fn run_all_checks(report: &mut Report) {
     run_filtered_checks(report, None);
 }
@@ -76,6 +78,7 @@ pub fn run_filtered_checks(report: &mut Report, filter: Option<&str>) {
     maybe!("servo_flap", check_servo_flap);
     maybe!("cyclic_phase_servo", check_cyclic_phase_servo);
     maybe!("flap_directional", check_flap_directional);
+    maybe!("vpm_hforce_directional", check_vpm_hforce_directional);
     maybe!("hover_empirical", check_hover_empirical);
     maybe!("hover_cq_empirical", check_hover_cq_empirical);
     maybe!("descent_empirical", check_descent_empirical);

validation_rs/src/checks/vpm_hforce_directional.rs

@@ -0,0 +1,95 @@
+// Directional check: the VPM free-wake rotor accumulates an in-plane hub
+// force ("H-force") from the blade-element loads, matching the BEM-family
+// convention (see checks/h_force.rs for the quasi-static-BEM analogue and
+// dynbem_rs/src/vpm/march.rs for the projection).
+//
+// Physical claim under test: with edgewise flow across a level disk the
+// rotor's in-plane hub force points along the incoming freestream
+// (downwind), the same directional convention as the BEM level-disk
+// crosswind check. FlightCondition::v_hub is the freestream air velocity in
+// the hub frame (relative wind = v_hub - blade velocity), so a +X freestream
+// yields a +X in-plane force. In hover the loading is axisymmetric and the
+// in-plane force vanishes.
+//
+// The VPM marches the flap DOF in the time domain, so the flapping-tilt part
+// of the H-force emerges from the instantaneous disk geometry -- no separate
+// harmonic flap solve is needed (unlike the BEM-family models). This check
+// runs with a rigid rotor to isolate the always-present profile/induced-drag
+// term; the flapping-tilt contribution is exercised in flap_directional.
+
+use crate::helpers::*;
+use crate::report::Report;
+use std::f64::consts::PI;
+
+pub fn check_vpm_hforce_directional(r: &mut Report) {
+    r.begin_module(
+        "vpm_hforce_directional",
+        "Directional: VPM in-plane hub force (H-force) points downwind in edgewise flow; ~0 in hover",
+    );
+
+    let dt = (2.0 * PI / OMEGA) / STEPS_PER_REV as f64;
+    let n_revs = 21 * STEPS_PER_REV;
+
+    let defn = theory_rotor(10, 0.0);
+    let rotor = make_fast_rotor(&defn);
+
+    // --- Hover baseline: axisymmetric loading -> no in-plane force ---
+    let (hover, _) = rotor.march(&hover_fc(8.0), None, dt, n_revs);
+    let h_hover = hover.fx_hub.hypot(hover.fy_hub);
+    r.info("hover", "fx_hub_N", hover.fx_hub, f64::NAN);
+    r.info("hover", "fy_hub_N", hover.fy_hub, f64::NAN);
+    r.info("hover", "h_force_N", h_hover, f64::NAN);
+    r.info("hover", "thrust_N", hover.thrust, f64::NAN);
+    r.assert_bool(
+        "hover",
+        "no_in_plane_force",
+        h_hover,
+        0.0,
+        h_hover < 0.02 * hover.thrust.abs().max(1.0),
+        &format!(
+            "hover loading is axisymmetric -> in-plane force ~0, got |H|={:.4} N (T={:.2} N)",
+            h_hover, hover.thrust
+        ),
+    );
+
+    // --- Edgewise flow along +X: H-force points downwind (+X) ---
+    let (fwd, _) = rotor.march(&forward_fc(8.0, 0.15), None, dt, n_revs);
+    r.info("edgewise_flow", "fx_hub_N", fwd.fx_hub, f64::NAN);
+    r.info("edgewise_flow", "fy_hub_N", fwd.fy_hub, f64::NAN);
+    r.info("edgewise_flow", "thrust_N", fwd.thrust, f64::NAN);
+    r.assert_bool(
+        "edgewise_flow",
+        "h_force_downwind",
+        fwd.fx_hub,
+        0.0,
+        fwd.fx_hub > 0.05,
+        &format!(
+            "+X freestream should give a downwind (+X) H-force, got Fx_hub={:.4} N",
+            fwd.fx_hub
+        ),
+    );
+    r.assert_bool(
+        "edgewise_flow",
+        "mostly_aligned_with_flow",
+        fwd.fx_hub.abs(),
+        fwd.fy_hub.abs(),
+        fwd.fx_hub.abs() > fwd.fy_hub.abs(),
+        &format!(
+            "H-force should be mostly along the flow axis: |Fx|={:.4} N vs |Fy|={:.4} N",
+            fwd.fx_hub.abs(),
+            fwd.fy_hub.abs()
+        ),
+    );
+    r.assert_bool(
+        "edgewise_flow",
+        "larger_than_hover",
+        fwd.fx_hub.hypot(fwd.fy_hub),
+        h_hover,
+        fwd.fx_hub.hypot(fwd.fy_hub) > h_hover,
+        &format!(
+            "edgewise flow should raise the H-force above the hover baseline: {:.4} N vs {:.4} N",
+            fwd.fx_hub.hypot(fwd.fy_hub),
+            h_hover
+        ),
+    );
+}