diff --git a/AGENTS.md b/AGENTS.md index 2679e26..ce3f51f 100644 --- a/AGENTS.md +++ b/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`. diff --git a/docs/VPM_DESIGN.md b/docs/VPM_DESIGN.md index b9efb5d..503ee5f 100644 --- a/docs/VPM_DESIGN.md +++ b/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 | diff --git a/dynbem_rs/src/vpm/march.rs b/dynbem_rs/src/vpm/march.rs index dbb3311..621893d 100644 --- a/dynbem_rs/src/vpm/march.rs +++ b/dynbem_rs/src/vpm/march.rs @@ -152,6 +152,8 @@ impl VpmRotor

{ 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 VpmRotor

{ 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 VpmRotor

{ 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 VpmRotor

{ 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 VpmRotor

{ 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, }; diff --git a/dynbem_rs/src/vpm/mod.rs b/dynbem_rs/src/vpm/mod.rs index 78a263f..bf04279 100644 --- a/dynbem_rs/src/vpm/mod.rs +++ b/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 AeroModel for VpmRotor

{ 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, ); diff --git a/validation_rs/src/checks/cyclic_phase_servo.rs b/validation_rs/src/checks/cyclic_phase_servo.rs index cb5c8e0..c0d02e7 100644 --- a/validation_rs/src/checks/cyclic_phase_servo.rs +++ b/validation_rs/src/checks/cyclic_phase_servo.rs @@ -150,6 +150,8 @@ fn march_nan_debug( torque: results.iter().map(|r| r.torque).sum::() / n, mx_hub: results.iter().map(|r| r.mx_hub).sum::() / n, my_hub: results.iter().map(|r| r.my_hub).sum::() / n, + fx_hub: results.iter().map(|r| r.fx_hub).sum::() / n, + fy_hub: results.iter().map(|r| r.fy_hub).sum::() / 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]), } diff --git a/validation_rs/src/checks/mod.rs b/validation_rs/src/checks/mod.rs index 51fc31a..f821b1a 100644 --- a/validation_rs/src/checks/mod.rs +++ b/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); diff --git a/validation_rs/src/checks/vpm_hforce_directional.rs b/validation_rs/src/checks/vpm_hforce_directional.rs new file mode 100644 index 0000000..ff57f59 --- /dev/null +++ b/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 + ), + ); +}