diff --git a/Cargo.toml b/Cargo.toml index e643268..30e668f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,7 +4,7 @@ members = ["dynbem", "dynbem_rs", "validation_rs"] # Fields below are inherited by member crates via `field.workspace = true`. [workspace.package] -version = "0.6.0" +version = "0.7.0" edition = "2021" license = "MIT" authors = ["Kristof"] diff --git a/README.md b/README.md index 90358dd..9718c5c 100644 --- a/README.md +++ b/README.md @@ -90,6 +90,7 @@ model = dynbem.create_aero(defn, model="pitt_peters") # or "oye", "bem", "vpm" state = model.initial_rotor_state() omega = 125.7 # rad/s -- caller owns mechanical state +spin_angle = 0.0 # rad -- caller owns rotor azimuth, if tracked inputs = dynbem.RotorInputs( collective_rad=0.14, tilt_lon=0.0, tilt_lat=0.0, # swashplate (helicopter-standard signs) @@ -104,9 +105,11 @@ result, state = model.step(inputs, state, dt, integration_method="semi_implicit" # integration_method: "semi_implicit" | "explicit" | "exponential" (ignored for quasi static and VPM) # result.F_world, result.m_hub_world, result.M_spin, result.Q_spin -# Mechanical ODE lives in the caller: -from dynbem.mechanical import omega_derivative -omega += dt * omega_derivative(result.Q_spin, motor_torque_Nm, I_ode_kgm2) +# Mechanical ODE lives in the caller. Coulomb bearing friction is always +# a parameter (default 0.0). step_omega is the single canonical, recommended +# integrator (semi-implicit, unconditionally stable in the aero term). +from dynbem.mechanical import step_omega +omega, spin_angle = step_omega(omega, spin_angle, result.Q_spin, motor_torque_Nm, I_ode_kgm2, dt) # Level-3 VPM: same step() API; the free-wake particle cloud advances each call. # dt typically T_rev/36-72; average over several revolutions for steady-state loads. diff --git a/docs/API.md b/docs/API.md index 860571d..0265e00 100644 --- a/docs/API.md +++ b/docs/API.md @@ -7,7 +7,7 @@ import dynbem from dynbem import ( create_aero, RotorInputs, AeroResult, solve_trim_cyclic, relax_inflow, TrimResult, - omega_derivative, euler_step_omega, + omega_derivative, step_omega, QuasiStaticRotorState, PittPetersRotorState, OyeRotorState, LinearPolar, TabulatedPolar, ) @@ -76,8 +76,12 @@ result, dstate = model.compute_forces(inputs, state) arr = state.to_array() + dt * dstate.to_array() state = state.from_array(arr) -# Advance rotor speed (caller owns the mechanical ODE) -omega += dt * omega_derivative(result.Q_spin, motor_torque_Nm, I_kgm2) +# Advance rotor speed (caller owns the mechanical ODE). bearing_friction_Nm +# is always a parameter (Coulomb friction); pass 0.0 for a frictionless +# bearing. step_omega is the single canonical, recommended integrator +# (semi-implicit, unconditionally stable in the aero term). spin_angle is +# caller-owned state, advanced alongside omega. +omega, spin_angle = step_omega(omega, spin_angle, result.Q_spin, motor_torque_Nm, I_kgm2, dt) inputs = RotorInputs(..., omega_rad_s=omega, ...) ``` @@ -109,7 +113,7 @@ Returned by `compute_forces`. |---|---|---| | `F_world` | ndarray (3,) | Total rotor force in NED world frame [N]. `F_world[2] < 0` for upward lift | | `m_hub_world` | ndarray (3,) | Hub moments in NED world frame [N*m]. Roll about X, pitch about Y | -| `Q_spin` | float | Aerodynamic reaction torque on the shaft [N*m]. Positive opposes rotation. Feed directly to `omega_derivative` | +| `Q_spin` | float | Aerodynamic reaction torque on the shaft [N*m]. Positive opposes rotation. Feed directly to `omega_derivative`/`step_omega` as `Q_aero` | | `M_spin` | ndarray (3,) | Full spin moment vector in NED world frame [N*m] | --- @@ -146,14 +150,28 @@ model). `to_array()` concatenates `[W_int..., W...]` (length ## Mechanical ODE helpers -### `omega_derivative(Q_aero, motor_torque_Nm, I_ode_kgm2) -> float` +This is the single canonical way to advance the mechanical (rigid-body +spin) ODE -- exactly two functions. Bearing (Coulomb/dry) friction is +always a parameter (`bearing_friction_Nm`, default `0.0`) -- a +constant-magnitude torque that opposes the current rotation direction, +and is exactly zero at `omega == 0` (no static breakaway threshold +modelled). `step_omega` applies a zero-crossing clamp so friction (or any +purely-decelerating explicit term) cannot flip the sign of `omega` within +one step -- the rotor comes to rest instead of reversing direction. -Return `d(omega)/dt = (motor_torque - Q_aero) / I` [rad/s^2]. -Pass `AeroResult.Q_spin` directly for `Q_aero`. +### `omega_derivative(omega, Q_aero, motor_torque_Nm, I_ode_kgm2, bearing_friction_Nm=0.0) -> float` -### `euler_step_omega(omega, spin_angle, Q_aero, motor_torque_Nm, I_ode_kgm2, dt) -> (omega_new, spin_angle_new)` +Return `d(omega)/dt = (motor_torque - Q_aero - Q_friction(omega)) / I` +[rad/s^2]. Pass `AeroResult.Q_spin` directly for `Q_aero`. Use this for +diagnostics or a custom integrator -- do not re-derive the formula. -Forward-Euler step for both `omega` and the spin angle. +### `step_omega(omega, spin_angle, Q_aero, motor_torque_Nm, I_ode_kgm2, dt, bearing_friction_Nm=0.0) -> (omega_new, spin_angle_new)` + +Semi-implicit (locally-frozen relaxation) step -- the one recommended +integrator for both `omega` and the spin angle. Unconditionally stable in +the aerodynamic term: pure aero drag can never send `omega` past zero or +oscillate for any `dt > 0`. Use this for all production and test +time-stepping. --- @@ -266,6 +284,7 @@ model = dynbem.create_aero(defn, model="pitt_peters") state = model.initial_rotor_state() omega = 28.0 # rad/s +spin_angle = 0.0 # rad dt = 0.005 # s I = 0.08 # kg*m^2 (rotor inertia) @@ -287,8 +306,8 @@ for step in range(500): arr = state.to_array() + dt * dstate.to_array() state = state.from_array(arr) - # advance rotor speed (free-spin, no motor torque) - omega += dt * dynbem.omega_derivative(result.Q_spin, 0.0, I) + # advance rotor speed (free-spin, no motor torque, no bearing friction) + omega, spin_angle = dynbem.step_omega(omega, spin_angle, result.Q_spin, 0.0, I, dt) inputs = dynbem.RotorInputs( collective_rad=0.14, tilt_lon=0.0, tilt_lat=0.0, R_hub=np.eye(3), v_hub_world=np.zeros(3), wind_world=np.zeros(3), diff --git a/dynbem/README.md b/dynbem/README.md index 66f25a9..b61661a 100644 --- a/dynbem/README.md +++ b/dynbem/README.md @@ -23,7 +23,7 @@ and Pitt-Peters L-matrix see [`../AGENTS.md`](../AGENTS.md). +-- bem.py compat shim re-exporting dynbem.bem.* +-- cyclic.py compat shim +-- factory.py create_aero(), build_polar(), load_tabulated_polar() - +-- mechanical.py omega_derivative(), euler_step_omega() -- caller-owned spin ODE + +-- mechanical.py omega_derivative(), step_omega() -- single canonical spin ODE (caller-owned), Coulomb bearing friction always a parameter (default 0.0) +-- oye.py compat shim +-- pitt_peters.py compat shim +-- polar.py compat shim + AirfoilPolar tuple alias @@ -48,10 +48,14 @@ loading). Implemented as Python subclasses of the Rust pyclasses in `__init__.py`. Requires `subclass=True` on the three Rust model pyclasses in `wrappers.rs`. -**Mechanical ODE**: `mechanical.py` provides `omega_derivative(Q_aero, -motor_torque_Nm, I_ode_kgm2)` and `euler_step_omega`. The aero models are -pure aerodynamic -- they take `omega_rad_s` via `RotorInputs` and return -loads only. The caller owns the spin ODE. +**Mechanical ODE**: `mechanical.py` provides the single canonical way to +advance the spin ODE -- `omega_derivative(omega, Q_aero, motor_torque_Nm, +I_ode_kgm2, bearing_friction_Nm=0.0)` (raw right-hand side) and +`step_omega` (the one recommended integrator: semi-implicit, unconditionally +stable in the aerodynamic term). Coulomb (dry) bearing friction is always a +parameter, defaulting to `0.0`. The aero +models are pure aerodynamic -- they take `omega_rad_s` via `RotorInputs` +and return loads only. The caller owns the spin ODE. **Virtual ABCs** `dynbem.AeroBase` and `dynbem.RotorState` are declared in `__init__.py` with `ABC.register(...)` so `isinstance(model, AeroBase)` diff --git a/dynbem/pyproject.toml b/dynbem/pyproject.toml index 4c94fd3..7d9c539 100644 --- a/dynbem/pyproject.toml +++ b/dynbem/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "dynbem" -version = "0.6.0" +version = "0.7.0" description = "Rotor aerodynamics: BEM, Pitt-Peters, Oye dynamic inflow, and VPM free-wake -- all in one library" requires-python = ">=3.10" license = "MIT" diff --git a/dynbem/python/dynbem/__init__.py b/dynbem/python/dynbem/__init__.py index e540622..95743a7 100644 --- a/dynbem/python/dynbem/__init__.py +++ b/dynbem/python/dynbem/__init__.py @@ -34,7 +34,7 @@ ) from .factory import create_aero, build_polar, load_tabulated_polar, _build_polar_from_defn # noqa: F401 from .trim import solve_trim_cyclic, relax_inflow # noqa: F401 -from .mechanical import omega_derivative, euler_step_omega # noqa: F401 +from .mechanical import omega_derivative, step_omega # noqa: F401 def cyclic_coeffs(tilt_lon, tilt_lat, control=None): # noqa: F401 diff --git a/dynbem/python/dynbem/mechanical.py b/dynbem/python/dynbem/mechanical.py index 29097ac..7f6123c 100644 --- a/dynbem/python/dynbem/mechanical.py +++ b/dynbem/python/dynbem/mechanical.py @@ -4,28 +4,57 @@ pure aerodynamic: they take omega as part of RotorInputs and return only inflow state derivatives. The caller owns the mechanical ODE. -These helpers implement the standard rigid-body spin-up equation: +These helpers implement the standard rigid-body spin-up equation, with +Coulomb (dry) bearing friction always a first-class parameter (pass +``bearing_friction_Nm=0.0`` for a frictionless bearing): - I * d(omega)/dt = motor_torque - Q_aero + I * d(omega)/dt = motor_torque - Q_aero - Q_friction(omega) d(spin_angle)/dt = omega +The math is implemented in Rust (``dynbem_rs::mechanical``); this module +is a thin wrapper that also carries the ``spin_angle`` bookkeeping. + +This is the single canonical way to advance the mechanical ODE -- there +are exactly two functions: ``omega_derivative`` (the raw ODE right-hand +side, for diagnostics or custom integrators) and ``step_omega`` (the one +recommended integrator, semi-implicit and unconditionally stable in the +aerodynamic term: it never overshoots omega past zero from pure aero +drag). ``step_omega`` also applies a zero-crossing clamp so that Coulomb +friction (or any other purely-decelerating explicit term) cannot flip the +sign of omega within a single step -- the rotor comes to rest instead of +reversing direction. + Usage example:: result, inflow_deriv = aero.compute_forces(inputs, state) - d_omega = omega_derivative(result.Q_spin, motor_torque_Nm, I_ode_kgm2) - omega += dt * d_omega - spin_angle += dt * omega + omega, spin_angle = step_omega( + omega, spin_angle, result.Q_spin, motor_torque_Nm, I_ode_kgm2, dt, + ) inputs.omega_rad_s = omega """ -__all__ = ["omega_derivative", "euler_step_omega"] +from ._dynbem import ( + omega_derivative as _omega_derivative_rust, + step_omega as _step_omega_rust, +) + +__all__ = ["omega_derivative", "step_omega"] -def omega_derivative(Q_aero: float, motor_torque_Nm: float, I_ode_kgm2: float) -> float: +def omega_derivative( + omega: float, + Q_aero: float, + motor_torque_Nm: float, + I_ode_kgm2: float, + bearing_friction_Nm: float = 0.0, +) -> float: """Return d(omega)/dt for the rotor rigid-body spin ODE. Parameters ---------- + omega: + Current rotor speed [rad/s] (needed only to get the sign of the + bearing friction torque right). Q_aero: Aerodynamic reaction torque on the rotor shaft [N.m]. Positive Q_aero opposes rotation (drag convention): use @@ -35,24 +64,41 @@ def omega_derivative(Q_aero: float, motor_torque_Nm: float, I_ode_kgm2: float) - Positive drives rotation in the direction of omega. I_ode_kgm2: Rotor polar moment of inertia about the spin axis [kg.m^2]. + bearing_friction_Nm: + Coulomb (dry) bearing friction torque magnitude [N.m]. Always + opposes the current rotation direction; exactly zero at + ``omega == 0`` (no static breakaway threshold modelled). Pass + ``0.0`` for a frictionless bearing. Returns ------- float d(omega)/dt [rad/s^2]. """ - return (motor_torque_Nm - Q_aero) / I_ode_kgm2 + return _omega_derivative_rust(omega, Q_aero, motor_torque_Nm, I_ode_kgm2, bearing_friction_Nm) -def euler_step_omega( +def step_omega( omega: float, spin_angle: float, Q_aero: float, motor_torque_Nm: float, I_ode_kgm2: float, dt: float, + bearing_friction_Nm: float = 0.0, ) -> tuple: - """Forward-Euler step for the rigid-body spin ODE. + """Semi-implicit (locally-frozen relaxation) step for the spin ODE. + + This is the single canonical integrator for the mechanical ODE -- use + it for all production and test time-stepping. Freezes the local + aerodynamic damping coefficient over the step (the same "freeze the + coefficient/target over dt" idea used for inflow states' semi-implicit + integration) and treats it implicitly, so it is unconditionally stable + in the aerodynamic term: pure aero drag can never send omega past zero + or oscillate for any ``dt > 0``. Motor torque and Coulomb friction are + treated explicitly and are subject to a zero-crossing clamp (the rotor + comes to rest, it doesn't reverse direction). See + ``dynbem_rs::mechanical::step_omega`` for the full derivation. Parameters ---------- @@ -68,12 +114,16 @@ def euler_step_omega( Rotor polar moment of inertia [kg.m^2]. dt: Timestep [s]. + bearing_friction_Nm: + Coulomb (dry) bearing friction torque magnitude [N.m]. Pass + ``0.0`` for a frictionless bearing. Returns ------- (omega_new, spin_angle_new) : (float, float) """ - d_omega = omega_derivative(Q_aero, motor_torque_Nm, I_ode_kgm2) - omega_new = omega + dt * d_omega + omega_new = _step_omega_rust( + omega, Q_aero, motor_torque_Nm, I_ode_kgm2, dt, bearing_friction_Nm + ) spin_angle_new = spin_angle + dt * omega return omega_new, spin_angle_new diff --git a/dynbem/src/lib.rs b/dynbem/src/lib.rs index 3cf73b9..69b74c4 100644 --- a/dynbem/src/lib.rs +++ b/dynbem/src/lib.rs @@ -30,9 +30,55 @@ fn cyclic_coeffs(tilt_lon: f64, tilt_lat: f64, control: Option f64 { + dynbem_rs::mechanical::omega_derivative( + omega, + q_aero, + motor_torque_nm, + i_ode_kgm2, + bearing_friction_nm, + ) +} + +/// Semi-implicit (locally-frozen relaxation) step for the rigid-body spin +/// ODE. This is the single canonical, recommended integrator -- see +/// dynbem_rs::mechanical for the full derivation and stability properties. +#[pyfunction] +#[pyo3(signature = (omega, q_aero, motor_torque_nm, i_ode_kgm2, dt, bearing_friction_nm = 0.0))] +fn step_omega( + omega: f64, + q_aero: f64, + motor_torque_nm: f64, + i_ode_kgm2: f64, + dt: f64, + bearing_friction_nm: f64, +) -> f64 { + dynbem_rs::mechanical::step_omega( + omega, + q_aero, + motor_torque_nm, + i_ode_kgm2, + dt, + bearing_friction_nm, + ) +} + #[pymodule] fn _dynbem(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(cyclic_coeffs, m)?)?; + m.add_function(wrap_pyfunction!(omega_derivative, m)?)?; + m.add_function(wrap_pyfunction!(step_omega, m)?)?; m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/dynbem_rs/README.md b/dynbem_rs/README.md index 16a09e6..c7c3f72 100644 --- a/dynbem_rs/README.md +++ b/dynbem_rs/README.md @@ -115,7 +115,7 @@ Add to `Cargo.toml`: ```toml [dependencies] -dynbem_rs = "0.5" +dynbem_rs = "0.7" ``` ```rust @@ -180,7 +180,7 @@ All quantities are in SI units. The frame is NED (North-East-Down) throughout. | `R_hub` | `Mat3` | Rotation matrix from hub frame to world (NED) frame. `R_hub * [0,0,1]` gives the hub spin axis in world coordinates. Pass `Mat3::eye()` for a level rotor with +Z spin axis pointing down. | | `v_hub_world` | `Vec3` | Hub-centre velocity in world (NED) frame, m/s. Positive +Z is downward. | | `wind_world` | `Vec3` | Ambient wind velocity in world (NED) frame, m/s. Only the relative velocity `wind_world - v_hub_world` matters -- equal and opposite values produce identical aerodynamics. | -| `omega_rad_s` | `f64` | Rotor angular velocity (rad/s). Positive for the correct spin direction (CCW from above, American convention). The caller advances this externally: `omega += dt * (motor_torque - Q_spin) / I_kgm2`. | +| `omega_rad_s` | `f64` | Rotor angular velocity (rad/s). Positive for the correct spin direction (CCW from above, American convention). The caller advances this externally via the single canonical mechanical ODE stepper: `omega = crate::mechanical::step_omega(omega, Q_spin, motor_torque, I_kgm2, dt, bearing_friction)`. | | `rho_kg_m3` | `f64` | Air density (kg/m^3). Use 1.225 for ISA sea level. | | `t` | `f64` | *(Python API only)* Simulation time (s) carried as a convenience field; not present in the Rust `RotorInputs` struct and not read by any model. | @@ -201,7 +201,7 @@ All quantities are in SI units in the world (NED) frame. |---|---|---| | `F_world` | `Vec3` | Net aerodynamic force on the rotor hub, N, in NED world frame. In hover with the rotor level, thrust appears as a negative Z component (i.e. `F_world.0[2] < 0` lifts the aircraft). | | `M_hub_world` | `Vec3` | Aerodynamic hub moment (pitching/rolling) about the hub centre, N·m, in NED world frame. Corresponds to `M_x`/`M_y` in standard helicopter notation. This is what the airframe sees after blade flap reduction (if `FlapProperties` is set). In NED, positive Mx is roll-right, positive My is pitch-up. | -| `Q_spin` | `f64` | Aerodynamic reaction torque opposing rotation, N·m. Always positive for a powered rotor (torque opposes spin). Use in the angular acceleration equation: `omega += dt * (motor_torque - Q_spin) / I_kgm2`. | +| `Q_spin` | `f64` | Aerodynamic reaction torque opposing rotation, N·m. Always positive for a powered rotor (torque opposes spin). Pass straight through to `crate::mechanical::step_omega` -- no sign flip needed. | | `M_spin` | `Vec3` | `Q_spin` expressed as a vector along the rotor spin axis in world frame, N·m. Direction is `R_hub * [0,0,1]` (the hub +Z axis in world coordinates). Convenience for rigid-body integrators that accept a 3-vector torque input. | ### State Stepping @@ -300,6 +300,9 @@ are equivalent. | apply_flap_reduction() +-- common.rs numerical floors (EPS_*), vrs_lambda1 VRS polynomial +-- cyclic.rs swashplate -> theta_1c, theta_1s mapping + +-- mechanical.rs omega_derivative(), step_omega() -- single canonical + | rigid-body spin ODE (caller-owned), Coulomb bearing + | friction always a parameter (default 0.0) +-- oye.rs OyeBEMModel (annular 2-stage filter) +-- pitt_peters.rs PittPetersModel (3-state L-matrix ODE) +-- polar.rs LinearPolar, TabulatedPolar, Polar trait diff --git a/dynbem_rs/src/aero_io.rs b/dynbem_rs/src/aero_io.rs index da9633d..9dd7384 100644 --- a/dynbem_rs/src/aero_io.rs +++ b/dynbem_rs/src/aero_io.rs @@ -114,13 +114,17 @@ impl std::ops::Mul for Mat3 { // --------------------------------------------------------------------------- /// Inputs supplied by the caller on every compute_forces call. -/// The caller owns omega_rad_s and advances it externally each step: +/// The caller owns omega_rad_s and advances it externally each step via the +/// single canonical mechanical ODE stepper: /// -/// omega += dt * (motor_torque - result.Q_spin) / I_kgm2; +/// omega = crate::mechanical::step_omega( +/// omega, result.Q_spin, motor_torque, I_kgm2, dt, bearing_friction, +/// ); /// /// Use the same dt as the inflow-state integration loop. The returned /// `Q_spin` is the aerodynamic reaction torque on the shaft (positive -/// opposes rotation), so no sign flip is needed in the formula above. +/// opposes rotation), so no sign flip is needed -- pass it straight through +/// to `step_omega`. #[derive(Clone, Debug)] #[allow(non_snake_case)] pub struct RotorInputs { diff --git a/dynbem_rs/src/aero_model.rs b/dynbem_rs/src/aero_model.rs index 796cabc..492c434 100644 --- a/dynbem_rs/src/aero_model.rs +++ b/dynbem_rs/src/aero_model.rs @@ -106,15 +106,20 @@ pub trait AeroModel { // supplied by the caller on every compute_forces call; the mechanical ODE // is NOT part of this state -- the caller owns and integrates omega. // -// Canonical integration pattern (explicit Euler, same dt as inflow loop): +// Canonical integration pattern (same dt as inflow loop): // // let (result, dstate) = aero.compute_forces(&inputs, &state); // state = step_state(state, dstate, dt); // inflow -// omega += dt * (motor_torque - result.Q_spin) / I_kgm2; // spin +// omega = crate::mechanical::step_omega( +// omega, result.Q_spin, motor_torque, i_kgm2, dt, bearing_friction, +// ); // spin // inputs.omega_rad_s = omega; // feed back // -// The Python helper `dynbem.mechanical.omega_derivative` and -// `euler_step_omega` implement the scalar spin ODE. +// `crate::mechanical` (Rust) / `dynbem.mechanical` (Python) is the SINGLE +// canonical place the scalar spin ODE is evaluated and stepped: +// `omega_derivative` for the raw right-hand side, `step_omega` for the +// (semi-implicit, unconditionally stable) integrator. Callers should not +// re-derive `(motor_torque - Q_spin) / I` locally. // The inflow states may be stiff (tau << dt); use the semi-implicit stepper // in `dynbem.mechanical` / `envelope.point_mass._step_state_semi_implicit` // when dt is large relative to the inflow time constants. diff --git a/dynbem_rs/src/lib.rs b/dynbem_rs/src/lib.rs index 2657e9f..894dd3b 100644 --- a/dynbem_rs/src/lib.rs +++ b/dynbem_rs/src/lib.rs @@ -10,6 +10,7 @@ pub mod aero_model; pub mod bem_common; pub mod common; pub mod cyclic; +pub mod mechanical; pub mod oye; pub mod pitt_peters; pub mod polar; diff --git a/dynbem_rs/src/mechanical.rs b/dynbem_rs/src/mechanical.rs new file mode 100644 index 0000000..fe9ea48 --- /dev/null +++ b/dynbem_rs/src/mechanical.rs @@ -0,0 +1,623 @@ +// Rigid-body rotor spin ODE: +// +// I * d(omega)/dt = motor_torque - Q_aero(omega) - Q_friction(omega) +// +// The aero models (QuasiStaticBEM, PittPetersModel, OyeBEMModel, VpmRotor) +// only ever return Q_aero (AeroResult::Q_spin) for the omega they were +// handed -- integrating omega forward in time is caller-owned (see the +// module doc comment on AeroModel::step in aero_model.rs). This module is +// the SINGLE canonical place that ODE is evaluated and stepped: every +// caller in this repo (production code, benches, tests) must call into +// these functions rather than re-deriving `(motor_torque - Q_spin) / I` +// locally -- that duplication is exactly how the sign/physics bugs this +// module fixed crept in before. +// +// Exactly two functions: +// - `omega_derivative` -- the canonical ODE right-hand side. Anything +// that needs d(omega)/dt (diagnostics, custom integrators, tests) +// calls this instead of reimplementing the formula. +// - `step_omega` -- the single canonical integrator (semi-implicit, +// unconditionally stable in the aero term). Use this for all +// production and test time-stepping; there is no separate +// plain-forward-Euler stepper to choose between. +// +// Mirrors the Python API in dynbem.mechanical (omega_derivative, +// step_omega). +// +// Bearing friction is modelled as simple Coulomb (dry) friction: a +// constant-magnitude torque that always opposes the current rotation +// direction. It is a first-class parameter of every function here (not a +// separate optional add-on) -- callers that don't want it simply pass +// `bearing_friction_nm = 0.0`. + +/// Coulomb (dry) bearing friction torque: constant magnitude, opposing the +/// current rotation direction. Zero at exactly `omega == 0` -- this ODE +/// does not model static breakaway torque (a rotor already at rest with no +/// aero/motor torque driving it stays at rest; see the zero-crossing clamp +/// in the stepper functions below for how a decelerating rotor comes to +/// rest in finite time despite friction being explicitly integrated). +#[inline] +fn friction_torque(omega: f64, bearing_friction_nm: f64) -> f64 { + if omega > 0.0 { + bearing_friction_nm + } else if omega < 0.0 { + -bearing_friction_nm + } else { + 0.0 + } +} + +/// Clamp a single integration step so that Coulomb friction (or any other +/// purely-decelerating term integrated explicitly) cannot flip the sign of +/// omega within one step. A constant-magnitude decelerating torque will +/// always eventually overshoot past exactly zero in a finite-step +/// integrator; physically the rotor simply comes to rest sometime during +/// that step and stays there (static friction), it does not reverse +/// direction. This clamp is the standard fix for that. +#[inline] +fn clamp_zero_crossing(old_omega: f64, new_omega: f64) -> f64 { + if old_omega > 0.0 && new_omega < 0.0 { + 0.0 + } else if old_omega < 0.0 && new_omega > 0.0 { + 0.0 + } else { + new_omega + } +} + +/// d(omega)/dt for the rotor rigid-body spin ODE. +/// +/// `omega` is the current rotor speed [rad/s] (needed only to get the sign +/// of the friction torque right). `q_aero` is the aerodynamic reaction +/// torque (`AeroResult::Q_spin`; positive opposes rotation -- use it +/// directly, it is already the shaft reaction). `motor_torque_nm` is the +/// applied shaft torque (positive drives rotation in the direction of +/// omega). `i_ode_kgm2` is the rotor's polar moment of inertia about the +/// spin axis. `bearing_friction_nm` is the Coulomb friction torque +/// magnitude [N.m] (pass `0.0` for a frictionless bearing). +#[inline] +pub fn omega_derivative( + omega: f64, + q_aero: f64, + motor_torque_nm: f64, + i_ode_kgm2: f64, + bearing_friction_nm: f64, +) -> f64 { + let q_friction = friction_torque(omega, bearing_friction_nm); + (motor_torque_nm - q_aero - q_friction) / i_ode_kgm2 +} + +/// Semi-implicit (locally-frozen relaxation) step for the spin ODE. +/// +/// **This is the single canonical integrator for this ODE.** All +/// production code and tests should step omega through this function -- +/// there is no separate plain-forward-Euler stepper to choose between. +/// +/// Freezes the local aerodynamic damping coefficient over the step as +/// `tau = I * omega / Q_aero` (the same "freeze the coefficient/target over +/// dt" assumption used by [`crate::aero_model::IntegrationMethod::SemiImplicitEuler`] +/// for inflow states) and treats it implicitly, while motor torque and +/// Coulomb friction (which doesn't scale with omega, so it can't be +/// linearized the same way) are treated explicitly: +/// +/// ```text +/// explicit_term = omega + dt * (motor_torque - Q_friction) / I +/// omega_new = explicit_term / (1 + dt / tau) +/// ``` +/// +/// This is unconditionally stable in the aero term for any `dt > 0`: when +/// `Q_aero` opposes the current spin direction (the normal no-wind / +/// powered-flight case), the aero contribution alone can never overshoot +/// past zero or oscillate. The explicit friction/motor contribution is +/// subject to a zero-crossing clamp (a rotor decelerating under friction +/// alone comes to rest, it doesn't reverse). Reduces to plain explicit +/// Euler as `dt -> 0` or when there's no meaningful local aero damping to +/// linearize around (`omega == 0`, `Q_aero == 0`, or -- defensively -- +/// `Q_aero` and `omega` have opposite signs, e.g. autorotation driving +/// torque, where there is no local *damping* to treat implicitly). +/// +/// Windmill / autorotation mode (`Q_aero` and `omega` opposite signs, i.e. +/// the aerodynamic torque is adding power to the shaft rather than +/// resisting it) always falls into that last explicit-Euler branch: `tau` +/// would come out negative there, and it isn't safe to treat a *driving* +/// term implicitly with this frozen-coefficient trick (it would divide by +/// a shrinking-toward-zero-or-negative `1 + dt/tau` and blow up). Plain +/// explicit Euler is only conditionally stable (`dt` less than about twice +/// the true local spin-up timescale), same as any other explicit method, +/// but in a properly-coupled simulation `Q_aero` is recomputed from the +/// real aero model at the new `omega` every step -- so as `omega` +/// approaches the equilibrium where the driving and resisting torques +/// balance, the fallback's own inputs shrink and it settles smoothly. See +/// the `test_step_omega_windmill_*` tests below for multi-step +/// confirmation of this at a realistic fixed timestep. +#[inline] +pub fn step_omega( + omega: f64, + q_aero: f64, + motor_torque_nm: f64, + i_ode_kgm2: f64, + dt: f64, + bearing_friction_nm: f64, +) -> f64 { + let q_friction = friction_torque(omega, bearing_friction_nm); + let non_aero_term = omega + dt * (motor_torque_nm - q_friction) / i_ode_kgm2; + + let full_explicit_step = || { + // Full explicit Euler step, including the aero term -- must not be + // silently dropped just because there's no local damping + // coefficient that can be treated implicitly this step. + omega + + dt * omega_derivative( + omega, + q_aero, + motor_torque_nm, + i_ode_kgm2, + bearing_friction_nm, + ) + }; + + let raw = if omega.abs() < 1e-12 { + // At omega == 0 there's no rotation to define a local damping + // timescale from (`tau = I*omega/Q_aero` would be exactly zero, or + // NaN if Q_aero is also zero) -- but Q_aero can still be nonzero + // here (e.g. a windmill/autorotation torque acting on a rotor at + // rest) and must not be dropped from the step. + full_explicit_step() + } else { + let tau = i_ode_kgm2 * omega / q_aero; + if !tau.is_finite() || tau <= 0.0 { + // No positive local damping timescale to linearize around + // (Q_aero == 0, or Q_aero/omega have opposite signs -- e.g. a + // driving/autorotation torque). Fall back to the *full* + // explicit Euler step (including the aero term) -- it must + // not be silently dropped just because it can't be treated + // implicitly this step. + full_explicit_step() + } else { + non_aero_term / (1.0 + dt / tau) + } + }; + + clamp_zero_crossing(omega, raw) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_omega_derivative_matches_formula_no_friction() { + // d_omega = (motor_torque - Q_aero) / I + assert_eq!(omega_derivative(5.0, 10.0, 0.0, 2.0, 0.0), -5.0); + assert_eq!(omega_derivative(5.0, 0.0, 4.0, 2.0, 0.0), 2.0); + assert_eq!(omega_derivative(5.0, 10.0, 10.0, 2.0, 0.0), 0.0); + } + + #[test] + fn test_omega_derivative_friction_opposes_positive_omega() { + // omega > 0: friction subtracts (extra deceleration). + let d = omega_derivative(5.0, 0.0, 0.0, 1.0, 3.0); + assert_eq!(d, -3.0); + } + + #[test] + fn test_omega_derivative_friction_opposes_negative_omega() { + // omega < 0: friction acts in +direction (still decelerating + // toward zero, i.e. opposing the current rotation). + let d = omega_derivative(-5.0, 0.0, 0.0, 1.0, 3.0); + assert_eq!(d, 3.0); + } + + #[test] + fn test_omega_derivative_friction_zero_at_rest() { + // omega == 0: no friction torque (no static breakaway modelled). + let d = omega_derivative(0.0, 0.0, 7.0, 1.0, 3.0); + assert_eq!(d, 7.0); // pure motor torque, no friction subtracted + } + + #[test] + fn test_step_omega_never_overshoots_past_zero() { + // With Q_aero ~ k*omega^2, a huge dt in plain Euler would send + // omega deeply negative. Semi-implicit must stay >= 0. + let omega = 10.0; + let k = 0.5; // Q_aero = k * omega^2 (quadratic drag model) + let q_aero = k * omega * omega; + let i_ode = 1.0; + + for &dt in &[0.01, 0.1, 1.0, 10.0, 1000.0] { + let new_omega = step_omega(omega, q_aero, 0.0, i_ode, dt, 0.0); + assert!( + new_omega >= 0.0, + "semi-implicit step went negative at dt={dt}: {new_omega}" + ); + assert!( + new_omega <= omega, + "semi-implicit step increased omega at dt={dt}: {omega} -> {new_omega}" + ); + } + } + + #[test] + fn test_step_omega_matches_exact_riccati_solution_for_frozen_k() { + // For the frozen-coefficient ODE d(omega)/dt = -k*omega^2 (no motor + // torque, no friction), the exact analytic solution over a step is + // omega_new = omega / (1 + k*omega*dt). step_omega should + // reproduce this exactly (it's the same algebraic form). + let omega = 50.0; + let k = 0.02; + let q_aero = k * omega * omega; + let i_ode = 1.0; + let dt = 5.0; + + let exact = omega / (1.0 + k * omega * dt); + let got = step_omega(omega, q_aero, 0.0, i_ode, dt, 0.0); + assert!( + (got - exact).abs() < 1e-9, + "semi-implicit step {got} != exact frozen-k solution {exact}" + ); + } + + #[test] + fn test_step_omega_reduces_to_explicit_for_small_dt() { + // For tiny dt, step_omega should match plain forward Euler on the + // ODE right-hand side (omega_derivative is the single canonical + // formula for that -- no separate stepper to compare against). + let omega = 20.0; + let q_aero = 3.0; + let i_ode = 1.0; + let dt = 1e-6; + + let explicit = omega + dt * omega_derivative(omega, q_aero, 0.0, i_ode, 0.0); + let semi = step_omega(omega, q_aero, 0.0, i_ode, dt, 0.0); + assert!( + (explicit - semi).abs() < 1e-9, + "step_omega ({semi}) should match explicit Euler ({explicit}) for tiny dt" + ); + } + + #[test] + fn test_step_omega_handles_zero_omega_and_zero_torque() { + // Stopped rotor, no aero torque, no motor torque, no friction: + // stays at rest. + assert_eq!(step_omega(0.0, 0.0, 0.0, 1.0, 0.1, 0.0), 0.0); + } + + #[test] + fn test_step_omega_falls_back_for_driving_torque() { + // Autorotation: Q_aero < 0 while omega > 0 (driving, not damping). + // No positive local damping timescale exists; must fall back to + // the full explicit-Euler derivative rather than blow up or flip + // sign from a negative tau. + let omega = 10.0; + let q_aero = -2.0; + let i_ode = 1.0; + let dt = 0.1; + + let expected = omega + dt * omega_derivative(omega, q_aero, 0.0, i_ode, 0.0); + let got = step_omega(omega, q_aero, 0.0, i_ode, dt, 0.0); + assert_eq!(got, expected); + } + + #[test] + fn test_step_omega_applies_motor_torque_explicitly() { + // Pure motor spin-up from rest, zero aero torque, zero friction: + // reduces to explicit Euler on the forcing term. + let omega = 0.0; + let q_aero = 0.0; + let motor_torque = 5.0; + let i_ode = 2.0; + let dt = 0.5; + + let expected = omega + dt * motor_torque / i_ode; + let got = step_omega(omega, q_aero, motor_torque, i_ode, dt, 0.0); + assert_eq!(got, expected); + } + + #[test] + fn test_step_omega_clamps_zero_crossing_from_friction() { + // Near-zero omega, aero torque tiny, large friction + dt: must not + // flip sign. + let omega = 0.02; + let q_aero = 1e-6; + let bearing_friction_nm = 5.0; + let i_ode = 1.0; + let dt = 1.0; + let new_omega = step_omega(omega, q_aero, 0.0, i_ode, dt, bearing_friction_nm); + assert_eq!(new_omega, 0.0); + } + + #[test] + fn test_friction_allows_restart_from_rest_with_motor_torque() { + // Once clamped to rest, a nonzero motor torque should still be + // able to spin the rotor back up on the next step (friction is + // exactly zero at omega == 0 -- no static breakaway threshold + // modelled). + let omega = 0.0; + let motor_torque = 2.0; + let bearing_friction_nm = 5.0; + let i_ode = 1.0; + let dt = 0.1; + let new_omega = step_omega(omega, 0.0, motor_torque, i_ode, dt, bearing_friction_nm); + assert!( + new_omega > 0.0, + "motor torque should restart rotor from rest" + ); + } + + #[test] + fn test_friction_adds_extra_deceleration_beyond_aero_alone() { + // With friction, decay over a step should be at least as fast as + // with zero friction (monotone comparison). + let omega = 20.0; + let q_aero = 4.0; + let i_ode = 1.0; + let dt = 0.1; + + let no_friction = step_omega(omega, q_aero, 0.0, i_ode, dt, 0.0); + let with_friction = step_omega(omega, q_aero, 0.0, i_ode, dt, 2.0); + assert!( + with_friction < no_friction, + "friction should add extra deceleration: {with_friction} should be < {no_friction}" + ); + } + + #[test] + fn test_omega_derivative_windmill_q_aero_adds_power() { + // Q_aero < 0 while omega > 0 means the aerodynamic torque is + // driving (adding power to) the shaft rather than resisting it -- + // the classic windmill/autorotation regime. domega/dt must be + // positive (the rotor speeds up). + let d = omega_derivative(10.0, -5.0, 0.0, 1.0, 0.0); + assert!( + d > 0.0, + "windmill mode: negative Q_aero should accelerate positive omega, got d={d}" + ); + assert_eq!(d, 5.0); + } + + #[test] + fn test_step_omega_windmill_multistep_converges_without_overshoot() { + // Synthetic "windmill/autorotation" torque curve + // Q_aero(omega) = k*(omega - omega_eq): negative (aero ADDS power, + // spinning the rotor up) below the equilibrium speed, positive + // (normal aerodynamic drag) above it -- the generic shape of a + // real autorotation/windmill torque curve (it crosses zero at the + // equilibrium tip speed) without needing a full BEM aero model. A + // real caller recomputes Q_aero from the actual aero model at the + // *new* omega every step -- this test does the same via the + // closure below, so it exercises step_omega the way production + // code actually drives it, not just a single frozen-Q_aero step. + let k = 2.0; // N.m per rad/s of speed error + let omega_eq = 50.0; // rad/s, equilibrium windmill speed + let i_ode = 1.0; + let dt = 0.0025; // 400 Hz -- the fixed timestep this sim uses + let q_aero_of = |omega: f64| k * (omega - omega_eq); + + let mut omega = 5.0; // start well below equilibrium: aero is driving + assert!( + q_aero_of(omega) < 0.0, + "test setup: must start in the windmill/driving regime" + ); + + let mut prev = omega; + for _ in 0..4000 { + omega = step_omega(omega, q_aero_of(omega), 0.0, i_ode, dt, 0.0); + assert!( + omega >= prev - 1e-9, + "windmill spin-up must be monotonically increasing: {prev} -> {omega}" + ); + assert!( + omega <= omega_eq + 1e-6, + "windmill spin-up must not overshoot equilibrium: {omega} > {omega_eq}" + ); + prev = omega; + } + + assert!( + (omega - omega_eq).abs() < 1e-3, + "windmill spin-up should converge to equilibrium, got {omega}" + ); + } + + #[test] + fn test_step_omega_windmill_spinup_from_rest() { + // omega == 0 exactly hits the `omega.abs() < 1e-12` branch in + // step_omega (no local damping coefficient to freeze at rest), so + // the very first step is explicit -- confirm the rotor still + // starts moving in the windmill regime and keeps converging. + let k = 2.0; + let omega_eq = 50.0; + let i_ode = 1.0; + let dt = 0.0025; + let q_aero_of = |omega: f64| k * (omega - omega_eq); + + let mut omega = 0.0; + let mut prev = omega; + for _ in 0..4000 { + omega = step_omega(omega, q_aero_of(omega), 0.0, i_ode, dt, 0.0); + assert!( + omega >= prev - 1e-9, + "must not decelerate while windmilling from rest: {prev} -> {omega}" + ); + prev = omega; + } + + assert!( + omega > 0.0, + "windmill torque should spin the rotor up from rest" + ); + assert!( + (omega - omega_eq).abs() < 1e-3, + "should converge near equilibrium, got {omega}" + ); + } + + #[test] + fn test_step_omega_windmill_equilibrium_shifts_with_friction() { + // With Coulomb friction opposing rotation, the windmill can only + // spin up to the point where the aero driving torque exactly + // balances friction (Q_aero == -bearing_friction), which is BELOW + // the frictionless zero-crossing equilibrium `omega_eq`. + let k = 2.0; + let omega_eq = 50.0; + let bearing_friction_nm = 4.0; + let i_ode = 1.0; + let dt = 0.0025; + let q_aero_of = |omega: f64| k * (omega - omega_eq); + + let mut omega = 5.0; + let mut prev = omega; + for _ in 0..8000 { + omega = step_omega(omega, q_aero_of(omega), 0.0, i_ode, dt, bearing_friction_nm); + assert!( + omega >= prev - 1e-9, + "windmill spin-up with friction must still be monotonic: {prev} -> {omega}" + ); + prev = omega; + } + + let expected_equilibrium = omega_eq - bearing_friction_nm / k; + assert!( + (omega - expected_equilibrium).abs() < 1e-2, + "friction should shift the windmill equilibrium down to {expected_equilibrium}, got {omega}" + ); + assert!( + omega < omega_eq, + "windmill equilibrium with friction must be below the frictionless equilibrium" + ); + } + + #[test] + fn test_step_omega_continuous_near_q_aero_zero_crossing() { + // Directly probes continuity of step_omega as Q_aero crosses zero + // at fixed omega. This is exactly the boundary between the + // semi-implicit "damping" branch (tau valid) and the explicit + // "driving/fallback" branch (tau invalid) -- a realistic worry + // given this repo has previously had windmill/helicopter branch + // discontinuity bugs in the BEM solver itself (see the + // lambda_climb boundary fix). A tiny change in Q_aero must not + // produce a large jump in the result. + let omega = 35.0; // rad/s, realistic autorotation-range rotor speed + let i_ode = 5.0; + let dt = 0.0025; // 400 Hz + let eps = 1e-4; + + let just_above = step_omega(omega, eps, 0.0, i_ode, dt, 0.0); // damping branch + let at_zero = step_omega(omega, 0.0, 0.0, i_ode, dt, 0.0); // fallback branch + let just_below = step_omega(omega, -eps, 0.0, i_ode, dt, 0.0); // fallback branch + + assert!( + (just_above - at_zero).abs() < 1e-6, + "branch switch at Q_aero=0 should be continuous: {just_above} vs {at_zero}" + ); + assert!( + (at_zero - just_below).abs() < 1e-6, + "fallback branch should be continuous in Q_aero near zero: {at_zero} vs {just_below}" + ); + } + + #[test] + fn test_step_omega_no_chattering_from_sign_flipping_aero_torque_near_equilibrium() { + // Realistic corner case: near the autorotative equilibrium, gusts + // or BEM-solver noise can flip the sign of Q_aero step-to-step + // even though omega itself barely moves. The branch switch + // between semi-implicit and explicit-fallback must not cause + // chattering or a blow-up. + let i_ode = 5.0; + let dt = 0.0025; + let q_sequence = [0.05, -0.05, 0.05, -0.05, 0.02, -0.02, 0.0, 0.05, -0.05]; + + let mut omega = 35.0; + for &q in &q_sequence { + let next = step_omega(omega, q, 0.0, i_ode, dt, 0.0); + assert!( + next.is_finite(), + "chattering aero torque must not produce a non-finite result: {next}" + ); + assert!( + (next - omega).abs() < 0.01, + "chattering aero torque should not cause a large single-step jump: {omega} -> {next}" + ); + omega = next; + } + } + + #[test] + fn test_step_omega_autorotation_entry_crosses_damping_to_driving_branch_smoothly() { + // Realistic engine-failure-in-forward-flight scenario: the rotor + // starts ABOVE the true (zero-torque) autorotative equilibrium + // speed, so it is initially decelerated by ordinary aerodynamic + // drag (Q_aero > 0, the semi-implicit "damping" branch). As it + // slows, Q_aero crosses zero and goes negative (the windmill + // torque now drives the rotor), landing in the explicit-fallback + // "driving" branch, and it settles at the friction-shifted + // equilibrium below the zero-torque point. Unlike the windmill + // tests above (which start already in the driving branch), this + // exercises the actual mid-simulation branch switch. + let i_ode = 5.0; + let c = 20.0; // N.m per rad/s of speed above/below the autorotative eq. + let omega_auto = 35.0; // rad/s, true (zero-torque) autorotation speed + let bearing_friction_nm = 2.0; + let dt = 0.0025; // 400 Hz + let q_aero_of = |omega: f64| c * (omega - omega_auto); + + let mut omega = 45.0; // powered-flight NR, above eq.: engine just failed + assert!( + q_aero_of(omega) > 0.0, + "test setup: must start in the damping regime" + ); + + let mut prev = omega; + let mut crossed_into_driving_branch = false; + for _ in 0..4000 { + let q = q_aero_of(omega); + if q < 0.0 { + crossed_into_driving_branch = true; + } + omega = step_omega(omega, q, 0.0, i_ode, dt, bearing_friction_nm); + assert!( + omega <= prev + 1e-9, + "autorotation entry must decelerate monotonically: {prev} -> {omega}" + ); + assert!( + (prev - omega).abs() < 1.0, + "must not take an unreasonably large single-step jump across the branch boundary: {prev} -> {omega}" + ); + prev = omega; + } + + assert!( + crossed_into_driving_branch, + "test setup: simulation should reach the windmill/driving regime below omega_auto" + ); + + let expected_equilibrium = omega_auto - bearing_friction_nm / c; + assert!( + (omega - expected_equilibrium).abs() < 1e-2, + "should settle at the friction-shifted autorotative equilibrium {expected_equilibrium}, got {omega}" + ); + } + + #[test] + fn test_step_omega_settles_exactly_at_rest_and_stays_there() { + // Realistic rotor-shutdown scenario: motor off, negligible + // residual aero torque, only bearing friction decelerating the + // rotor. Once it reaches exactly rest it must stay there + // indefinitely -- no creeping negative from floating-point + // residue, and no spurious restart from the omega==0 special + // case. + let i_ode = 5.0; + let bearing_friction_nm = 3.0; + let dt = 0.0025; + let mut omega = 2.0; // rad/s, rotor coasting down after shutdown + + for _ in 0..2000 { + omega = step_omega(omega, 0.0, 0.0, i_ode, dt, bearing_friction_nm); + } + + assert_eq!( + omega, 0.0, + "rotor should come to rest and stay there, got {omega}" + ); + } +} diff --git a/dynbem_rs/src/quasi_static_bem.rs b/dynbem_rs/src/quasi_static_bem.rs index 002524e..507ec81 100644 --- a/dynbem_rs/src/quasi_static_bem.rs +++ b/dynbem_rs/src/quasi_static_bem.rs @@ -1257,4 +1257,269 @@ mod tests { rel_err * 100.0 ); } + + /// Omega decay continuity: rotor spinning down to zero with no wind/climb + /// must be smooth (no force/torque jumps), and torque must always resist + /// rotation (Q_spin > 0) since there's no wind to drive it -- second law + /// of thermodynamics: with no energy input from the air, the rotor can + /// only ever lose rotational energy, never gain it. This must hold for + /// every collective setting (positive AND negative pitch), not just one. + /// + /// Q_spin convention (see AeroResult docs / aero_io.rs): positive opposes + /// rotation (drag, motor must supply +torque to sustain); negative means + /// the airflow is driving the rotor (autorotation). With no wind, only + /// Q_spin > 0 is physically possible. + #[test] + fn test_omega_decay_continuity_to_zero() { + use crate::rotor_definition::BladeGeometry; + + // Caradonna-Tung rotor: 2 blades, R=1.143 m, NACA 0012, no twist + let defn = RotorDefinition { + blade: BladeGeometry { + n_blades: 2, + radius_m: 1.143, + root_cutout_m: 0.1, + chord_m: 0.1905, + twist_deg: 0.0, + n_elements: 20, + tip_loss: true, + r_stations_m: Vec::new(), + chord_stations_m: Vec::new(), + twist_stations_deg: Vec::new(), + }, + airfoil: LinearPolarParameters { + CL0: 0.0, + CL_alpha_per_rad: 2.0 * PI, + CD0: 0.008, + alpha_stall_deg: 15.0, + }, + control: None, + pitch_actuation: PitchActuation::DirectMechanical, + flap: None, + name: "Caradonna-Tung".to_string(), + description: String::new(), + }; + + let polar = crate::polar::LinearPolar::from_properties(&defn.airfoil); + let bem = QuasiStaticBEM::build(defn, 36, polar); + + let rho = 1.225; + let base_omega = 1250.0 * PI / 30.0; + + // Sweep across positive and negative collectives -- no wind means + // the rotor must decelerate (Q_spin > 0) regardless of blade pitch + // sign. + for &collective_deg in &[-10.0_f64, -5.0, -2.0, 0.0, 2.0, 5.0, 8.0, 12.0] { + // Geometric decay (constant ratio per step) so that thrust/torque + // coefficients -- which are what should stay continuous, not raw + // magnitudes (those naturally scale like omega^2 by momentum + // theory) -- can be compared step-to-step on an equal footing. + // Stop at a small but nonzero omega; exact omega=0 is covered by + // test_zero_omega_gives_zero_forces. + let mut omega = base_omega; + let mut prev_ct = None; + let mut prev_cq = None; + + for _ in 0..30 { + let inputs = RotorInputs { + collective_rad: collective_deg.to_radians(), + tilt_lon: 0.0, + tilt_lat: 0.0, + R_hub: Mat3::eye(), + v_hub_world: Vec3::zero(), + wind_world: Vec3::zero(), + omega_rad_s: omega, + rho_kg_m3: rho, + }; + + let (res, _) = bem.compute_forces(&inputs, &QuasiStaticRotorState::default()); + + let thrust = res.F_world[2].abs(); + let torque_raw = res.Q_spin; + + // No wind driving the rotor: torque must always resist + // rotation (never accelerate it) for every collective. + assert!( + torque_raw >= 0.0, + "collective={collective_deg} deg, omega={omega:.4}: Q_spin must be \ + >= 0 with no wind (rotor can only decelerate), got {torque_raw:.3e} N.m" + ); + + // Non-dimensionalize by omega^2 (momentum theory: T, Q ~ + // omega^2 at fixed collective/geometry) so we're checking + // for genuine regime-switch discontinuities, not the + // expected quadratic falloff of raw magnitude. + let ct = thrust / (omega * omega); + let cq = torque_raw / (omega * omega); + + if let Some(prev) = prev_ct { + if prev > 1e-6 { + let ratio = ct / prev; + assert!( + 0.5 < ratio && ratio < 2.0, + "collective={collective_deg} deg: thrust coefficient \ + (T/omega^2) jumped at omega={omega:.4}: {prev:.3e} to \ + {ct:.3e} (ratio {ratio:.2})" + ); + } + } + if let Some(prev) = prev_cq { + if prev > 1e-8 { + let ratio = cq / prev; + assert!( + 0.5 < ratio && ratio < 2.0, + "collective={collective_deg} deg: torque coefficient \ + (Q/omega^2) jumped at omega={omega:.4}: {prev:.3e} to \ + {cq:.3e} (ratio {ratio:.2})" + ); + } + } + + prev_ct = Some(ct); + prev_cq = Some(cq); + + omega *= 0.7; // ~30 steps from 1250 RPM down to ~1e-4 RPM + } + } + } + + /// True time-integration spindown: uses the built-in `AeroModel::step()` + /// trait method (aero_model.rs) for the inflow-state update, plus the + /// caller-owned mechanical spin ODE via `crate::mechanical::omega_derivative` + /// / `crate::mechanical::step_omega` -- the single canonical place this + /// ODE is evaluated and stepped (see mechanical.rs module docs). + /// Integrates omega from a nonzero start, with no wind and no motor + /// torque, and checks it decays toward zero monotonically and smoothly -- + /// exactly the caller pattern documented in AeroModel::step's doc + /// comment. + /// + /// Note: with pure aerodynamic (quadratic) drag and zero other losses, + /// the continuous-time ODE only reaches exactly omega=0 asymptotically + /// (1/t tail) -- so this checks decay to a small fraction of the start + /// speed, not literal zero (test_zero_omega_gives_zero_forces already + /// covers the omega==0 endpoint directly). + #[test] + fn test_omega_spindown_time_integration() { + use crate::aero_model::IntegrationMethod; + use crate::rotor_definition::BladeGeometry; + + let defn = RotorDefinition { + blade: BladeGeometry { + n_blades: 2, + radius_m: 1.143, + root_cutout_m: 0.1, + chord_m: 0.1905, + twist_deg: 0.0, + n_elements: 20, + tip_loss: true, + r_stations_m: Vec::new(), + chord_stations_m: Vec::new(), + twist_stations_deg: Vec::new(), + }, + airfoil: LinearPolarParameters { + CL0: 0.0, + CL_alpha_per_rad: 2.0 * PI, + CD0: 0.008, + alpha_stall_deg: 15.0, + }, + control: None, + pitch_actuation: PitchActuation::DirectMechanical, + flap: None, + name: "Caradonna-Tung".to_string(), + description: String::new(), + }; + + let polar = crate::polar::LinearPolar::from_properties(&defn.airfoil); + let bem = QuasiStaticBEM::build(defn, 36, polar); + + let rho = 1.225; + let i_ode_kgm2 = 1.0_f64; + let omega_start = 1250.0 * PI / 30.0; + let max_steps = 2000; + + for &collective_deg in &[-10.0_f64, -5.0, 5.0, 8.0] { + let mut omega = omega_start; + let mut state = bem.initial_state(); + let mut dt = 1e-3_f64; + let mut reached_target = false; + + for step in 0..max_steps { + let inputs = RotorInputs { + collective_rad: collective_deg.to_radians(), + tilt_lon: 0.0, + tilt_lat: 0.0, + R_hub: Mat3::eye(), + v_hub_world: Vec3::zero(), + wind_world: Vec3::zero(), + omega_rad_s: omega, + rho_kg_m3: rho, + }; + + // Built-in trait method: advances the inflow state (a + // no-op here -- QuasiStaticBEM carries zero inflow DOF). + let (res, new_state) = + bem.step(&inputs, &state, dt, IntegrationMethod::ExplicitEuler); + state = new_state; + + // Mechanical spin ODE is caller-owned by design (see + // AeroModel::step doc comment) -- call into the single + // canonical mechanical ODE functions rather than + // re-deriving the formula locally. + let motor_torque_nm = 0.0; + let d_omega = crate::mechanical::omega_derivative( + omega, + res.Q_spin, + motor_torque_nm, + i_ode_kgm2, + 0.0, + ); + + // No wind/motor driving the rotor: it must never + // accelerate (Q_spin >= 0 always, per the physics + // established in test_omega_decay_continuity_to_zero). + assert!( + d_omega <= 1e-9, + "collective={collective_deg} deg, step {step}: omega must \ + not accelerate with no wind/motor torque, d_omega={d_omega:.3e}" + ); + + // Adapt dt to keep the per-step relative change bounded + // (~5%) -- necessary because the aerodynamic torque + // shrinks as omega^2, so a fixed dt would either take + // forever near omega_start or blow past zero near the + // end. + if d_omega.abs() > 1e-12 { + dt = (0.05 * omega / d_omega.abs()).clamp(1e-6, 5.0); + } + + let new_omega = crate::mechanical::step_omega( + omega, + res.Q_spin, + motor_torque_nm, + i_ode_kgm2, + dt, + 0.0, + ); + assert!( + new_omega <= omega + 1e-9, + "collective={collective_deg} deg, step {step}: omega \ + increased from {omega:.6} to {new_omega:.6} during spindown" + ); + + omega = new_omega; + + if omega <= 1e-3 * omega_start { + reached_target = true; + break; + } + } + + assert!( + reached_target, + "collective={collective_deg} deg: omega failed to decay to \ + 0.1% of its starting value within {max_steps} steps (stuck \ + at {omega:.4} rad/s -- possible spindown regression)" + ); + } + } } diff --git a/envelope/point_mass.py b/envelope/point_mass.py index f0967c0..d70098f 100644 --- a/envelope/point_mass.py +++ b/envelope/point_mass.py @@ -38,7 +38,7 @@ import numpy as np from dynbem import RotorInputs, create_aero -from dynbem.mechanical import omega_derivative +from dynbem.mechanical import step_omega from dynbem.rotor_definition import RotorDefinition G = 9.81 # m/s² @@ -222,8 +222,8 @@ def simulate_point( # Clip dynamic-inflow states (guards against divergence at very # short time constants or model start-up transients). state = state.from_array(_clip_state(arr)) - omega = max(0.5, min(300.0, - omega + dt * omega_derivative(aero.Q_spin, 0.0, I_ode))) + omega_new, _ = step_omega(omega, 0.0, aero.Q_spin, 0.0, I_ode, dt) + omega = max(0.5, min(300.0, omega_new)) # 1-D force integration along tether (explicit Euler). The previous # semi-implicit form used a finite-difference probe of ∂F/∂v_along @@ -395,7 +395,7 @@ def _step(tension_now: float, sim_t: float, dt_in: float) -> None: if not np.array_equal(arr_raw, arr_clipped): aero_clamped = True state = state.from_array(arr_clipped) - omega_new = omega + dt_in * omega_derivative(aero.Q_spin, 0.0, I_ode) + omega_new, _ = step_omega(omega, 0.0, aero.Q_spin, 0.0, I_ode, dt_in) omega_clipped = max(0.5, min(300.0, omega_new)) if omega_clipped != omega_new: aero_clamped = True diff --git a/tests/test_oye.py b/tests/test_oye.py index 7992d78..81db206 100644 --- a/tests/test_oye.py +++ b/tests/test_oye.py @@ -287,7 +287,7 @@ def test_descent_with_sidewind_converges(self, defn): tether_hat, _build_r_hub, G, ) - from dynbem.mechanical import omega_derivative + from dynbem.mechanical import step_omega m = create_aero(defn, model="oye") # Use the beaupoil rotor -- the one that exhibited the regression. @@ -324,8 +324,8 @@ def test_descent_with_sidewind_converges(self, defn): arr = _step_state_semi_implicit(m, state, dst, dt, inp) arr = _clip_state(arr) state = state.from_array(arr) - omega = max(0.5, min(300.0, - omega + dt * omega_derivative(aero_r.Q_spin, 0.0, I_ode))) + omega_new, _ = step_omega(omega, 0.0, aero_r.Q_spin, 0.0, I_ode, dt) + omega = max(0.5, min(300.0, omega_new)) f_thrust = float(np.dot(aero_r.F_world, t_hat)) f_along = f_thrust + mass * G * float(t_hat[2]) + T_tether v_along = max(-30.0, min(30.0, v_along + dt / mass * f_along)) diff --git a/vpm_viz/src/main.rs b/vpm_viz/src/main.rs index 30ad79c..d837d43 100644 --- a/vpm_viz/src/main.rs +++ b/vpm_viz/src/main.rs @@ -376,9 +376,18 @@ impl eframe::App for VpmVizApp { self.last_result = result; self.step_count += 1; - // Integrate omega: constant drive torque minus aerodynamic resistive torque. - let d_omega = (Q_DRIVE - self.last_result.torque) / I_ROTOR * DT; - self.omega = (self.omega + d_omega).clamp(OMEGA_MIN, OMEGA_MAX); + // Integrate omega: constant drive torque minus aerodynamic resistive + // torque, via the single canonical mechanical ODE stepper (no + // bearing friction modelled in this viz). + self.omega = dynbem_rs::mechanical::step_omega( + self.omega, + self.last_result.torque, + Q_DRIVE, + I_ROTOR, + DT, + 0.0, + ) + .clamp(OMEGA_MIN, OMEGA_MAX); // Ramp wind from 0 to V_WIND_MAX over WIND_RAMP_STEPS steps. self.v_wind = (V_WIND_MAX * (self.step_count as f64 / WIND_RAMP_STEPS as f64))