Metadata-Version: 2.4
Name: dynbem
Version: 0.1.0
Summary: Blade-element momentum rotor aerodynamics with dynamic inflow (Pitt-Peters and Øye), covering the full operating envelope from helicopter hover through climb, descent, VRS, autorotation, and wind-turbine extraction in a single code path.
Author: Kristof
License-Expression: MIT
Project-URL: Homepage, https://github.com/mcroomp/dynbem
Project-URL: Repository, https://github.com/mcroomp/dynbem
Project-URL: Issues, https://github.com/mcroomp/dynbem/issues
Keywords: rotor,aerodynamics,BEM,blade-element-momentum,dynamic-inflow,Pitt-Peters,Oye,helicopter,wind-turbine,autorotation,VRS
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.11
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: numba>=0.58
Provides-Extra: yaml
Requires-Dist: PyYAML>=6; extra == "yaml"
Provides-Extra: dev
Requires-Dist: pytest; extra == "dev"
Requires-Dist: pytest-timeout; extra == "dev"
Requires-Dist: PyYAML>=6; extra == "dev"
Requires-Dist: build; extra == "dev"
Requires-Dist: twine; extra == "dev"
Dynamic: license-file

# dynbem

**Dynamic blade-element momentum rotor aerodynamics — helicopter and
wind-turbine modes in one code path.**

`dynbem` is a Python rotor-aerodynamics library built around a multi-element
blade-element-momentum (BEM) solver coupled to dynamic-inflow models. It is
designed to be numerically valid across the **full operating envelope** —
helicopter hover, axial climb, axial descent, vortex-ring state (VRS),
windmill-brake state (WBS), autorotation, and wind-turbine power extraction
— without switching equations or sign conventions between regimes.

Two dynamic-inflow models are provided:

- **Pitt-Peters** (three-state global ν₀/ν_s/ν_c) — Numba-JIT-compiled,
  with the Peters L-matrix, Glauert wake-skew via the off-diagonal
  coupling, and the Leishman empirical VRS polynomial baked into the
  uniform-inflow state.
- **Øye 2-stage annular** — per-annulus filtered momentum inflow (the
  OpenFAST DBEMT formulation), independent across radii and numerically
  stable at high advance ratios where Pitt-Peters becomes stiff.

Both models share a tabulated polar interpolator and a common BEM ψ-loop
kernel, and they plug into the same `AeroBase` interface. The package also
includes a flight-envelope sweep driver (`envelope/compute_map.py`), a
cyclic-trim solver, a point-mass + cyclic-pitch attitude simulator, and a
test suite that validates against NACA TN-2474 (Castles & Gray) and Peters'
own Nikolsky-lecture data.

Coordinates are NED throughout; rotor rotation is CCW-from-above
(American helicopter convention).

## Install

The package is set up with a standard `pyproject.toml` and can be installed
into any environment:

```
pip install -e .
```

For the bundled `.venv` on Windows, run `setup.cmd` from the repo root —
it creates `.venv\`, upgrades pip, and installs `requirements.txt`.
Activate with `.venv\Scripts\activate`, or invoke directly via
`.venv\Scripts\python` / `.venv\Scripts\pytest`.

## Usage

```python
import numpy as np
import dynbem

defn   = dynbem.rotor_definition.load("rotors/castles_gray_6ft/rotor.yaml")
model  = dynbem.create_aero(defn, model="pitt_peters_jit")  # or "pitt_peters", "oye", "bem"
state  = model.initial_rotor_state()

