Theory

What pybmodes models, what assumptions it carries, what references it cross-validates against, and the non-obvious implementation decisions that shape the public API.

The structural model

pybmodes solves the linearised modal eigenproblem for a slender beam with coupled flap, lag, torsion, and axial deformation. The element is a 15-DOF Bernoulli-Euler beam with cubic-Hermite shape functions on the bending DOFs and linear shape functions on the axial and torsion DOFs.

Discrete eigenproblem

After assembly + boundary-condition reduction the generalised eigenproblem is:

\[\bigl[\, \mathbf{K} - \omega^2\, \mathbf{M} \,\bigr]\, \boldsymbol{\varphi} \;=\; \mathbf{0}\]

with \(\mathbf{K}\) the global stiffness matrix (elastic + geometric, including centrifugal stiffening on a rotating blade), \(\mathbf{M}\) the global mass matrix (including added mass and parallel-axis contributions for a floating platform), \(\boldsymbol{\varphi}\) the mode-shape eigenvector, and \(\omega\) the natural circular frequency. pyBmodes reports frequencies in Hz: \(f = \omega / (2\pi)\).

Element matrices

The element-level stiffness and mass matrices are built from the distributed section properties:

  • mass density \(\rho A(s)\) (kg / m)

  • axial stiffness \(EA(s)\) (N)

  • bending stiffness \(EI_\text{flap}(s)\), \(EI_\text{edge}(s)\) (N · m²)

  • torsional stiffness \(GJ(s)\) (N · m²)

with \(s\) the dimensionless span coordinate. Mass and elastic-centre offsets, pre-twist, and concentrated tip mass are added through the same element framework.

Element matrix assembly is vectorised over Gauss points and over elements via numpy.einsum — the original per-element loop became a single tensor contraction with ~ 2–3× speedup on representative towers without changing the numerics.

Sign convention + normalisation

Mode shapes are mass-normalised:

\[\boldsymbol{\varphi}_i^{\mathrm{T}}\, \mathbf{M}\, \boldsymbol{\varphi}_i \;=\; 1 \qquad \forall\, i\]

The sign of each mode shape is canonicalised by the rule “maximum-amplitude DOF is positive”, so two independent solves of the same problem produce sign-stable shapes that can be MAC- compared without a sign-flip ambiguity.

Boundary conditions

Four hub_conn modes are supported. Each is cross-verified against either the BModes Fortran reference solver or an analytical closed form:

hub_conn

Meaning

Verified against

Worst-case error

1

Cantilever. All six base DOFs locked.

BModes CertTest Test03 / Test04

< 0.005 %

2

Free-free. All six base DOFs released; reactions through the PlatformSupport block (6 × 6 mooring + hydrodynamic + inertial matrices).

OC3 Hywind floating spar (Jonkman 2010)

≤ 0.0003 % across the first 9 modes

3

Soft monopile. Axial + torsion locked; lateral + rocking free. Optional distributed Winkler soil stiffness along the embedded section.

CS_Monopile (Jonkman & Musial 2010)

< 0.005 % across 10 modes

4

Pinned-free cable. Axial + transverse deflections + twist locked; bending slopes free — Bir 2009’s cable BC.

Bir 2009 Eq. 8 (analytical Legendre polynomial)

< 0.5 % across flap modes 1–3, Ω ∈ {2..30} rad/s

See Validation matrix for the per-case test-file references.

Platform-attached frame (hub_conn = 2)

For hub_conn = 2 the tower lives in the platform-attached frame: the rigid-arm transform brings PlatformSupport’s 6 × 6 mass and stiffness matrices to the tower base before assembly. The structural-inertia transform uses the full 3-D arm r = (PtfmCMxt, PtfmCMyt, cm_pform − draft); the hydro and mooring matrices are referenced at the platform reference point (PtfmRefxt = PtfmRefyt = 0 for every standard HydroDyn / WAMIT deck), so a zero horizontal arm applies to those — see pybmodes.fem.nondim._rigid_arm_T().

Asymmetric platform support routes the eigenproblem through scipy.linalg.eig() (general dense) instead of scipy.linalg.eigh() (symmetric) — see Solver dispatch below. The OC3 Hywind cross-coupled hydro_K + mooring_K exercises this exact branch and is the benchmark behind the 0.0003 % regression.

ElastoDyn-compatible polynomial ansatz

