feat(mechanical): single canonical mechanical spin ODE + bearing friction

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"]

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.

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),

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)`

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"

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

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

dynbem/src/lib.rs

@@ -30,9 +30,55 @@ fn cyclic_coeffs(tilt_lon: f64, tilt_lat: f64, control: Option<PyControlProperti
     dynbem_rs::cyclic::cyclic_coeffs(tilt_lon, tilt_lat, gains)
 }
 
+/// Rigid-body rotor spin ODE derivative. See dynbem_rs::mechanical for the
+/// full physics (Coulomb bearing friction is always a parameter; pass
+/// bearing_friction_nm=0.0 for a frictionless bearing).
+#[pyfunction]
+#[pyo3(signature = (omega, q_aero, motor_torque_nm, i_ode_kgm2, bearing_friction_nm = 0.0))]
+fn omega_derivative(
+    omega: f64,
+    q_aero: f64,
+    motor_torque_nm: f64,
+    i_ode_kgm2: f64,
+    bearing_friction_nm: f64,
+) -> 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::<PyLinearPolar>()?;
     m.add_class::<PyTabulatedPolar>()?;
     m.add_class::<PyBladeGeometry>()?;

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