inputs = dynbem.RotorInputs(
    collective_rad=0.14,
    tilt_lon=0.0, tilt_lat=0.0,                # swashplate (helicopter-standard signs)
    R_hub=np.eye(3),
    v_hub_world=np.zeros(3),
    wind_world=np.zeros(3),
    t=0.0,
)
result, derivative = model.compute_forces(inputs, state)
# result.F_world, result.M_orbital, result.M_spin, result.Q_spin
# derivative is a RotorState with d/dt of every state field
```

`dynbem/__init__.py` lists the full public surface — models (`BEMModel`,
`PittPetersModel`, `PittPetersModelJIT`, the `create_aero` factory),
inputs/outputs (`RotorInputs`, `AeroResult`), state types
(`QuasiStaticRotorState`, `PittPetersRotorState`), polars (`AirfoilPolar`,
`LinearPolar`), rotor-definition types (`RotorDefinition`,
`BladeGeometry`, `AirfoilProperties`, `ControlProperties`, etc.), and the
`vrs_lambda1` helper. Cyclic mapping (`tilt_lon`/`tilt_lat` →
blade-pitch coefficients) lives in [dynbem/cyclic.py](dynbem/cyclic.py).

## Flight envelope sweep

```
run_map.cmd                  # quick grid, saves to out\map.npz, plots to out\
run_map.cmd --full           # full grid
python -m envelope.compute_map --help
```

## Tests

```
.venv\Scripts\pytest
```

The `tests/` directory contains unit tests, validation scripts against
published rotor data, and end-to-end force-balance / frame-transform
checks. Key files:

- `test_bem.py`, `test_bem_components.py` — Level-1 BEM unit tests
  (Prandtl tip/hub loss, momentum-BEM convergence, interface).
- `test_castles_gray.py` — full Castles-Gray (NACA TN-2474) hover CT/CQ,
  WBS, and autorotation sweep against the paper's Table I data.
- `test_pitt_peters.py`, `test_pitt_peters_jit.py` — Pitt-Peters
  hover/VRS/WBS validation; JIT-vs-reference numerical agreement.
- `test_cyclic.py` — helicopter-standard cyclic sign convention on all
  three models, plus Pitt-Peters inflow-dynamics response to cyclic.
- `test_force_balance.py` — hub-frame transform (`R_hub` tilt → force
  vector tilt), hover collective sweep against vehicle weight,
  translating-flight + wind, and end-to-end swashplate phase mapping.
- `val_step*.py`, `validate_table_i.py` — standalone validation
  scripts (not part of the pytest run; invoked manually).

---

## Coordinate system — NED

This project uses **NED (North-East-Down)** throughout, without exception:

- X = North, Y = East, Z = Down
- Gravity acts in the **+Z** direction
- Rotor thrust (upward lift) is **negative Z** in world frame: `F_world[2] < 0`
- Wind blowing upward (driving a flying turbine) is **negative Z** in world frame
- `R_hub` rotates from hub frame → NED world frame

### Reading literature — coordinate trap

Most helicopter and wind-turbine literature uses one of:
- **SAE / helicopter**: X forward, Y right, Z down (body frame, not world NED)
- **Wind-turbine (IEC 61400)**: X downwind, Y lateral, Z up (**ENU-like**)
- **Aeronautics (NED)**: X North, Y East, Z Down

When adapting equations or sign conventions from papers, always check
which frame the authors use. Windmill-brake-state and axial-induction
literature (Glauert, Buhl) often defines positive inflow **upward**
(opposing thrust), which is **negative Z** here. Flip signs accordingly.

### Inflow sign convention (NED)

For a rotor disk lying in the XY-plane (hub pointing down):

- `lambda` (inflow ratio) is positive when flow passes through the disk
  from above (downward, +Z direction), i.e. in **normal rotor mode**
  (helicopter hover).
- In **windmill / autorotation mode** the wind drives flow upward (−Z),
  so `lambda` is **negative** when the rotor is in energy-harvesting mode.
- Collective pitch `theta_0 > 0` pitches blade leading edge up
  (toward −Z thrust).

---

## Implementation roadmap

The model is built in phases from simple to state-of-the-art, each a
drop-in upgrade behind the same `AeroBase` interface.

### Level 1 — Multi-element quasi-static BEM ✅ DONE

- Multi-element BEM loop (radial quadrature over `n_elements` annuli)
- Hover-safe inflow iteration on `λ_r` (not wind-turbine induction
  factor `a`)
- Per-element Prandtl **tip + hub** loss `F = F_tip · F_hub` (both
  factors exported from `dynbem.bem`)
- Glauert / Buhl Windmill Brake State correction (quadratic root
  selection)
- Forward-flight ψ-loop: per-azimuth blade pitch (cyclic), tangential
  wind projection (advancing/retreating), and in-plane hub moment
  accumulation (M_orbital)
- `dω/dt = (Q_aero + Q_motor) / I_ode` — rotor speed integrated as ODE
  state
- `dψ/dt = ω` — spin angle integrated as ODE state
- Returns `QuasiStaticRotorState` derivative
- **Validation target**: Caradonna-Tung rotor (NASA TM-81232, 1981) —
  2-blade NACA 0012, CT vs collective at 5°/8°/12°
  (see `Research/CaradonnaTung/` for CT tables and test notes)

### Level 2 — Pitt-Peters 3-state dynamic inflow ✅ DONE

- `PittPetersModel` (numpy) and `PittPetersModelJIT` (Numba-compiled, same
  physics) in `dynbem/pitt_peters.py` / `dynbem/pitt_peters_jit.py`
- Prescribed-inflow blade element loop with per-element Prandtl tip + hub
  loss; blade sees `λ_total = λ_0 + v_climb/ΩR` (induced state +
  freestream), so WBS and autorotation work correctly
- Pitt-Peters ODE in matrix form (Peters 2009 Eq 7, hub axes):
  `[M] dλ/dt + V·[L]⁻¹ λ = forcing` with
  `M = diag(8/(3π), 16/(45π), 16/(45π))` → `τ_0 = 8R/(3πV_T)`,
  `τ_cs = 16R/(45πV_T)`.
- Steady-state targets follow the canonical L matrix (Peters Eq 10) with
  X = tan(χ/2), translated to our ψ=0-at-+X convention:
  ```
  λ_0_ss = C_T/(2·µ_T)        + (15π·X/64) · C_M_hub / µ_T
  λ_c_ss = −(15π·X/64) · C_T  + 4·cos(χ)/(1+cos χ) · C_M_hub) / µ_T
  λ_s_ss =                      4/(1+cos χ)         · C_L_hub  / µ_T
  ```
  The `−(15π·X/64)·C_T` cross-coupling in λ_c_ss is the Pitt-Peters term
  that produces Glauert wake-skew naturally from thrust forcing — no
  closed-form Glauert tilt needed.
- Cyclic input (`tilt_lon`, `tilt_lat`) wired through both models:
  blade pitch `θ(ψ) = collective + θ_1c·cos(ψ) + θ_1s·sin(ψ)` with
  helicopter-standard signs (`tilt_lon > 0` → nose-down,
  `tilt_lat > 0` → roll right). See `dynbem/cyclic.py` and CLAUDE.md.
- In-plane hub moments returned via `AeroResult.M_orbital`
  (`Mx_hub, My_hub` accumulated in the ψ-loop) — needed for cyclic to
  produce vehicle attitude response in the outer loop.
- **VRS empirical correction** (Leishman 2000, fit to Castles-Gray
  data): in 0 < λ₂ < 2,
  `λ_0_ss` comes from the polynomial
  `λ₁/V_h = 1 + 1.125λ₂ − 1.372λ₂² + 1.718λ₂³ − 0.655λ₂⁴`
  rather than momentum theory, preventing the Level-1 CT blow-up in VRS.
  Cross-coupling is also skipped in the VRS regime.
- **Canonical reference**: Peters, D.A. (2009), "How Dynamic Inflow
  Survives in the Competitive World of Rotorcraft Aerodynamics: The
  Alexander Nikolsky Honorary Lecture," *JAHS* 54(1):011001. PDF and
  extraction notes in `Research/Peters_Nikolsky_2008/`.
- **Validation**: `tests/test_pitt_peters.py` — hover CT vs Level-1
  BEM, VRS no-blow-up, WBS autorotation sign, first-order inflow lag.
  `tests/test_cyclic.py` — helicopter-standard cyclic signs on all
  three models, inflow dynamics under hover + cyclic.
- **Known limitations**:
  - VRS CT still rises to ~2× nominal in deep VRS (λ₂ ≈ 1.5–2) at
    fixed θ; real rotor stays near nominal (paper: θ barely adjusts).
    The Leishman polynomial shifts the operating point but doesn't
    fully suppress it.
  - Autorotation torque crossing at V/ΩR ≈ 0.14 vs paper's 0.083.
  - Mass-flow scaling uses `µ_T = √(µ²+λ²)` (classical Glauert) rather
    than Peters' Eq 8 `V = (µ²+(λ+ν)(λ+2ν))/√(µ²+(λ+ν)²)`. They agree
    in high-speed forward flight but differ by 2× in hover. Switching
    would need validation against a hover dataset.
  - Wind-axis rotation of the L-matrix is NOT applied; oblique flight
    `µ_y ≠ 0` is approximate.  Exact for axial and pure-longitudinal
    flight.  A previous implementation was reverted because it
    destabilised the tethered-rotor envelope — see CLAUDE.md.

### Level 2 alt — Øye 2-stage annular dynamic inflow ✅ DONE

- `OyeBEMModel` in `dynbem/oye.py` (Numba-compiled ψ-loop)
- **Annulus-local** inflow: each radial annulus has its own pair of
  first-order lag filters `(W_int, W)` chasing the quasi-steady
  momentum target `W_qs`.  No global L-matrix; no λ_c/λ_s harmonic
  states.
- Two time constants per annulus (Øye 1990, OpenFAST AD Theory §6.3.4):
  ```
  τ₁ = 1.1 / (1 − 1.3·min(a, 0.5)) · R / V_∞
  τ₂(r) = (0.39 − 0.26·(r/R)²) · τ₁
  τ₁·dW_int/dt + W_int = W_qs + k·τ₁·dW_qs/dt
  τ₂·dW/dt + W = W_int
  ```
  with empirical coupling `k = 0.6`.  DBEMT_Mod=1 equivalent
  (`dW_qs/dt = 0` across each outer step — exact for envelope sweeps).
- W_qs per annulus from Glauert momentum balance using rotor-mean
  `µ_T = V_T / Ω·R`:  `W_qs[i] = dCT/dx[i] / (4·x[i]·µ_T)`.
- Same VRS override (Leishman polynomial) as Pitt-Peters for
  `0 < V_descent/V_h < 2` — applied uniformly across annuli.
- Same cyclic-pitch wiring (`tilt_lon` / `tilt_lat` → per-ψ blade
  pitch) and same in-plane hub moments returned via `M_orbital`.
- **Why this alongside Pitt-Peters**: Pitt-Peters' L-matrix couples
  thrust + hub moments back into all three inflow harmonics
  globally, which produces a stiff BEM-driven feedback at high
  advance ratios and in descent + edgewise wind.  Øye's annulus-local
  filters are independent → no feedback loop → numerically stable in
  the same regimes that needed adaptive time-stepping with
  Pitt-Peters.  OpenFAST's DBEMT uses the same Øye-style formulation
  for this reason.
- **Trade-off**: no harmonic inflow states means the inflow doesn't
  develop a `λ_c`-like tilt in response to cyclic pitching moments,
  so `tests/test_cyclic.py::test_cyclic_inflow_reduces_hub_moment`
  (which checks PP's specific feedback mechanism) doesn't apply.
  Cyclic *control* still works (hub moments respond correctly to
  swashplate inputs), but cyclic *inflow feedback* is absent.
- **Validation**: `tests/test_oye.py` — hover CT vs Pitt-Peters (5%
  match at moderate thrust), climb-vs-hover induction sign, finite-τ
  lag response, descent + edgewise wind convergence where Pitt-Peters
  was numerically stiff.
- **References**:
  - Øye, S. (1990).  A simple vortex model.  IEA Symposium.
  - Snel, H. & Schepers, J.G. (1995).  Joint investigation of dynamic
    inflow effects.  ECN.
  - OpenFAST AeroDyn Theory v3.5, §6.3.4 (DBEMT).

### Level 3 — Peters-He finite-state dynamic inflow (state of the art)

- 9-state (or higher-order) Peters-He inflow model
- Requires new `PetersHeRotorState` dataclass
- Captures higher harmonics of the inflow distribution
- Best accuracy for maneuvering flight and aeroelastic coupling
- **Validation target**: Caradonna-Tung unsteady / forward-flight data

### Forward flight (applies to all levels) — implemented

- Oblique inflow: advance ratio `µ = V_edge / (Ω·R)` ≠ 0
- Blade azimuth-dependent velocity in the BEM loop (`n_psi=36` stations
  by default, triggered when `µ > 0.01`, cyclic input is nonzero, or
  cyclic inflow state is nonzero)
- In-plane hub moments `Mx_hub`, `My_hub` returned in `AeroResult.M_orbital`
- Pitt-Peters L matrix off-diagonal `−L_off·C_T` produces Glauert
  wake-skew naturally from thrust forcing (exact for axial and
  pure-longitudinal flight; approximate for oblique `µ_y ≠ 0`)

---

## BEM solver design — critical notes

### Hover-safe inflow iteration

The standard wind-turbine BEM uses the induction factor `a = v_i / V_inf`,
which **collapses to zero in hover** (`V_inf = 0`). This code instead
iterates on the **total inflow ratio** `λ_r = v_a / (Ω·R)`, where `v_a`
is the total axial velocity at the disk (external freestream + induced).

The combined momentum-BEM equation at each annulus is:

    k·(λ_r² + x²) = λ_r·(λ_r − λ_c)

where `k = σ_r·cn / (8·F)`, `x = r/R`, and `λ_c = v_climb / (Ω·R)`.

This quadratic is solved per iteration step; `v_climb = 0` in hover is
handled naturally (gives the standard hover solution
`λ_r = x·sqrt(k/(1−k))`).

### v_climb sign convention (internal BEM)

`v_climb = dot(v_rel_world, hub_axis_ned)` (no negation):

- `v_climb > 0`: air flows **downward** through disk (helicopter climb /
  normal inflow)
- `v_climb = 0`: hover
- `v_climb < 0`: air flows **upward** through disk (autorotation /
  flying wind turbine)

### Root selection in the momentum-BEM quadratic

The quadratic has two roots. Selection is by operating mode:

- Helicopter / hover (`λ_c ≥ 0`): take the **positive** root (`λ_r > 0`)
- Turbine / autorotation (`λ_c < 0`): take the **negative** root
  (`λ_r < 0`)

### Autorotation torque sign

In autorotation (upward wind, `λ_c < 0`):
- `λ_r < 0` → `φ < 0` → `ct = cl·sin(φ) − cd·cos(φ) < 0` → `Q_total < 0`
- `d_omega = (−Q_total + Q_motor) / I` → positive angular acceleration ✓

In powered/hover mode (`λ_c ≥ 0`):
- `Q_total > 0` (aerodynamic drag on rotor) → `d_omega < 0` without
  motor torque ✓

### Force direction

`F_world = −T_total · hub_axis_ned`

`T_total` is always positive for a rotor generating lift (cn > 0 in
both modes). With `hub_axis_ned = [0, 0, 1]` for a level rotor:
`F_world[2] = −T_total < 0` (upward). ✓

---

## Pitt-Peters design notes (`dynbem/pitt_peters.py`)

### State interpretation

`λ_0` (and `λ_c`, `λ_s`) is the **induced** inflow ratio `v_i / (ΩR)`,
not the total inflow. The total axial flow seen by each blade element
is:

    λ_total = λ_0 + λ_climb    where  λ_climb = v_climb / (ΩR) < 0 in descent

This must be computed inside the blade element loop — **do not pass
only `λ_0`**. Without the freestream term the blade never sees
net-upward flow in WBS, so CQ never goes negative and autorotation is
suppressed entirely.

### VRS polynomial sign convention

The Leishman (2000) polynomial uses descent-positive
λ₂ = V_descent / V_h:

    λ₁/V_h = 1 + 1.125·λ₂ − 1.372·λ₂² + 1.718·λ₂³ − 0.655·λ₂⁴

This is NOT the form with coefficients
(−1.125, −1.372, −1.718, −0.655), which applies when the argument is
V_climb/V_h (negative for descent). The two forms are equivalent;
this code uses descent-positive throughout.

### V_T floor

`V_T = |v_climb + v_0|` → 0 in the middle of VRS (upward freestream ≈
downward induced). A floor of `1e-2 · max(ΩR, 1)` prevents
`τ_0 → ∞` and division by zero. This is physically reasonable:
`τ_0 → large` in VRS is correct (slow, unsteady response), and the
exact floor value doesn't matter for stability.

### Why CT still rises in deep VRS

At λ₂ ≈ 1.5, the Leishman polynomial gives
`λ_0_ss ≈ 2 · V_h/ΩR`. Combined with `λ_climb ≈ −1.5 · V_h/ΩR`,
the net blade inflow `λ_total ≈ 0.5 · V_h/ΩR` is less than hover, so
AoA increases and CT rises. The real VRS has recirculating wakes that
further restrict net throughflow; the 1-D polynomial captures the mean
induced velocity but not the 3-D blockage. This is a known limitation
of all momentum-based VRS models.

---

## Øye design notes (`dynbem/oye.py`)

### State interpretation

`W[i]` and `W_int[i]` are induced inflow ratios `v_i / (Ω·R)` **per
annulus**, not global harmonics. The total axial flow at annulus `i`
seen by the blade is `λ_total[i] = λ_climb + W[i]` (compare with
Pitt-Peters' `λ_total = λ_climb + λ_0 + x·(λ_c·cos ψ + λ_s·sin ψ)`).

`W` is what the blade actually reads in the ψ-loop. `W_int` is the
intermediate filter stage between the quasi-steady target `W_qs[i]`
and `W`. Both arrays have length `n_elements`.

### Quasi-steady target

`W_qs[i]` is solved per annulus from Glauert momentum balance using
the rotor-mean `µ_T = V_T / Ω·R`:

    W_qs[i] = dCT/dx[i] / (4·x[i]·µ_T)
    where  V_T = √(v_edge² + (v_climb + v_0_mean)²)

This linear (in `W_qs`) form is what Pitt-Peters effectively uses in
its aggregate `λ_0_ss = T / (2ρA·V_T·ΩR)`. The pure axial-momentum
form `4·x·λ_r·W = dCT/dx` is unstable in forward flight (small λ_r in
descent makes W blow up) and was rejected during development.

### Why no L matrix

Annulus-local: each `W[i]` evolves independently, driven only by
`W_qs[i]` from its own annulus. Cross-annulus coupling happens only
through the rotor-mean `µ_T` in the τ formulas and `V_h` in the VRS
override. There's no analogue of Pitt-Peters' `−L_off·C_T` term that
feeds total thrust into the cyclic harmonics, so no BEM-driven
feedback loop and no associated stiffness — at the cost of not
modelling cyclic inflow harmonics at all.

### Time constants

`τ₁` is rotor-mean (depends on `a_avg`, not per-annulus); `τ₂(r)`
varies with radius. With `dt = 5 ms` and a 1 m rotor at `V_∞ ~ 10 m/s`,
`τ₁ ~ 0.1 s` and `τ₂ ~ 0.04 s` — both well above the envelope's outer
`dt`, so the semi-implicit Euler in `envelope/point_mass.py` is gentle
damping at most.

### Cyclic input

Cyclic pitch flows through the same `cyclic_coeffs` → `θ(ψ) =
collective + θ_1c·cos ψ + θ_1s·sin ψ` path as Pitt-Peters; the ψ-loop
produces correct hub moments. What's *missing* compared to
Pitt-Peters: the cyclic-driven hub moment doesn't develop a
counter-acting inflow harmonic (no `λ_c`/`λ_s` states), so the
steady-state moment is over-predicted vs Pitt-Peters at hover.
Cyclic *control* (sign and order-of-magnitude) is right; cyclic
*inflow damping* is absent.

---

## Research sources

Extracted tables and figures from primary literature live under
`Research/`. Each paper subfolder uses the convention
`page_NN_<description>.md` so extractions trace back to their source
page image.

- **CaradonnaTung/** — NASA TM-81232 (1981). 2-blade NACA 0012 hover
  CT data at θc = 5°/8°/12°. Primary BEM validation source. No CP /
  torque data.
- **Buhl_NREL_TP500_36834/** — NREL TP-500-36834 (2005). Windmill
  Brake State correction extending Glauert. Used for the WBS
  quadratic.
- **Castles_TN2474/** — NACA TN-2474 (Castles & Gray, 1951). Induced
  velocity in hover/descent — experimental basis for the Leishman VRS
  polynomial.
- **Harrington_TN2318/** — NACA TN-2318 (Harrington, 1951). Hover
  CT vs CP polars for two full-scale rotors. Candidate dataset for
  CP-CT polar validation.

## License

MIT — see `LICENSE`.