ElastoDyn represents tower and blade mode shapes as a constrained 6th-order polynomial in the dimensionless span coordinate \(s = h / H\):

\[\mathrm{SHP}(s) \;=\; \sum_{i=2}^{6} c_i\, s^{\,i}\]

The constraint \(i \geq 2\) enforces \(\mathrm{SHP}(0) = \mathrm{SHP}'(0) = 0\) algebraically — a clamped-base condition. pybmodes fits FEM mode shapes to this ansatz under a least-squares constraint with design-matrix condition-number reporting. Public entry points:

  • pybmodes.elastodyn.compute_blade_params()

  • pybmodes.elastodyn.compute_tower_params()

  • pybmodes.elastodyn.compute_tower_params_report() — same plus a TowerSelectionReport exposing FA / SS family scoring and rejected-mode lists.

Two non-obvious consequences shape every floating workflow:

Floating-deck polynomials use the cantilever basis

The polynomial form can only express clamped-base modes. A floating-deck workflow that wants ElastoDyn-compatible polynomials must therefore use Tower.from_elastodyn(...) (hub_conn = 1), not Tower.from_bmi(...) with hub_conn = 2.

The audit behind this — OpenFAST source-code citations showing that ElastoDyn integrates only the tower beam plus TwrTpMass (no platform / hydro / mooring matrices) and adds platform 6-DOF motion at runtime via the rigid-body sum — lives in cases/ECOSYSTEM_FINDING.md and the FLOATING_CASES.md next to each bundled reference deck.

Torsion-contamination filter

Some FEM modes are hybrid bending + twist; the polynomial form cannot represent twist content. _select_tower_family applies a two-gate filter:

  1. Polynomial-fit residual gate. Candidates whose clamped-base polynomial fit has rms_residual <= 0.09 survive.

  2. Torsion-contamination gate. Candidates whose modal- kinetic-energy torsion fraction \(T_\text{tor} = \sum \varphi_\text{tor}^2 / \sum \varphi_\text{total}^2\) is \(\geq 0.10\) are dropped.

Per-mode \((T_\text{FA}, T_\text{SS}, T_\text{tor})\) participations travel through TowerSelectionReport.rejected_fa_modes / rejected_ss_modes and CoeffBlockResult.rejected_modes so the user sees what was dropped.

Improved Direct Method (root rigid-body subtraction)

For a free-free floating tower the FEM mode shape carries non-zero deflection and slope at the root (rigid-body platform motion). The polynomial form can’t represent this directly; pybmodes subtracts the tangent line at the root — \(\widetilde{y}(x) = y(x) - y(0) - y'(0) \cdot x\) — before fitting. That’s BModes’ Improved Direct Method (Bir 2010 §3); the residual is what the polynomial sees. The alternative Projection Method (rotate by \(-\arctan(\text{slope} \cdot \text{scale})\)) is not implemented — for small base slopes the two are equivalent, and the validation matrix shows the IDM suffices.

Solver dispatch

pybmodes chooses between three SciPy eigensolvers based on problem size + structure:

Solver

When used

Why

scipy.linalg.eigh()

Small symmetric problems (default)

Stable + cheap; returns eigenvalues sorted and orthonormal eigenvectors.

scipy.sparse.linalg.eigsh() (shift-invert)

Symmetric, ngd > 500, subset of modes requested

5–18× faster than dense eigh on real-tower meshes. mode='normal' with sigma=0 (not 'buckling', which degenerates to OP = K⁻¹ K = I there).

scipy.linalg.eig()

Genuinely asymmetric \(\mathbf{K}\) (after platform- support assembly)

Cross-coupled hydro_K + mooring_K on OC3 Hywind is not symmetric. Symmetrising would bias the platform rigid-body modes; the general dense path matches BModes.

The benchmark (scripts/benchmark_sparse_solver.py) reports speedups at n_elements ∈ {20, 50, 100, 200, 500} and asserts sparse beats dense above 100 elements.

Pre-solve sanity checks

pybmodes.checks.check_model() runs nine cheap, deterministic gates before a solve:

  1. Section properties all finite (no NaN / ±Inf) — runs first so the per-field checks below need not be NaN-aware

  2. Span monotonicity (no out-of-order nodes)

  3. Mass non-negativity at every node

  4. Stiffness-jump detection (> 5× between adjacent stations)

  5. FA / SS bending-ratio sanity (inside [0.1, 10])

  6. RNA mass vs tower mass ratio (warn when RNA > tower)

  7. Support-matrix singularity (for free-free without PlatformSupport)

  8. n_modes exceeding the DOF cap

  9. Polynomial-fit design-matrix condition number > 1e4 (WARN) / > 1e6 (FAIL)

Auto-runs in .run() (suppress with check_model=False); WARN and ERROR findings route through UserWarning. INFO findings are suppressed on the auto-path — they’re contextual, not actionable, and noisy at scale — but check_model(model) called directly surfaces every gate.

Mooring + hydrodynamics

Quasi-static elastic catenary

pybmodes.mooring.MooringSystem implements Jonkman 2007 Appendix B:

  • B-1 / B-2 for the fully-suspended branch (line entirely off the seabed).

  • B-7 / B-8 with \(C_B = 0\) for the anchor-on-seabed branch (one end resting on the bottom).

The non-linear catenary \((H, V_F)\) system is solved by a damped Newton iteration with an analytical 2 × 2 Jacobian. MooringSystem.stiffness_matrix(body_r6=None) returns the 6 × 6 linearised stiffness in OpenFAST DOF order [surge, sway, heave, roll, pitch, yaw]. The OC3 Hywind cross-coupling sign convention is pinned against Jonkman 2010 Table 5-1 by tests.test_mooring.test_oc3hywind_bmi_dof_order_matches_jonkman_2010().

WAMIT / HydroDyn

pybmodes.io.wamit_reader.HydroDynReader follows the PotFile pointer in a HydroDyn .dat to the WAMIT .1 / .hst outputs and returns a WamitData with SI A_inf (infinite-frequency added mass), A_0 (zero-frequency added mass), and C_hst (hydrostatic restoring) 6 × 6 matrices.

WAMIT’s v7 re-dimensionalisation (ρ · L^k / ρ · g · L^k per DOF-pair type, \(k \in \{2..5\}\)) is handled by the reader; Fortran-style D / d exponent notation is normalised before parsing; upper-triangle-only outputs are mirrored after read.

Citable references

The reference set pybmodes validates against (NREL reports link to the canonical docs.nrel.gov PDF; papers by DOI; textbooks have no DOI):

  • NREL 5MW Reference Turbine — Jonkman, Butterfield, Musial, Scott (2009), Definition of a 5-MW Reference Wind Turbine for Offshore System Development, NREL/TP-500-38060 (PDF).

  • OC3 Monopile and OC3 Hywind (floating spar) — Jonkman & Musial (2010), NREL/TP-5000-48191 (PDF); Jonkman (2010), Definition of the Floating System for Phase IV of OC3, NREL/TP-500-47535 (PDF).

  • IEA-3.4-130-RWT, IEA-10-198-RWT, IEA-15-240-RWT, IEA-22-280-RWT — Bortolotti, Tarrés, Dykes et al. (2019), IEA Wind TCP Task 37: Systems Engineering in Wind Energy — WP2.1 Reference Wind Turbines, NREL/TP-5000-73492 (PDF), and the follow-on IEA Wind Task 37 reports for the larger sizes — notably the IEA-15-240-RWT in Gaertner et al. (2020), NREL/TP-5000-75698 (PDF), on the UMaine VolturnUS-S platform of Allen et al. (2020), NREL/TP-5000-76773 (PDF).

  • BModes — Bir (2005), User’s Guide to BModes, NREL/TP-500-39133 (PDF); verification in Bir (2010), NREL/CP-500-47953 (PDF).

  • Rotating-blade closed forms — Wright (1982), J. Appl. Mech. 49(1), 197–202 (DOI 10.1115/1.3161966); Bir (2009), AIAA 2009-1035 (DOI 10.2514/6.2009-1035); Bir (2010) Table 5 (tip-mass rotating blade); Bir 2009 Eq. 8 (rotating cable, analytical Legendre solution).

  • Beam tip-mass formulas — Blevins (1979 / 2016); Karnovsky & Lebed (2001). (Textbooks; no DOI.)

  • Catenary mooring — Jonkman (2007), Dynamics Modeling and Loads Analysis of an Offshore Floating Wind Turbine, NREL/TP-500-41958, Appendix B (PDF).

  • WAMIT — Lee & Newman (1991/2006), WAMIT User Manual.

  • Cable structures — Irvine (1981), Cable Structures, MIT Press, §2.4.

See Validation matrix for the per-case mapping of reference → quantity-being-checked → test file → worst-observed margin.