From c5655fff7d8b4555f74efcc1d2d8664c6873974d Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sun, 5 Jul 2026 22:05:39 +0900 Subject: [PATCH 01/17] fix: include rotor inertia from spanwise blade mass in auto-RNA (#130) read_windio_rna lumped the blades as a bare point mass at the rotor apex, which captured the hub-to-tower-top translation but dropped the rotor's own diametral inertia from the blade mass spread along the span. On IEA-22 that dropped term is ~2.5x the tower-top rotary inertia the point-mass lump produced, so rigid-RNA tower FA/SS frequencies came out too high versus a rigid-rotor reference (thanks @TheMercer). Assemble the rotor as a rigid body instead: diag([I_polar, I_polar/2, I_polar/2]) about the hub, with I_polar = N_bl * integral (dm/ds) * r^2 ds and r = hub_radius + span*cos(cone) (hub diameter and cone angle used when present, else 0 / no cone). Only each blade's own sectional spin inertia is still excluded. Mass and CM are unchanged; the inertia (and coupled tower frequencies) change relative to 1.16.0. Updates the hand-computed assembly test and adds a rotor-inertia test (hub-radius and span-length sensitivity). --- CHANGELOG.md | 20 ++++++++++- src/pybmodes/io/windio.py | 76 +++++++++++++++++++++++++++++++-------- tests/test_windio_rna.py | 48 ++++++++++++++++++++++--- 3 files changed, 125 insertions(+), 19 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 64537ac..9f998d8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,25 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). ## [Unreleased] -(nothing yet) +### Fixed + +- **WindIO auto-RNA now carries the rotor inertia from the spanwise blade + mass (#130).** `read_windio_rna` (and `lumped_rna_cal=True`) previously + lumped the blades as a bare point mass at the rotor apex, which captured + the hub-to-tower-top translation but dropped the rotor's own diametral + inertia from the blade mass being spread along the span. That term is the + dominant part of the rotor's contribution to the tower-top rotary + inertia (on IEA-22 it is ~2.5x the value the point-mass lump produced), + so tower fore-aft / side-side frequencies from a rigid-RNA lump came out + too high versus a rigid-rotor reference. The rotor is now assembled as a + rigid body, `diag([I_polar, I_polar/2, I_polar/2])` about the hub with + `I_polar = N_bl · ∫ (dm/ds) · r² ds` (`r = hub_radius + span·cos(cone)`), + using the hub diameter and cone angle when present. Only each blade's own + sectional spin inertia is still excluded. This changes the auto-RNA + inertia (and hence coupled tower frequencies) relative to 1.16.0; mass + and centre of mass are unchanged. This intentionally goes beyond the + ElastoDyn deck path's point-mass lumping, which the WindIO per-station + blade mass makes possible. ## [1.16.0] — 2026-07-03 diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 21741cc..65d1af1 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -783,8 +783,10 @@ def _blade_reference_axis(six: dict, comp: dict, blade_component: str) -> dict: ) -def _blade_single_mass(six: dict, ref: dict, blade_component: str) -> float: - """Integrate a WindIO blade's mass per unit length over its span. +def _blade_span_mass_inertia( + six: dict, ref: dict, blade_component: str, hub_r: float, cone_cos: float, +) -> tuple[float, float]: + """Integrate a WindIO blade's span into ``(mass, polar_second_moment)``. ``six`` is ``components..elastic_properties_mb.six_x_six`` and ``ref`` is @@ -794,6 +796,11 @@ def _blade_single_mass(six: dict, ref: dict, blade_component: str) -> float: ``inertia_matrix.values``, integrated over the arc length of the reference axis. Only the ``z`` curve is required; ``x`` / ``y`` (prebend / sweep) are optional and default to zero (a straight span). + + Returns the single-blade mass and the single-blade polar second moment + about the rotor axis, ``∫ (dm/ds) · r(s)² ds`` with the radial distance + ``r(s) = hub_r + z(s)·cone_cos``, so the caller can build the rotor's + diametral inertia from the spanwise mass distribution (issue #130). """ im = _require_mapping(six.get("inertia_matrix"), "blade six_x_six.inertia_matrix") grid = np.asarray(_require_key(im, "grid", "blade inertia_matrix"), dtype=float) @@ -864,7 +871,12 @@ def _blade_single_mass(six: dict, ref: dict, blade_component: str) -> float: xyz = np.vstack(coords).T seg = np.linalg.norm(np.diff(xyz, axis=0), axis=1) s = np.concatenate([[0.0], np.cumsum(seg)]) - return _trapezoid(mpl, s) + mass = _trapezoid(mpl, s) + # Radial distance from the rotor axis: the hub radius plus the spanwise + # position projected onto the rotor plane by the cone angle. + radial = hub_r + coords[2] * cone_cos + polar_second_moment = _trapezoid(mpl * radial * radial, s) + return mass, polar_second_moment def read_windio_rna( @@ -892,11 +904,11 @@ def read_windio_rna( Frame: tower-top ``x = downwind, y = lateral, z = up``. The nacelle inertia is the ontology's own tensor about the nacelle CM; the hub is - its tensor about the hub centre; the blades are a point mass at the - rotor apex (matching the ElastoDyn rigid-RNA convention). The result is - expressed at the tower top with ``cm_offset = 0`` and the vertical CM - lever in ``cm_axial`` so the FEM nondimensionaliser does not re-apply - parallel-axis. + its tensor about the hub centre; the rotor is a rigid body at the apex + (total blade mass plus the diametral inertia from the spanwise blade + mass). The result is expressed at the tower top with ``cm_offset = 0`` + and the vertical CM lever in ``cm_axial`` so the FEM nondimensionaliser + does not re-apply parallel-axis. ``cm_offset`` is intentionally ``0``, not the horizontal CM / overhang. The auto-RNA is consumed only by the clamped-base (``hub_conn = 1``) @@ -906,10 +918,16 @@ def read_windio_rna( returned inertia tensor: the parallel-axis shift to the tower top puts the ``m·overhang²`` and ``m·height²`` terms into ``ixx`` / ``iyy`` / ``izz``. Setting ``cm_offset = cm[0]`` on top of the already-shifted - tensor would double-count the overhang. Likewise the blades' own - spanwise spin inertia is left at zero (their parallel-axis contribution - is still in the tensor via the point mass at the apex) — the ElastoDyn - convention that reproduces the Jonkman (2009) tower frequencies. + tensor would double-count the overhang. + + Blade inertia (issue #130): the rotor carries the diametral inertia + ``N_bl · ∫ (dm/ds) · r² ds`` from the blade mass spread along the span + (``r = hub_radius + span·cos(cone)``), as ``diag([I_polar, I_polar/2, + I_polar/2])`` about the hub (the same perpendicular-axis form as the + hub tensor). Only each blade's own *sectional* spin inertia is left out + (its parallel-axis / span contribution is captured). This intentionally + goes beyond the ElastoDyn deck path's bare point-mass lumping, which the + WindIO ontology's per-station blade mass makes possible. Requires the optional ``[windio]`` extra (PyYAML). """ @@ -989,6 +1007,21 @@ def _nac_geom(key: str) -> float: _require_key(hub_ep, "system_inertia", "hub elastic_properties_mb"), "hub", ) + # Rotor geometry for the blade-span inertia (issue #130): the hub radius + # offsets the blade root from the rotor axis, and the cone angle + # projects the coned span onto the rotor plane. Both optional (default + # no hub offset / no cone). + hub_diam = hub.get("diameter") + hub_r = ( + 0.5 * _finite_float(hub_diam, "hub diameter") if hub_diam is not None else 0.0 + ) + if hub_r < 0.0: + raise ValueError(f"hub diameter must be non-negative; got {hub_diam!r}.") + cone = hub.get("cone_angle") + cone_cos = ( + float(np.cos(_finite_float(cone, "hub cone_angle"))) + if cone is not None else 1.0 + ) # --- Blades: assembly.number_of_blades x integrated span mass --- n_blades = _require_key(assembly, "number_of_blades", "assembly") @@ -1004,6 +1037,7 @@ def _nac_geom(key: str) -> float: f"{assembly.get('rotor_orientation')!r}." ) m_blade_each = 0.0 + i_polar_each = 0.0 if n_blades > 0: blade = _require_mapping( comps.get(component_blade), f"components.{component_blade}" @@ -1017,9 +1051,23 @@ def _nac_geom(key: str) -> float: f"components.{component_blade}.elastic_properties_mb.six_x_six", ) ref = _blade_reference_axis(six, blade, component_blade) - m_blade_each = _blade_single_mass(six, ref, component_blade) + m_blade_each, i_polar_each = _blade_span_mass_inertia( + six, ref, component_blade, hub_r, cone_cos, + ) m_blades = n_blades * m_blade_each + # Rotor inertia from the spanwise blade mass (issue #130). For N >= 3 + # symmetric blades the rotor tensor about the hub is + # diag([I_polar, I_polar/2, I_polar/2]) in the shaft frame (polar about + # the shaft, half on each transverse axis by the perpendicular-axis + # theorem), the same form as the hub. Lumping the blades as a bare point + # mass at the apex would drop this, which is the dominant part of the + # rotor's contribution to the tower-top rotary inertia. + i_polar_rotor = n_blades * i_polar_each + i_blades = np.diag( + [i_polar_rotor, 0.5 * i_polar_rotor, 0.5 * i_polar_rotor] + ) + # Rotor apex relative to the tower top. An upwind rotor sits at negative # x (downwind-positive frame); the vertical hub position is # distance_tt_hub directly (= Twr2Shft + overhang*sin(uptilt)). @@ -1032,7 +1080,7 @@ def _nac_geom(key: str) -> float: bodies = [ (m_nac, r_nac, i_nac), (m_hub, apex, i_hub), - (m_blades, apex.copy(), np.zeros((3, 3))), + (m_blades, apex.copy(), i_blades), ] m_total = float(sum(m for m, _, _ in bodies)) diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index 88ebf61..c52166d 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -96,15 +96,55 @@ def test_rna_assembly_matches_hand_computation(tmp_path: pathlib.Path) -> None: # cm_z = (500000*3 + 100000*4 + 15000*4) / 615000 assert rna.cm_axial == pytest.approx(1.96e6 / 615000.0) assert rna.cm_offset == 0.0 - # Parallel-axis tensor at the tower top (hand-computed). - assert rna.ixx == pytest.approx(1.834e7) - assert rna.iyy == pytest.approx(3.884e7) - assert rna.izz == pytest.approx(3.35e7) + # Parallel-axis tensor at the tower top (hand-computed), plus the rotor + # diametral inertia from the spanwise blade mass (issue #130): + # I_polar_rotor = 3 · ∫100·z² dz over [0, 50] (trapezoid, 2 stations) + # = 3 · 6.25e6 = 1.875e7, added about the apex as + # diag([1.875e7, 9.375e6, 9.375e6]). + assert rna.ixx == pytest.approx(1.834e7 + 1.875e7) + assert rna.iyy == pytest.approx(3.884e7 + 0.9375e7) + assert rna.izz == pytest.approx(3.35e7 + 0.9375e7) assert rna.izx == pytest.approx(1.06e7) assert rna.ixy == pytest.approx(0.0, abs=1e-6) assert rna.iyz == pytest.approx(0.0, abs=1e-6) +def test_rna_rotor_inertia_from_span(tmp_path: pathlib.Path) -> None: + """The rotor diametral inertia from the spanwise blade mass is included, + and the hub radius increases it (issue #130). Concentrating the same + blade mass near the rotor axis (short span) carries far less rotary + inertia than spreading it along a long span.""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + base = read_windio_rna(_write(_base_ontology(), tmp_path, "base_r.yaml")) + # The pre-#130 value for izz was 3.35e7 (no blade inertia); it must now + # be strictly larger because the rotor diametral term is included. + assert base.izz > 3.35e7 + + # A hub radius offsets every blade section further from the rotor axis, + # so the polar / diametral inertia grows. + with_hub = _base_ontology() + with_hub["components"]["hub"]["diameter"] = 8.0 # hub_R = 4 m + wh = read_windio_rna(_write(with_hub, tmp_path, "hub_r.yaml")) + assert wh.izz > base.izz + assert wh.iyy > base.iyy + assert wh.mass == pytest.approx(base.mass) # mass unchanged + + # Spreading the same blade mass over a longer span carries much more + # rotary inertia than keeping it near the root. + longer = _base_ontology() + longer["components"]["blade"]["elastic_properties_mb"]["six_x_six"][ + "inertia_matrix" + ]["values"] = [[50.0] + [0.0] * 20, [50.0] + [0.0] * 20] # half kg/m ... + longer["components"]["blade"]["elastic_properties_mb"]["six_x_six"][ + "reference_axis" + ]["z"] = {"grid": [0.0, 1.0], "values": [0.0, 100.0]} # ... over 100 m + lg = read_windio_rna(_write(longer, tmp_path, "long_r.yaml")) + assert lg.mass == pytest.approx(base.mass) # same 5000 kg/blade + assert lg.izz > base.izz # longer span -> more inertia + + def test_rna_diagonal_inertia_accepted(tmp_path: pathlib.Path) -> None: """A diagonal 3-vector system_inertia [Ixx, Iyy, Izz] is accepted as a diagonal tensor, matching the equivalent zero-product 6-vector (Codex From bf63315f65112191618cd4d774a5591279a2fb0d Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sun, 5 Jul 2026 22:19:58 +0900 Subject: [PATCH 02/17] fix: rotate shaft-frame RNA inertia by the uptilt (#130 review) Codex flagged that the rotor (and hub) inertia tensors are in the shaft frame but were summed as tower-top tensors without rotating by the shaft tilt, so a tilted rotor lost the tensor izx product and misallocated the large polar term between ixx and izz. Rotate both shaft-frame tensors into the tower frame by the uptilt (shaft axis [sign_x*cos, 0, sin]) before the parallel-axis assembly; the nacelle tensor is already tower-frame and is left alone. Adds a tilt-rotation test (izx appears, trace preserved). --- src/pybmodes/io/windio.py | 14 ++++++++++++++ tests/test_windio_rna.py | 26 ++++++++++++++++++++++++++ 2 files changed, 40 insertions(+) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 65d1af1..6b8cb5f 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -1077,6 +1077,20 @@ def _nac_geom(key: str) -> float: dtype=float, ) + # The hub and rotor inertia tensors are in the shaft frame; rotate them + # into the tower-top frame by the shaft tilt before they are summed as + # tower-top tensors (a tilted rotor otherwise loses the tensor izx + # product and misallocates the large polar term between ixx and izz). + # The shaft unit axis is ``[sign_x*cos(uptilt), 0, sin(uptilt)]``; R maps + # the shaft x-axis onto it and fixes the lateral y-axis. The nacelle + # tensor is already in the tower frame and is left alone. + c_t, s_t = float(np.cos(uptilt)), float(np.sin(uptilt)) + r_tilt = np.array( + [[sign_x * c_t, 0.0, -s_t], [0.0, 1.0, 0.0], [s_t, 0.0, sign_x * c_t]] + ) + i_hub = r_tilt @ i_hub @ r_tilt.T + i_blades = r_tilt @ i_blades @ r_tilt.T + bodies = [ (m_nac, r_nac, i_nac), (m_hub, apex, i_hub), diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index c52166d..b1436b0 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -145,6 +145,32 @@ def test_rna_rotor_inertia_from_span(tmp_path: pathlib.Path) -> None: assert lg.izz > base.izz # longer span -> more inertia +def test_rna_rotor_inertia_tilt_rotation(tmp_path: pathlib.Path) -> None: + """The shaft-frame rotor / hub inertia is rotated by the uptilt into the + tower frame, so a tilted rotor gains an izx product while preserving the + inertia trace (Codex review on #130).""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + def _mk(uptilt: float): + o = _base_ontology() + dt = o["components"]["nacelle"]["drivetrain"] + dt["overhang"] = 0.0 # apex purely vertical + dt["uptilt"] = uptilt + dt["elastic_properties_mb"]["system_center_mass"] = [0.0, 0.0, 3.0] + return read_windio_rna(_write(o, tmp_path, f"tilt_{uptilt}.yaml")) + + flat = _mk(0.0) + tilt = _mk(0.3) # ~17 deg + # overhang 0 + CM x 0 -> all izx comes from the shaft-frame tensor + # rotation; a flat rotor/hub are diagonal (izx 0), a tilted one gains izx. + assert flat.izx == pytest.approx(0.0, abs=1.0) + assert abs(tilt.izx) > 1.0e6 + # the rotation preserves the inertia trace (apex identical, overhang 0). + assert (tilt.ixx + tilt.iyy + tilt.izz) == pytest.approx( + flat.ixx + flat.iyy + flat.izz, rel=1e-9) + + def test_rna_diagonal_inertia_accepted(tmp_path: pathlib.Path) -> None: """A diagonal 3-vector system_inertia [Ixx, Iyy, Izz] is accepted as a diagonal tensor, matching the equivalent zero-product 6-vector (Codex From bb5c7ba6b9ab43d2e1cc186f4f32b6a89e206627 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sun, 5 Jul 2026 22:33:05 +0900 Subject: [PATCH 03/17] fix: account for coned-rotor axial inertia in transverse RNA terms The rigid-rotor diametral moment used I_polar/2, which is the flat-disc (perpendicular-axis) result and only holds for a planar rotor. A coned rotor sits its section masses off the rotor plane by span*sin(cone), so each blade adds N_bl * integral(dm * (span*sin(cone))**2) to both transverse moments about the hub. _blade_span_mass_inertia now returns the axial second moment alongside the radial one, and read_windio_rna folds it into the diametral term. The radial reach itself already used span*cos(cone), so the two are consistent. Behaviour is unchanged for an unconed rotor (sin(cone) = 0). --- src/pybmodes/io/windio.py | 63 ++++++++++++++++++++++----------------- tests/test_windio_rna.py | 24 +++++++++++++++ 2 files changed, 59 insertions(+), 28 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 6b8cb5f..88e21bd 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -784,9 +784,10 @@ def _blade_reference_axis(six: dict, comp: dict, blade_component: str) -> dict: def _blade_span_mass_inertia( - six: dict, ref: dict, blade_component: str, hub_r: float, cone_cos: float, -) -> tuple[float, float]: - """Integrate a WindIO blade's span into ``(mass, polar_second_moment)``. + six: dict, ref: dict, blade_component: str, + hub_r: float, cone_cos: float, cone_sin: float, +) -> tuple[float, float, float]: + """Integrate a blade's span into ``(mass, polar, axial)`` second moments. ``six`` is ``components..elastic_properties_mb.six_x_six`` and ``ref`` is @@ -797,10 +798,12 @@ def _blade_span_mass_inertia( reference axis. Only the ``z`` curve is required; ``x`` / ``y`` (prebend / sweep) are optional and default to zero (a straight span). - Returns the single-blade mass and the single-blade polar second moment - about the rotor axis, ``∫ (dm/ds) · r(s)² ds`` with the radial distance - ``r(s) = hub_r + z(s)·cone_cos``, so the caller can build the rotor's - diametral inertia from the spanwise mass distribution (issue #130). + Returns the single-blade mass, the polar second moment about the rotor + axis ``∫ (dm/ds) · r(s)² ds`` (radial ``r = hub_r + z·cone_cos``), and + the axial second moment ``∫ (dm/ds) · a(s)² ds`` (axial offset + ``a = z·cone_sin``). The caller builds the rotor's polar inertia from + the former and the coned transverse (diametral) inertia + ``polar/2 + N_bl·axial`` from both (issue #130). """ im = _require_mapping(six.get("inertia_matrix"), "blade six_x_six.inertia_matrix") grid = np.asarray(_require_key(im, "grid", "blade inertia_matrix"), dtype=float) @@ -872,11 +875,16 @@ def _blade_span_mass_inertia( seg = np.linalg.norm(np.diff(xyz, axis=0), axis=1) s = np.concatenate([[0.0], np.cumsum(seg)]) mass = _trapezoid(mpl, s) - # Radial distance from the rotor axis: the hub radius plus the spanwise - # position projected onto the rotor plane by the cone angle. - radial = hub_r + coords[2] * cone_cos + # The spanwise position resolves into a radial distance from the rotor + # axis (hub radius + span*cos(cone)) and, for a coned rotor, an axial + # offset along the shaft (span*sin(cone)). The radial part gives the + # polar second moment; the axial part adds to the transverse moment. + span = coords[2] + radial = hub_r + span * cone_cos + axial = span * cone_sin polar_second_moment = _trapezoid(mpl * radial * radial, s) - return mass, polar_second_moment + axial_second_moment = _trapezoid(mpl * axial * axial, s) + return mass, polar_second_moment, axial_second_moment def read_windio_rna( @@ -1017,11 +1025,9 @@ def _nac_geom(key: str) -> float: ) if hub_r < 0.0: raise ValueError(f"hub diameter must be non-negative; got {hub_diam!r}.") - cone = hub.get("cone_angle") - cone_cos = ( - float(np.cos(_finite_float(cone, "hub cone_angle"))) - if cone is not None else 1.0 - ) + cone_ang = _finite_float(hub.get("cone_angle", 0.0), "hub cone_angle") + cone_cos = float(np.cos(cone_ang)) + cone_sin = float(np.sin(cone_ang)) # --- Blades: assembly.number_of_blades x integrated span mass --- n_blades = _require_key(assembly, "number_of_blades", "assembly") @@ -1038,6 +1044,7 @@ def _nac_geom(key: str) -> float: ) m_blade_each = 0.0 i_polar_each = 0.0 + i_axial_each = 0.0 if n_blades > 0: blade = _require_mapping( comps.get(component_blade), f"components.{component_blade}" @@ -1051,22 +1058,22 @@ def _nac_geom(key: str) -> float: f"components.{component_blade}.elastic_properties_mb.six_x_six", ) ref = _blade_reference_axis(six, blade, component_blade) - m_blade_each, i_polar_each = _blade_span_mass_inertia( - six, ref, component_blade, hub_r, cone_cos, + m_blade_each, i_polar_each, i_axial_each = _blade_span_mass_inertia( + six, ref, component_blade, hub_r, cone_cos, cone_sin, ) m_blades = n_blades * m_blade_each - # Rotor inertia from the spanwise blade mass (issue #130). For N >= 3 - # symmetric blades the rotor tensor about the hub is - # diag([I_polar, I_polar/2, I_polar/2]) in the shaft frame (polar about - # the shaft, half on each transverse axis by the perpendicular-axis - # theorem), the same form as the hub. Lumping the blades as a bare point - # mass at the apex would drop this, which is the dominant part of the - # rotor's contribution to the tower-top rotary inertia. + # Rotor inertia from the spanwise blade mass (issue #130). In the shaft + # frame the rotor tensor about the hub is diag([I_polar, I_diam, I_diam]): + # the polar term ``I_polar = N_bl·∫dm·r²`` (radial), and the transverse + # (diametral) term ``I_diam = I_polar/2 + N_bl·∫dm·a²``. The ``I_polar/2`` + # is the perpendicular-axis split for a planar rotor; the ``N_bl·∫dm·a²`` + # adds the coned rotor's shaft-axis (``a = span·sin(cone)``) offset, which + # a flat-disc assumption would drop. Lumping the blades as a bare point + # mass at the apex would drop the whole tensor. i_polar_rotor = n_blades * i_polar_each - i_blades = np.diag( - [i_polar_rotor, 0.5 * i_polar_rotor, 0.5 * i_polar_rotor] - ) + i_diam_rotor = 0.5 * i_polar_rotor + n_blades * i_axial_each + i_blades = np.diag([i_polar_rotor, i_diam_rotor, i_diam_rotor]) # Rotor apex relative to the tower top. An upwind rotor sits at negative # x (downwind-positive frame); the vertical hub position is diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index b1436b0..a781edb 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -171,6 +171,30 @@ def _mk(uptilt: float): flat.ixx + flat.iyy + flat.izz, rel=1e-9) +def test_rna_rotor_inertia_cone_axial(tmp_path: pathlib.Path) -> None: + """A coned rotor adds the axial (span·sin(cone))² term to the transverse + (diametral) moment, so the transverse inertia grows with cone even as + the radial reach shrinks (Codex review on #130).""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + def _mk(cone: float): + o = _base_ontology() + o["components"]["hub"]["cone_angle"] = cone + dt = o["components"]["nacelle"]["drivetrain"] + dt["uptilt"] = 0.0 # no shaft-frame rotation + dt["overhang"] = 0.0 # apex purely vertical + dt["elastic_properties_mb"]["system_center_mass"] = [0.0, 0.0, 3.0] + return read_windio_rna(_write(o, tmp_path, f"cone_{cone}.yaml")) + + flat = _mk(0.0) + coned = _mk(0.4) # ~23 deg (exaggerated) + # I_diam = I_polar/2 + N_bl·axial: a flat-disc assumption (I_polar/2 only) + # would drop the axial term and give a smaller transverse moment. + assert coned.iyy > flat.iyy + assert coned.mass == pytest.approx(flat.mass) + + def test_rna_diagonal_inertia_accepted(tmp_path: pathlib.Path) -> None: """A diagonal 3-vector system_inertia [Ixx, Iyy, Izz] is accepted as a diagonal tensor, matching the equivalent zero-product 6-vector (Codex From 4be1d83137ec7a6728898a475ce83b792285ee21 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sun, 5 Jul 2026 23:00:11 +0900 Subject: [PATCH 04/17] fix: integrate spanwise rotor inertia with exact per-segment quadrature The second-moment integrals used the trapezoid rule, which is only exact for a linear integrand. mass*r**2 is cubic on each piecewise-linear segment (linear mass per length times a linear coordinate squared), so trapezoid chords over it and over-estimates on a coarse grid: a straight two-station 0-50 m blade came out 6.25e6 instead of the exact 4.17e6 per blade, a 50 % bias on the added rotor inertia. Replace it with _second_moment, a Simpson evaluation (endpoints plus the segment midpoint, which is the endpoint average for a linear coordinate). Simpson is exact through degree three, so it integrates the cubic without error. The mass integral keeps the trapezoid rule, which is exact for the linear mass-per-length integrand. --- src/pybmodes/io/windio.py | 30 ++++++++++++++++++++++++++++-- tests/test_windio_rna.py | 17 +++++++++++------ 2 files changed, 39 insertions(+), 8 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 88e21bd..d85f282 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -67,6 +67,32 @@ def _trapezoid(y: np.ndarray, x: np.ndarray) -> float: return float(np.sum(0.5 * (y[1:] + y[:-1]) * np.diff(x))) +def _second_moment(w: np.ndarray, c: np.ndarray, x: np.ndarray) -> float: + """Exact ``integral w(x)*c(x)**2 dx`` for piecewise-linear ``w`` and ``c``. + + On each ``[x_i, x_{i+1}]`` segment both the weight ``w`` (mass per unit + length) and the coordinate ``c`` vary linearly, so the integrand is a + cubic. Simpson's rule is exact through degree three, so evaluating the + endpoints plus the segment midpoint (``w`` / ``c`` averaged, since both + are linear) integrates it without error. The trapezoid rule chords over + the quadratic ``c**2`` and inflates the second moment on a coarse grid + (Codex review on #130 measured a 50 % bias for a straight two-station + blade). + """ + w = np.asarray(w, dtype=float) + c = np.asarray(c, dtype=float) + x = np.asarray(x, dtype=float) + dx = np.diff(x) + w0, w1 = w[:-1], w[1:] + c0, c1 = c[:-1], c[1:] + wm = 0.5 * (w0 + w1) + cm = 0.5 * (c0 + c1) + f0 = w0 * c0 * c0 + fm = wm * cm * cm + f1 = w1 * c1 * c1 + return float(np.sum(dx / 6.0 * (f0 + 4.0 * fm + f1))) + + def _require_yaml(): """Import PyYAML or raise the documented friendly error. @@ -882,8 +908,8 @@ def _blade_span_mass_inertia( span = coords[2] radial = hub_r + span * cone_cos axial = span * cone_sin - polar_second_moment = _trapezoid(mpl * radial * radial, s) - axial_second_moment = _trapezoid(mpl * axial * axial, s) + polar_second_moment = _second_moment(mpl, radial, s) + axial_second_moment = _second_moment(mpl, axial, s) return mass, polar_second_moment, axial_second_moment diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index a781edb..be2bedc 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -98,12 +98,13 @@ def test_rna_assembly_matches_hand_computation(tmp_path: pathlib.Path) -> None: assert rna.cm_offset == 0.0 # Parallel-axis tensor at the tower top (hand-computed), plus the rotor # diametral inertia from the spanwise blade mass (issue #130): - # I_polar_rotor = 3 · ∫100·z² dz over [0, 50] (trapezoid, 2 stations) - # = 3 · 6.25e6 = 1.875e7, added about the apex as - # diag([1.875e7, 9.375e6, 9.375e6]). - assert rna.ixx == pytest.approx(1.834e7 + 1.875e7) - assert rna.iyy == pytest.approx(3.884e7 + 0.9375e7) - assert rna.izz == pytest.approx(3.35e7 + 0.9375e7) + # I_polar_rotor = 3 · ∫100·z² dz over [0, 50] = 3 · (100·50³/3) + # = 3 · 4.1667e6 = 1.25e7, added about the apex as + # diag([1.25e7, 6.25e6, 6.25e6]). The exact span integral (Simpson, + # exact for the cubic mass·z² integrand) replaces the earlier trapezoid. + assert rna.ixx == pytest.approx(1.834e7 + 1.25e7) + assert rna.iyy == pytest.approx(3.884e7 + 0.625e7) + assert rna.izz == pytest.approx(3.35e7 + 0.625e7) assert rna.izx == pytest.approx(1.06e7) assert rna.ixy == pytest.approx(0.0, abs=1e-6) assert rna.iyz == pytest.approx(0.0, abs=1e-6) @@ -121,6 +122,10 @@ def test_rna_rotor_inertia_from_span(tmp_path: pathlib.Path) -> None: # The pre-#130 value for izz was 3.35e7 (no blade inertia); it must now # be strictly larger because the rotor diametral term is included. assert base.izz > 3.35e7 + # The diametral term is the exact span integral 3·(100·50³/3)/2 = 6.25e6, + # not the trapezoid over-estimate (9.375e6) a coarse two-station grid + # would give for the cubic mass·z² integrand (Codex review on #130). + assert base.izz == pytest.approx(3.35e7 + 6.25e6) # A hub radius offsets every blade section further from the rotor axis, # so the polar / diametral inertia grows. From 1d89bf1ca16aba9e77a46c1c590008bb4b855f92 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sun, 5 Jul 2026 23:21:39 +0900 Subject: [PATCH 05/17] fix: place coned rotor at its true CM before the parallel-axis shift The parallel-axis theorem I_O = I_C + m(|r|^2 I - r outer r) only holds when the tensor is taken about the body centre of mass. The rotor tensor was formed about the hub apex, but a coned rotor's mass sits off the hub plane along the shaft by span*sin(cone), so its CM is offset from the apex. Using the apex both as the tensor reference and as the shift origin double-counts that offset and drops the cross terms for a rotor with real overhang and shaft tilt. Compute the axial first moment (exact per-segment Simpson), take the shaft-axis CM offset as first-moment / mass, reduce the axial second moment to that CM (a mass-weighted variance, always non-negative), and place the rotor body at apex + offset along the shaft axis. A flat rotor (sin(cone) = 0) leaves the offset and axial term zero, so the planar I_polar/2 split about the apex is unchanged. For a coned rotor this raises the tower-top vertical CM by the precone term real ElastoDyn also carries (~0.06 m on IEA-22); the deck adapter approximates the blades as a point mass at the apex and omits it, so the cross-check test now documents that ~1 % divergence. --- src/pybmodes/io/windio.py | 78 ++++++++++++++++++++++++++++++--------- tests/test_windio_rna.py | 46 +++++++++++++++++++++-- 2 files changed, 103 insertions(+), 21 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index d85f282..558dc8f 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -93,6 +93,24 @@ def _second_moment(w: np.ndarray, c: np.ndarray, x: np.ndarray) -> float: return float(np.sum(dx / 6.0 * (f0 + 4.0 * fm + f1))) +def _first_moment(w: np.ndarray, c: np.ndarray, x: np.ndarray) -> float: + """Exact ``integral w(x)*c(x) dx`` for piecewise-linear ``w`` and ``c``. + + The integrand is quadratic on each segment (a product of two linear + factors), which Simpson's rule integrates exactly. Used for the rotor's + axial first moment, i.e. its shaft-axis CM offset, so a coned rotor body + is placed at its centre of mass before the parallel-axis shift. + """ + w = np.asarray(w, dtype=float) + c = np.asarray(c, dtype=float) + x = np.asarray(x, dtype=float) + dx = np.diff(x) + w0, w1 = w[:-1], w[1:] + c0, c1 = c[:-1], c[1:] + fm = 0.5 * (w0 + w1) * 0.5 * (c0 + c1) + return float(np.sum(dx / 6.0 * (w0 * c0 + 4.0 * fm + w1 * c1))) + + def _require_yaml(): """Import PyYAML or raise the documented friendly error. @@ -812,7 +830,7 @@ def _blade_reference_axis(six: dict, comp: dict, blade_component: str) -> dict: def _blade_span_mass_inertia( six: dict, ref: dict, blade_component: str, hub_r: float, cone_cos: float, cone_sin: float, -) -> tuple[float, float, float]: +) -> tuple[float, float, float, float]: """Integrate a blade's span into ``(mass, polar, axial)`` second moments. ``six`` is @@ -825,11 +843,14 @@ def _blade_span_mass_inertia( (prebend / sweep) are optional and default to zero (a straight span). Returns the single-blade mass, the polar second moment about the rotor - axis ``∫ (dm/ds) · r(s)² ds`` (radial ``r = hub_r + z·cone_cos``), and - the axial second moment ``∫ (dm/ds) · a(s)² ds`` (axial offset - ``a = z·cone_sin``). The caller builds the rotor's polar inertia from - the former and the coned transverse (diametral) inertia - ``polar/2 + N_bl·axial`` from both (issue #130). + axis ``∫ (dm/ds) · r(s)² ds`` (radial ``r = hub_r + z·cone_cos``), the + axial second moment ``∫ (dm/ds) · a(s)² ds`` (axial offset + ``a = z·cone_sin``), and the axial first moment ``∫ (dm/ds) · a(s) ds``. + The caller builds the rotor's polar inertia from the second radial + moment, the coned transverse (diametral) inertia from the radial and + axial second moments, and places the rotor body at its coned centre of + mass (shaft-axis offset ``axial_first / mass``) before the parallel-axis + shift (issue #130). """ im = _require_mapping(six.get("inertia_matrix"), "blade six_x_six.inertia_matrix") grid = np.asarray(_require_key(im, "grid", "blade inertia_matrix"), dtype=float) @@ -910,7 +931,8 @@ def _blade_span_mass_inertia( axial = span * cone_sin polar_second_moment = _second_moment(mpl, radial, s) axial_second_moment = _second_moment(mpl, axial, s) - return mass, polar_second_moment, axial_second_moment + axial_first_moment = _first_moment(mpl, axial, s) + return mass, polar_second_moment, axial_second_moment, axial_first_moment def read_windio_rna( @@ -1071,6 +1093,7 @@ def _nac_geom(key: str) -> float: m_blade_each = 0.0 i_polar_each = 0.0 i_axial_each = 0.0 + m_axial1_each = 0.0 if n_blades > 0: blade = _require_mapping( comps.get(component_blade), f"components.{component_blade}" @@ -1084,21 +1107,35 @@ def _nac_geom(key: str) -> float: f"components.{component_blade}.elastic_properties_mb.six_x_six", ) ref = _blade_reference_axis(six, blade, component_blade) - m_blade_each, i_polar_each, i_axial_each = _blade_span_mass_inertia( + ( + m_blade_each, + i_polar_each, + i_axial_each, + m_axial1_each, + ) = _blade_span_mass_inertia( six, ref, component_blade, hub_r, cone_cos, cone_sin, ) m_blades = n_blades * m_blade_each # Rotor inertia from the spanwise blade mass (issue #130). In the shaft - # frame the rotor tensor about the hub is diag([I_polar, I_diam, I_diam]): - # the polar term ``I_polar = N_bl·∫dm·r²`` (radial), and the transverse - # (diametral) term ``I_diam = I_polar/2 + N_bl·∫dm·a²``. The ``I_polar/2`` - # is the perpendicular-axis split for a planar rotor; the ``N_bl·∫dm·a²`` - # adds the coned rotor's shaft-axis (``a = span·sin(cone)``) offset, which - # a flat-disc assumption would drop. Lumping the blades as a bare point - # mass at the apex would drop the whole tensor. + # frame the rotor tensor is diag([I_polar, I_diam, I_diam]): the polar + # term ``I_polar = N_bl·∫dm·r²`` (radial ``r = hub_r + span·cos(cone)``), + # and the transverse (diametral) term. For a coned rotor the blade mass + # sits off the hub plane by ``a = span·sin(cone)``, so the rotor CM shifts + # along the shaft axis by ``offset = ∫dm·a / m_blades`` and the tensor is + # formed about that CM: ``I_diam = I_polar/2 + N_bl·∫dm·a² − m_blades· + # offset²`` (the axial second moment reduced to the CM by the parallel- + # axis theorem, a mass-weighted variance that stays non-negative). The + # body is placed at that CM below so the shift to the tower top is exact. + # A flat rotor (``sin(cone) = 0``) leaves the offset and the axial term + # zero, recovering the planar ``I_polar/2`` split about the apex. Lumping + # the blades as a bare point mass at the apex would drop the whole tensor. i_polar_rotor = n_blades * i_polar_each - i_diam_rotor = 0.5 * i_polar_rotor + n_blades * i_axial_each + cm_axial_offset = (m_axial1_each / m_blade_each) if m_blade_each > 0.0 else 0.0 + i_axial_cm = max( + 0.0, n_blades * i_axial_each - m_blades * cm_axial_offset * cm_axial_offset + ) + i_diam_rotor = 0.5 * i_polar_rotor + i_axial_cm i_blades = np.diag([i_polar_rotor, i_diam_rotor, i_diam_rotor]) # Rotor apex relative to the tower top. An upwind rotor sits at negative @@ -1124,10 +1161,17 @@ def _nac_geom(key: str) -> float: i_hub = r_tilt @ i_hub @ r_tilt.T i_blades = r_tilt @ i_blades @ r_tilt.T + # The coned rotor CM sits ``cm_axial_offset`` beyond the apex along the + # shaft axis (r_tilt's first column, the outboard unit vector the overhang + # already runs along), so a precone displaces the blade mass further from + # the tower in the same sense for both rotor orientations. + shaft_axis = r_tilt[:, 0] + rotor_cm = apex + cm_axial_offset * shaft_axis + bodies = [ (m_nac, r_nac, i_nac), (m_hub, apex, i_hub), - (m_blades, apex.copy(), i_blades), + (m_blades, rotor_cm, i_blades), ] m_total = float(sum(m for m, _, _ in bodies)) diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index be2bedc..cae49f7 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -200,6 +200,35 @@ def _mk(cone: float): assert coned.mass == pytest.approx(flat.mass) +def test_rna_coned_cm_shift_with_uptilt(tmp_path: pathlib.Path) -> None: + """A coned rotor's CM sits off the hub apex along the (tilted) shaft, so + with uptilt the tower-top vertical CM shifts by m_blades·offset·sin(uptilt) + over the total mass. Pins that the rotor body is placed at its coned CM, + not the apex, before the parallel-axis shift (Codex review on #130).""" + pytest.importorskip("yaml") + import numpy as np + + from pybmodes.io.windio import read_windio_rna + + uptilt = 0.10 + + def _mk(cone: float): + o = _base_ontology() + o["components"]["hub"]["cone_angle"] = cone + o["components"]["nacelle"]["drivetrain"]["uptilt"] = uptilt + return read_windio_rna(_write(o, tmp_path, f"tc_{cone}.yaml")) + + flat = _mk(0.0) + coned = _mk(0.15) + # base blade: uniform 100 kg/m over 0..50 m -> mass-weighted mean span 25 m, + # so the coned shaft-axis CM offset is sin(cone)·25. + offset = float(np.sin(0.15)) * 25.0 + m_blades = 3 * 5000.0 + expected = m_blades * offset * float(np.sin(uptilt)) / flat.mass + assert coned.mass == pytest.approx(flat.mass) + assert coned.cm_axial - flat.cm_axial == pytest.approx(expected, rel=1e-6) + + def test_rna_diagonal_inertia_accepted(tmp_path: pathlib.Path) -> None: """A diagonal 3-vector system_inertia [Ixx, Iyy, Izz] is accepted as a diagonal tensor, matching the equivalent zero-product 6-vector (Codex @@ -558,9 +587,15 @@ def test_from_windio_with_monopile_lumped_rna_cal(tmp_path: pathlib.Path) -> Non reason="IEA-22 WindIO ontology + ElastoDyn deck not present", ) def test_rna_iea22_matches_elastodyn_deck() -> None: - """The WindIO auto-RNA mass and CM match the IEA-22 ElastoDyn deck's - tower-top assembly to <0.5 % (issue #82). Inertia is WindIO-native and - is not expected to byte-match the deck's NacYIner (ecosystem drift).""" + """The WindIO auto-RNA mass matches the IEA-22 ElastoDyn deck's tower-top + assembly to <0.5 % (issue #82), and its CM agrees to ~1 %. The small CM + divergence is a documented, intentional one: the deck adapter lumps the + blades as a point mass at the rotor apex, whereas the WindIO rigid-rotor + path (issue #130) places the coned rotor at its true centre of mass, + which sits ~2.6 m outboard of the apex for the 4 deg precone and so + raises the tower-top vertical CM by ~0.06 m through the 6 deg shaft tilt + (real ElastoDyn carries the same precone term). Inertia is WindIO-native + and is not expected to byte-match the deck's NacYIner (ecosystem drift).""" pytest.importorskip("yaml") from pybmodes.io._elastodyn.adapter import _tower_top_assembly_mass from pybmodes.io.elastodyn_reader import ( @@ -576,7 +611,10 @@ def test_rna_iea22_matches_elastodyn_deck() -> None: ) assert rna.mass == pytest.approx(deck.mass, rel=5e-3) # <0.5 % - assert rna.cm_axial == pytest.approx(deck.cm_axial, rel=5e-3) + # The precone CM offset the deck adapter drops raises the WindIO CM + # slightly; it stays within ~1.5 % and lies above the apex-lumped value. + assert rna.cm_axial == pytest.approx(deck.cm_axial, rel=1.5e-2) + assert rna.cm_axial >= deck.cm_axial # WindIO-native inertia carries the fore-aft rotary term the same order # of magnitude as the deck, but not byte-identical (documented drift). assert rna.iyy > 0.0 From 90296794d2f3321786411a29885c69201f6ff99c Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sun, 5 Jul 2026 23:32:03 +0900 Subject: [PATCH 06/17] fix: include prebend and sweep in the rotor polar lever The rotor polar moment took the section lever from the z span alone, so a swept or prebent blade was undercounted: the mass integration follows the full 3-D reference-axis arc, but the lever ignored the in-plane tangential sweep and the prebend projection. Resolve each section onto the rotor axes instead: in-plane radial r = hub_r + z*cos(cone) - x*sin(cone), tangential t = y, axial a = z*sin(cone) + x*cos(cone). The polar lever is then the full in-plane distance from the shaft axis r**2 + t**2, and the axial term carries the prebend projection too. A straight blade (x = y = 0) is unchanged. --- src/pybmodes/io/windio.py | 31 ++++++++++++++++++++----------- tests/test_windio_rna.py | 24 ++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 11 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 558dc8f..1f4662b 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -843,9 +843,11 @@ def _blade_span_mass_inertia( (prebend / sweep) are optional and default to zero (a straight span). Returns the single-blade mass, the polar second moment about the rotor - axis ``∫ (dm/ds) · r(s)² ds`` (radial ``r = hub_r + z·cone_cos``), the - axial second moment ``∫ (dm/ds) · a(s)² ds`` (axial offset - ``a = z·cone_sin``), and the axial first moment ``∫ (dm/ds) · a(s) ds``. + axis ``∫ (dm/ds) · (r² + t²) ds`` (in-plane radial + ``r = hub_r + z·cone_cos − x·cone_sin`` and tangential ``t = y`` so + prebend / sweep are included), the axial second moment + ``∫ (dm/ds) · a(s)² ds`` (axial offset ``a = z·cone_sin + x·cone_cos``), + and the axial first moment ``∫ (dm/ds) · a(s) ds``. The caller builds the rotor's polar inertia from the second radial moment, the coned transverse (diametral) inertia from the radial and axial second moments, and places the rotor body at its coned centre of @@ -922,14 +924,21 @@ def _blade_span_mass_inertia( seg = np.linalg.norm(np.diff(xyz, axis=0), axis=1) s = np.concatenate([[0.0], np.cumsum(seg)]) mass = _trapezoid(mpl, s) - # The spanwise position resolves into a radial distance from the rotor - # axis (hub radius + span*cos(cone)) and, for a coned rotor, an axial - # offset along the shaft (span*sin(cone)). The radial part gives the - # polar second moment; the axial part adds to the transverse moment. - span = coords[2] - radial = hub_r + span * cone_cos - axial = span * cone_sin - polar_second_moment = _second_moment(mpl, radial, s) + # Resolve each section onto the rotor axes (r = in-plane radial from the + # shaft, t = in-plane tangential, a = along the shaft). The reference axis + # gives z spanwise, x prebend (flapwise, out of the rotor plane) and y + # sweep (edgewise, in the rotor plane). Coning rotates the span/prebend + # (radial/axial) plane by the cone angle about the edgewise axis; sweep is + # unaffected. The polar lever is the in-plane distance from the shaft axis + # ``r² + t²`` (so sweep and prebend are not dropped, issue #130); the + # axial offset feeds the coned transverse term. + x_off, y_off, span = coords[0], coords[1], coords[2] + radial = hub_r + span * cone_cos - x_off * cone_sin + tangential = y_off + axial = span * cone_sin + x_off * cone_cos + polar_second_moment = _second_moment(mpl, radial, s) + _second_moment( + mpl, tangential, s + ) axial_second_moment = _second_moment(mpl, axial, s) axial_first_moment = _first_moment(mpl, axial, s) return mass, polar_second_moment, axial_second_moment, axial_first_moment diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index cae49f7..90b4f30 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -200,6 +200,30 @@ def _mk(cone: float): assert coned.mass == pytest.approx(flat.mass) +def test_rna_rotor_inertia_includes_sweep(tmp_path: pathlib.Path) -> None: + """Sweep (reference_axis.y) sets an in-plane tangential distance from the + shaft axis, so it adds ``y²`` to the rotor polar lever. A constant + tangential offset leaves the span arc length (and thus the blade mass) + unchanged but raises the rotor inertia; a lever built from the z span + alone would drop it (Codex review on #130).""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + def _mk(y_off: float): + o = _base_ontology() + ra = o["components"]["blade"]["elastic_properties_mb"]["six_x_six"][ + "reference_axis" + ] + ra["y"] = {"grid": [0.0, 1.0], "values": [y_off, y_off]} + return read_windio_rna(_write(o, tmp_path, f"sw_{y_off}.yaml")) + + straight = _mk(0.0) + swept = _mk(10.0) + assert swept.mass == pytest.approx(straight.mass) # constant offset, same arc + assert swept.ixx > straight.ixx + assert swept.izz > straight.izz + + def test_rna_coned_cm_shift_with_uptilt(tmp_path: pathlib.Path) -> None: """A coned rotor's CM sits off the hub apex along the (tilted) shaft, so with uptilt the tower-top vertical CM shifts by m_blades·offset·sin(uptilt) From bf23bbe35afdc14cdc760f5f0879479d315d5034 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sun, 5 Jul 2026 23:40:39 +0900 Subject: [PATCH 07/17] fix: reject one- and two-bladed rotors in the auto-RNA inertia path The rotor tensor uses the axisymmetric transverse split I_diam = I_polar/2 + N*integral(dm*a**2), which is exact only for three or more evenly spaced blades, whose mass ring has isotropic in-plane inertia. A one- or two- bladed rotor lies on a single diameter, so its transverse inertia is azimuth-dependent and no single static tower-top lump represents it (one in-plane axis carries almost all of I_polar, the orthogonal one almost none). Rather than emit a corrupted fore-aft / side-side inertia, reject those blade counts with a message pointing at an explicit tip_mass. Zero blades (bare tower) and the usual three-bladed rotors are unaffected. --- src/pybmodes/io/windio.py | 15 +++++++++++++++ tests/test_windio_rna.py | 16 ++++++++++++++++ 2 files changed, 31 insertions(+) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 1f4662b..dd91731 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -1093,6 +1093,21 @@ def _nac_geom(key: str) -> float: f"assembly.number_of_blades must be a non-negative integer; got " f"{n_blades!r}." ) + if n_blades in (1, 2): + # The rotor tensor uses the axisymmetric transverse split + # ``I_diam = I_polar/2 + …``, which is exact only for three or more + # evenly spaced blades (a symmetric mass ring has isotropic in-plane + # inertia). A one- or two-bladed rotor lies on a single diameter, so + # its transverse inertia is azimuth-dependent and no single static + # tower-top lump represents it. Reject rather than emit a corrupted + # fore-aft / side-side inertia (issue #130 review). + raise ValueError( + f"assembly.number_of_blades = {n_blades}: the auto-RNA rotor " + f"inertia assembly assumes an azimuthally symmetric rotor (three " + f"or more evenly spaced blades). A one- or two-bladed rotor has an " + f"azimuth-dependent transverse inertia that a single rigid tower-" + f"top lump cannot represent; pass an explicit tip_mass instead." + ) orientation = str(assembly.get("rotor_orientation", "upwind")).strip().lower() if orientation not in {"upwind", "downwind"}: raise ValueError( diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index 90b4f30..7c2bbbb 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -200,6 +200,22 @@ def _mk(cone: float): assert coned.mass == pytest.approx(flat.mass) +@pytest.mark.parametrize("n_bl", [1, 2]) +def test_rna_rejects_one_or_two_bladed_rotor( + tmp_path: pathlib.Path, n_bl: int +) -> None: + """A one- or two-bladed rotor is not azimuthally symmetric, so the + ``I_polar/2`` transverse split is invalid and the auto-RNA rejects it + with a message pointing at an explicit tip_mass (Codex review on #130).""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + o = _base_ontology() + o["assembly"]["number_of_blades"] = n_bl + with pytest.raises(ValueError, match="azimuthally symmetric"): + read_windio_rna(_write(o, tmp_path, f"nb_{n_bl}.yaml")) + + def test_rna_rotor_inertia_includes_sweep(tmp_path: pathlib.Path) -> None: """Sweep (reference_axis.y) sets an in-plane tangential distance from the shaft axis, so it adds ``y²`` to the rotor polar lever. A constant From 6690b5d93f4d2d3c052fe7274d3b37dba50de1f8 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sun, 5 Jul 2026 23:50:48 +0900 Subject: [PATCH 08/17] harden: reject non-physical (degree-valued) cone and uptilt angles WindIO stores angles in radians (SI): the IEA references carry cone_angle 0.0698 (4 deg), uptilt 0.1047 (6 deg) and root twist 0.2094 (12 deg), and the auto-RNA reads them directly, which the IEA-22 deck cross-check confirms to <0.5 % on CM. A hand-authored file that leaves an angle in degrees (cone_angle: 4 meaning 4 deg) would otherwise be evaluated at 4 rad = 229 deg and silently corrupt the coned lever and axial offset. Guard cone_angle and uptilt against a physical +/-45 deg bound so that mistake fails loudly with a message naming the radians convention. --- src/pybmodes/io/windio.py | 25 +++++++++++++++++++++++-- tests/test_windio_rna.py | 24 ++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 2 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index dd91731..29206b8 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -732,6 +732,27 @@ def _finite_float(value: Any, what: str) -> float: return f +def _finite_angle(value: Any, what: str, max_abs: float = float(np.pi / 4)) -> float: + """Coerce ``value`` to a finite angle in radians within a physical bound. + + WindIO stores angles in radians (SI). A rotor precone or shaft tilt beyond + ~45 deg is non-physical for a horizontal-axis turbine and almost always + means degrees were supplied by mistake (e.g. ``cone_angle: 4`` read as + 4 rad = 229 deg instead of 4 deg), which would silently corrupt the coned + lever and axial offset. Reject it with a clear message rather than + evaluate the trig at the wrong angle. + """ + ang = _finite_float(value, what) + if abs(ang) > max_abs: + raise ValueError( + f"{what} = {ang:.4f} rad ({np.degrees(ang):.1f} deg) is outside the " + f"physical range +/-{np.degrees(max_abs):.0f} deg. WindIO angles are " + f"in radians; a value near a small integer usually means degrees " + f"were supplied." + ) + return ang + + def _positive_mass(value: Any, what: str) -> float: """Coerce ``value`` to a positive, finite mass (kg), else raise.""" m = _finite_float(value, what) @@ -1055,7 +1076,7 @@ def _nac_geom(key: str) -> float: return _finite_float(val, f"nacelle.drivetrain.{key}") overhang = _nac_geom("overhang") - uptilt = _nac_geom("uptilt") + uptilt = _finite_angle(_nac_geom("uptilt"), "nacelle.drivetrain.uptilt") dist_tt_hub = _nac_geom("distance_tt_hub") # --- Hub: components..elastic_properties_mb --- @@ -1082,7 +1103,7 @@ def _nac_geom(key: str) -> float: ) if hub_r < 0.0: raise ValueError(f"hub diameter must be non-negative; got {hub_diam!r}.") - cone_ang = _finite_float(hub.get("cone_angle", 0.0), "hub cone_angle") + cone_ang = _finite_angle(hub.get("cone_angle", 0.0), "hub cone_angle") cone_cos = float(np.cos(cone_ang)) cone_sin = float(np.sin(cone_ang)) diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index 7c2bbbb..f980b2d 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -200,6 +200,30 @@ def _mk(cone: float): assert coned.mass == pytest.approx(flat.mass) +def test_rna_rejects_degree_valued_cone_angle(tmp_path: pathlib.Path) -> None: + """WindIO angles are radians; a degrees-valued cone_angle (e.g. 4 read as + 4 rad = 229 deg) is non-physical and rejected rather than silently + evaluated at the wrong angle (Codex review on #130).""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + o = _base_ontology() + o["components"]["hub"]["cone_angle"] = 4.0 # 4 deg mistakenly left un-converted + with pytest.raises(ValueError, match="radians"): + read_windio_rna(_write(o, tmp_path, "cone_deg.yaml")) + + +def test_rna_rejects_degree_valued_uptilt(tmp_path: pathlib.Path) -> None: + """A degrees-valued uptilt is likewise rejected as out of physical range.""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + o = _base_ontology() + o["components"]["nacelle"]["drivetrain"]["uptilt"] = 6.0 # 6 deg un-converted + with pytest.raises(ValueError, match="radians"): + read_windio_rna(_write(o, tmp_path, "uptilt_deg.yaml")) + + @pytest.mark.parametrize("n_bl", [1, 2]) def test_rna_rejects_one_or_two_bladed_rotor( tmp_path: pathlib.Path, n_bl: int From 2900ff121516aa87e940115b9ef900757e216b25 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 6 Jul 2026 00:05:20 +0900 Subject: [PATCH 09/17] fix: disambiguate WindIO cone / uptilt rad vs deg by magnitude The WindIO v2 turbine schema annotates hub.cone_angle and drivetrain uptilt as degrees, but the IEA reference ontologies store radians (IEA-22 cone_angle 0.0698 = 4 deg, uptilt 0.1047 = 6 deg, matched to <0.5 % on CM against the ElastoDyn deck). A pure-radians reader rejects schema- conforming degree files; a pure-degrees reader mis-reads the reference files as ~0 precone. Physical precone and shaft tilt are small (well under ~20 deg), so the two encodings never overlap: interpret a stored magnitude at or below 0.5 as radians and a larger one as degrees, and reject a value non-physical under either (> 90). Both a 4 and a 0.0698 cone now assemble to the same 4 deg RNA. --- src/pybmodes/io/windio.py | 45 ++++++++++++++++++++++----------------- tests/test_windio_rna.py | 37 +++++++++++++++++++++----------- 2 files changed, 50 insertions(+), 32 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 29206b8..a61f14e 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -732,25 +732,32 @@ def _finite_float(value: Any, what: str) -> float: return f -def _finite_angle(value: Any, what: str, max_abs: float = float(np.pi / 4)) -> float: - """Coerce ``value`` to a finite angle in radians within a physical bound. - - WindIO stores angles in radians (SI). A rotor precone or shaft tilt beyond - ~45 deg is non-physical for a horizontal-axis turbine and almost always - means degrees were supplied by mistake (e.g. ``cone_angle: 4`` read as - 4 rad = 229 deg instead of 4 deg), which would silently corrupt the coned - lever and axial offset. Reject it with a clear message rather than - evaluate the trig at the wrong angle. +def _windio_angle(value: Any, what: str) -> float: + """Interpret a WindIO precone / shaft-tilt angle as radians, tolerant of + the rad/deg ambiguity across WindIO versions. + + The WindIO v2 turbine schema annotates ``hub.cone_angle`` and drivetrain + ``uptilt`` as degrees, but the IEA reference ontologies store radians + (IEA-22 ``cone_angle`` 0.0698 = 4 deg, ``uptilt`` 0.1047 = 6 deg). + Physical rotor precone and shaft tilt are small (well under ~20 deg), so + the two encodings never overlap: a stored magnitude at or below 0.5 is + radians (0.5 rad = 29 deg, already larger than any real precone / tilt), + a larger one is degrees. Convert to radians on that basis, and reject a + value that is non-physical under either reading (> 90). This keeps the + validated radians reference files working while also accepting a + schema-conforming degrees file (issue #130 review). """ ang = _finite_float(value, what) - if abs(ang) > max_abs: - raise ValueError( - f"{what} = {ang:.4f} rad ({np.degrees(ang):.1f} deg) is outside the " - f"physical range +/-{np.degrees(max_abs):.0f} deg. WindIO angles are " - f"in radians; a value near a small integer usually means degrees " - f"were supplied." - ) - return ang + mag = abs(ang) + if mag <= 0.5: + return ang + if mag <= 90.0: + return float(np.radians(ang)) + raise ValueError( + f"{what} = {ang!r} is non-physical as either radians " + f"({np.degrees(ang):.0f} deg) or degrees; expected a rotor precone / " + f"shaft tilt well under 90 deg." + ) def _positive_mass(value: Any, what: str) -> float: @@ -1076,7 +1083,7 @@ def _nac_geom(key: str) -> float: return _finite_float(val, f"nacelle.drivetrain.{key}") overhang = _nac_geom("overhang") - uptilt = _finite_angle(_nac_geom("uptilt"), "nacelle.drivetrain.uptilt") + uptilt = _windio_angle(_nac_geom("uptilt"), "nacelle.drivetrain.uptilt") dist_tt_hub = _nac_geom("distance_tt_hub") # --- Hub: components..elastic_properties_mb --- @@ -1103,7 +1110,7 @@ def _nac_geom(key: str) -> float: ) if hub_r < 0.0: raise ValueError(f"hub diameter must be non-negative; got {hub_diam!r}.") - cone_ang = _finite_angle(hub.get("cone_angle", 0.0), "hub cone_angle") + cone_ang = _windio_angle(hub.get("cone_angle", 0.0), "hub cone_angle") cone_cos = float(np.cos(cone_ang)) cone_sin = float(np.sin(cone_ang)) diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index f980b2d..d3ec933 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -200,28 +200,39 @@ def _mk(cone: float): assert coned.mass == pytest.approx(flat.mass) -def test_rna_rejects_degree_valued_cone_angle(tmp_path: pathlib.Path) -> None: - """WindIO angles are radians; a degrees-valued cone_angle (e.g. 4 read as - 4 rad = 229 deg) is non-physical and rejected rather than silently - evaluated at the wrong angle (Codex review on #130).""" +def test_rna_cone_and_uptilt_degrees_match_radians(tmp_path: pathlib.Path) -> None: + """WindIO's rad/deg ambiguity: the v2 schema annotates cone_angle / uptilt + as degrees while the IEA reference files store radians. The reader + disambiguates by magnitude, so a degrees encoding (4, 6) and the + equivalent radians encoding (0.0698, 0.1047) assemble to the same RNA + (Codex review on #130).""" pytest.importorskip("yaml") + import numpy as np + from pybmodes.io.windio import read_windio_rna - o = _base_ontology() - o["components"]["hub"]["cone_angle"] = 4.0 # 4 deg mistakenly left un-converted - with pytest.raises(ValueError, match="radians"): - read_windio_rna(_write(o, tmp_path, "cone_deg.yaml")) + def _mk(cone: float, uptilt: float): + o = _base_ontology() + o["components"]["hub"]["cone_angle"] = cone + o["components"]["nacelle"]["drivetrain"]["uptilt"] = uptilt + return read_windio_rna(_write(o, tmp_path, f"a_{cone}_{uptilt}.yaml")) + + rad = _mk(float(np.radians(4.0)), float(np.radians(6.0))) # 0.0698, 0.1047 + deg = _mk(4.0, 6.0) # degrees encoding of the same angles + for attr in ("mass", "cm_axial", "ixx", "iyy", "izz", "izx"): + assert getattr(deg, attr) == pytest.approx(getattr(rad, attr)) -def test_rna_rejects_degree_valued_uptilt(tmp_path: pathlib.Path) -> None: - """A degrees-valued uptilt is likewise rejected as out of physical range.""" +def test_rna_rejects_nonphysical_angle(tmp_path: pathlib.Path) -> None: + """A value non-physical as both radians and degrees (e.g. 200) is rejected + rather than silently evaluated at a meaningless angle.""" pytest.importorskip("yaml") from pybmodes.io.windio import read_windio_rna o = _base_ontology() - o["components"]["nacelle"]["drivetrain"]["uptilt"] = 6.0 # 6 deg un-converted - with pytest.raises(ValueError, match="radians"): - read_windio_rna(_write(o, tmp_path, "uptilt_deg.yaml")) + o["components"]["hub"]["cone_angle"] = 200.0 + with pytest.raises(ValueError, match="non-physical"): + read_windio_rna(_write(o, tmp_path, "cone_bad.yaml")) @pytest.mark.parametrize("n_bl", [1, 2]) From fb42ee500d9d5d528107c2d892b6684d4ef50c02 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 6 Jul 2026 00:13:43 +0900 Subject: [PATCH 10/17] fix: resolve WindIO rad/deg once per file from the larger rotor angle The per-value magnitude rule misread a small schema-conforming degrees angle: a WindIO v2 file with cone_angle 0.25 (deg) fell into the <=0.5 radians branch and became 0.25 rad = 14.3 deg. A file uses one angle convention throughout, so decide it once for the (cone_angle, uptilt) pair from the larger magnitude. A degrees file reliably carries a several-degree shaft tilt whose magnitude dwarfs any radians precone/tilt (which stay under ~0.3 rad), so the larger angle disambiguates reliably even when the cone is a fraction of a degree. Larger magnitude over 0.5 => degrees, else radians; over 90 rejected either way. --- src/pybmodes/io/windio.py | 50 ++++++++++++++++++++------------------- tests/test_windio_rna.py | 22 +++++++++++++++++ 2 files changed, 48 insertions(+), 24 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index a61f14e..3f6bc7c 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -732,32 +732,34 @@ def _finite_float(value: Any, what: str) -> float: return f -def _windio_angle(value: Any, what: str) -> float: - """Interpret a WindIO precone / shaft-tilt angle as radians, tolerant of - the rad/deg ambiguity across WindIO versions. +def _windio_rotor_angles(cone_value: Any, uptilt_value: Any) -> tuple[float, float]: + """Convert the WindIO precone and shaft tilt to radians, resolving the + rad/deg ambiguity once for the pair. The WindIO v2 turbine schema annotates ``hub.cone_angle`` and drivetrain ``uptilt`` as degrees, but the IEA reference ontologies store radians - (IEA-22 ``cone_angle`` 0.0698 = 4 deg, ``uptilt`` 0.1047 = 6 deg). - Physical rotor precone and shaft tilt are small (well under ~20 deg), so - the two encodings never overlap: a stored magnitude at or below 0.5 is - radians (0.5 rad = 29 deg, already larger than any real precone / tilt), - a larger one is degrees. Convert to radians on that basis, and reject a - value that is non-physical under either reading (> 90). This keeps the - validated radians reference files working while also accepting a - schema-conforming degrees file (issue #130 review). + (IEA-22 ``cone_angle`` 0.0698 = 4 deg, ``uptilt`` 0.1047 = 6 deg). A file + uses one convention throughout, so decide it once from the *larger* of the + two magnitudes rather than per value: a degrees file reliably carries a + several-degree shaft tilt (magnitude well above any radians precone / + tilt, which stay under ~0.3 rad), whereas judging each value alone would + misread a small schema-conforming degrees angle such as 0.25 deg as + 0.25 rad. If the larger magnitude exceeds 0.5 the pair is degrees; + otherwise radians. Anything above 90 is non-physical either way and is + rejected (issue #130 review). """ - ang = _finite_float(value, what) - mag = abs(ang) - if mag <= 0.5: - return ang - if mag <= 90.0: - return float(np.radians(ang)) - raise ValueError( - f"{what} = {ang!r} is non-physical as either radians " - f"({np.degrees(ang):.0f} deg) or degrees; expected a rotor precone / " - f"shaft tilt well under 90 deg." - ) + cone = _finite_float(cone_value, "hub cone_angle") + uptilt = _finite_float(uptilt_value, "nacelle.drivetrain.uptilt") + biggest = max(abs(cone), abs(uptilt)) + if biggest > 90.0: + raise ValueError( + f"rotor precone / shaft tilt (cone_angle={cone!r}, uptilt=" + f"{uptilt!r}) is non-physical as either radians or degrees; " + f"expected values well under 90 deg." + ) + if biggest > 0.5: # degrees: a real shaft tilt dwarfs any radians angle + return float(np.radians(cone)), float(np.radians(uptilt)) + return cone, uptilt def _positive_mass(value: Any, what: str) -> float: @@ -1083,7 +1085,7 @@ def _nac_geom(key: str) -> float: return _finite_float(val, f"nacelle.drivetrain.{key}") overhang = _nac_geom("overhang") - uptilt = _windio_angle(_nac_geom("uptilt"), "nacelle.drivetrain.uptilt") + uptilt_raw = _nac_geom("uptilt") # rad/deg resolved with cone_angle below dist_tt_hub = _nac_geom("distance_tt_hub") # --- Hub: components..elastic_properties_mb --- @@ -1110,7 +1112,7 @@ def _nac_geom(key: str) -> float: ) if hub_r < 0.0: raise ValueError(f"hub diameter must be non-negative; got {hub_diam!r}.") - cone_ang = _windio_angle(hub.get("cone_angle", 0.0), "hub cone_angle") + cone_ang, uptilt = _windio_rotor_angles(hub.get("cone_angle", 0.0), uptilt_raw) cone_cos = float(np.cos(cone_ang)) cone_sin = float(np.sin(cone_ang)) diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index d3ec933..8623909 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -223,6 +223,28 @@ def _mk(cone: float, uptilt: float): assert getattr(deg, attr) == pytest.approx(getattr(rad, attr)) +def test_rna_small_degree_cone_keyed_off_uptilt(tmp_path: pathlib.Path) -> None: + """A schema-conforming degrees file with a small cone (0.25 deg) and a + normal 6 deg uptilt: the unit decision keys off the larger (uptilt) + magnitude, so the small cone is read as 0.25 deg, not 0.25 rad = 14.3 deg + (Codex review on #130).""" + pytest.importorskip("yaml") + import numpy as np + + from pybmodes.io.windio import read_windio_rna + + def _mk(cone: float, uptilt: float): + o = _base_ontology() + o["components"]["hub"]["cone_angle"] = cone + o["components"]["nacelle"]["drivetrain"]["uptilt"] = uptilt + return read_windio_rna(_write(o, tmp_path, f"s_{cone}_{uptilt}.yaml")) + + deg = _mk(0.25, 6.0) # degrees, small cone + rad = _mk(float(np.radians(0.25)), float(np.radians(6.0))) # radians equivalent + for attr in ("cm_axial", "ixx", "iyy", "izz"): + assert getattr(deg, attr) == pytest.approx(getattr(rad, attr)) + + def test_rna_rejects_nonphysical_angle(tmp_path: pathlib.Path) -> None: """A value non-physical as both radians and degrees (e.g. 200) is rejected rather than silently evaluated at a meaningless angle.""" From e3d6d81759b42b3f8e632e603fe8ac24e72547ee Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 6 Jul 2026 00:24:33 +0900 Subject: [PATCH 11/17] feat: add explicit angle_units override for WindIO rotor angles The magnitude heuristic cannot resolve a degrees file whose cone and tilt are both sub-half-degree (a zero-tilt turbine with a fractional-degree precone reads as radians). Rejecting small values is not an option either, since a legitimate radians file routinely has a small precone with zero tilt. Add an explicit escape hatch instead: read_windio_rna gains angle_units ('auto' | 'rad' | 'deg'), forwarded from Tower.from_windio and from_windio_with_monopile as rna_angle_units. 'auto' (default) keeps the pair heuristic for the common case; 'rad' / 'deg' take the file at its word for the rare ambiguous one. Unknown values are rejected. --- src/pybmodes/io/windio.py | 58 +++++++++++++++++++++++++----------- src/pybmodes/models/tower.py | 14 +++++++-- tests/test_windio_rna.py | 36 ++++++++++++++++++++++ 3 files changed, 88 insertions(+), 20 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 3f6bc7c..b4347dd 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -732,32 +732,45 @@ def _finite_float(value: Any, what: str) -> float: return f -def _windio_rotor_angles(cone_value: Any, uptilt_value: Any) -> tuple[float, float]: - """Convert the WindIO precone and shaft tilt to radians, resolving the - rad/deg ambiguity once for the pair. +def _windio_rotor_angles( + cone_value: Any, uptilt_value: Any, units: str = "auto" +) -> tuple[float, float]: + """Convert the WindIO precone and shaft tilt to radians. The WindIO v2 turbine schema annotates ``hub.cone_angle`` and drivetrain ``uptilt`` as degrees, but the IEA reference ontologies store radians - (IEA-22 ``cone_angle`` 0.0698 = 4 deg, ``uptilt`` 0.1047 = 6 deg). A file - uses one convention throughout, so decide it once from the *larger* of the - two magnitudes rather than per value: a degrees file reliably carries a - several-degree shaft tilt (magnitude well above any radians precone / - tilt, which stay under ~0.3 rad), whereas judging each value alone would - misread a small schema-conforming degrees angle such as 0.25 deg as - 0.25 rad. If the larger magnitude exceeds 0.5 the pair is degrees; - otherwise radians. Anything above 90 is non-physical either way and is - rejected (issue #130 review). + (IEA-22 ``cone_angle`` 0.0698 = 4 deg, ``uptilt`` 0.1047 = 6 deg). + ``units`` selects how to reconcile that: + + - ``"rad"`` / ``"deg"`` — take the file at its word (no guessing). + - ``"auto"`` (default) — a file uses one convention throughout, so decide + it once from the *larger* of the two magnitudes: a degrees file reliably + carries a several-degree shaft tilt whose magnitude dwarfs any radians + precone / tilt (which stay under ~0.3 rad), so the larger angle + disambiguates the pair even when the cone is a fraction of a degree. + Larger magnitude over 0.5 means degrees, otherwise radians. The + residual blind spot is a degrees file whose angles are *both* below + ~0.5 deg (e.g. a zero-tilt turbine with a sub-degree precone); it reads + as radians. Pass ``units="deg"`` for such a file (issue #130 review). + + Anything above 90 (in the resolved unit) is non-physical and rejected. """ cone = _finite_float(cone_value, "hub cone_angle") uptilt = _finite_float(uptilt_value, "nacelle.drivetrain.uptilt") - biggest = max(abs(cone), abs(uptilt)) - if biggest > 90.0: + if units not in ("auto", "rad", "deg"): + raise ValueError( + f"angle_units must be 'auto', 'rad' or 'deg'; got {units!r}." + ) + if units == "auto": + biggest = max(abs(cone), abs(uptilt)) + units = "deg" if biggest > 0.5 else "rad" + if max(abs(cone), abs(uptilt)) > 90.0: raise ValueError( f"rotor precone / shaft tilt (cone_angle={cone!r}, uptilt=" - f"{uptilt!r}) is non-physical as either radians or degrees; " - f"expected values well under 90 deg." + f"{uptilt!r}) is non-physical in {units}; expected values well " + f"under 90 deg." ) - if biggest > 0.5: # degrees: a real shaft tilt dwarfs any radians angle + if units == "deg": return float(np.radians(cone)), float(np.radians(uptilt)) return cone, uptilt @@ -980,6 +993,7 @@ def read_windio_rna( component_hub: str = "hub", component_nacelle: str = "nacelle", component_blade: str = "blade", + angle_units: str = "auto", ) -> TipMassProps: """Lump the WindIO rotor-nacelle assembly into a tower-top ``TipMassProps``. @@ -1024,6 +1038,12 @@ def read_windio_rna( goes beyond the ElastoDyn deck path's bare point-mass lumping, which the WindIO ontology's per-station blade mass makes possible. + ``angle_units`` selects how ``cone_angle`` / ``uptilt`` are read. The + WindIO v2 schema annotates them as degrees, but the IEA reference files + store radians; ``"auto"`` (default) disambiguates by magnitude (see + :func:`_windio_rotor_angles`), while ``"rad"`` / ``"deg"`` take the file + at its word for the rare all-sub-degree case ``"auto"`` cannot resolve. + Requires the optional ``[windio]`` extra (PyYAML). """ from pybmodes.io.bmi import TipMassProps @@ -1112,7 +1132,9 @@ def _nac_geom(key: str) -> float: ) if hub_r < 0.0: raise ValueError(f"hub diameter must be non-negative; got {hub_diam!r}.") - cone_ang, uptilt = _windio_rotor_angles(hub.get("cone_angle", 0.0), uptilt_raw) + cone_ang, uptilt = _windio_rotor_angles( + hub.get("cone_angle", 0.0), uptilt_raw, angle_units + ) cone_cos = float(np.cos(cone_ang)) cone_sin = float(np.sin(cone_ang)) diff --git a/src/pybmodes/models/tower.py b/src/pybmodes/models/tower.py index 7eee652..c417923 100644 --- a/src/pybmodes/models/tower.py +++ b/src/pybmodes/models/tower.py @@ -325,6 +325,7 @@ def from_windio( tip_mass: TipMassProps | float | None = None, n_nodes: int | None = None, lumped_rna_cal: bool = False, + rna_angle_units: str = "auto", ) -> Tower: """Build a tower (or monopile) model directly from a **WindIO** ontology ``.yaml`` (issue #35). @@ -372,6 +373,11 @@ def from_windio( expressed at the tower top in the clamped-base convention, which a free-base / soil-flexible base interprets differently. Default ``False``. + rna_angle_units : how the auto-RNA reads ``cone_angle`` / ``uptilt`` + when ``lumped_rna_cal=True``. ``"auto"`` (default) disambiguates + the WindIO rad/deg ambiguity by magnitude; pass ``"rad"`` or + ``"deg"`` to take the file at its word. Ignored unless + ``lumped_rna_cal`` is set. Notes ----- @@ -411,7 +417,7 @@ def from_windio( ) from pybmodes.io.windio import read_windio_rna - tip_mass = read_windio_rna(yaml_path) + tip_mass = read_windio_rna(yaml_path, angle_units=rna_angle_units) g = read_windio_tubular( yaml_path, component=component, thickness_interp=thickness_interp, @@ -440,6 +446,7 @@ def from_windio_with_monopile( n_nodes: int | None = None, water_depth: float | None = None, lumped_rna_cal: bool = False, + rna_angle_units: str = "auto", ) -> Tower: """Build a combined **monopile + tower** fixed-bottom cantilever from a WindIO ontology ``.yaml`` (issue #92). @@ -484,6 +491,9 @@ def from_windio_with_monopile( blocks via :func:`pybmodes.io.windio.read_windio_rna` and use it as ``tip_mass``. Requires an IEA-22-class ontology; mutually exclusive with ``tip_mass``. Default ``False``. + rna_angle_units : how the auto-RNA reads ``cone_angle`` / ``uptilt`` + when ``lumped_rna_cal=True`` (``"auto"`` / ``"rad"`` / ``"deg"``); + see :meth:`from_windio`. Ignored unless ``lumped_rna_cal`` is set. Notes ----- @@ -509,7 +519,7 @@ def from_windio_with_monopile( ) from pybmodes.io.windio import read_windio_rna - tip_mass = read_windio_rna(yaml_path) + tip_mass = read_windio_rna(yaml_path, angle_units=rna_angle_units) mt = read_windio_monopile_tower( yaml_path, diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index 8623909..4346a48 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -245,6 +245,42 @@ def _mk(cone: float, uptilt: float): assert getattr(deg, attr) == pytest.approx(getattr(rad, attr)) +def test_rna_explicit_deg_units_override(tmp_path: pathlib.Path) -> None: + """The all-sub-degree case auto cannot resolve: a zero-tilt degrees file + with a 0.3 deg cone. angle_units='deg' takes the file at its word, so it + reads as 0.3 deg, not 0.3 rad (Codex review on #130).""" + pytest.importorskip("yaml") + import numpy as np + + from pybmodes.io.windio import read_windio_rna + + o = _base_ontology() + o["components"]["hub"]["cone_angle"] = 0.3 + o["components"]["nacelle"]["drivetrain"]["uptilt"] = 0.0 + p = _write(o, tmp_path, "small_deg.yaml") + + as_deg = read_windio_rna(p, angle_units="deg") + # equivalent radians file resolved by auto + o2 = _base_ontology() + o2["components"]["hub"]["cone_angle"] = float(np.radians(0.3)) + o2["components"]["nacelle"]["drivetrain"]["uptilt"] = 0.0 + as_rad = read_windio_rna(_write(o2, tmp_path, "small_rad.yaml")) + for attr in ("cm_axial", "ixx", "iyy", "izz"): + assert getattr(as_deg, attr) == pytest.approx(getattr(as_rad, attr)) + # and auto (default) would misread the 0.3 as radians -> different result + as_auto = read_windio_rna(p) + assert as_auto.izz != pytest.approx(as_deg.izz) + + +def test_rna_rejects_bad_angle_units(tmp_path: pathlib.Path) -> None: + """An unknown angle_units value is rejected.""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + with pytest.raises(ValueError, match="angle_units"): + read_windio_rna(_write(_base_ontology(), tmp_path, "u.yaml"), angle_units="grad") + + def test_rna_rejects_nonphysical_angle(tmp_path: pathlib.Path) -> None: """A value non-physical as both radians and degrees (e.g. 200) is rejected rather than silently evaluated at a meaningless angle.""" From 6871e4f686bbfd1b040e06ee125416217a9e0db8 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 6 Jul 2026 00:31:45 +0900 Subject: [PATCH 12/17] fix: bound the resolved radian rotor angle, not the raw input The 90 deg physical guard compared the raw cone / uptilt value to 90, which is only meaningful when the input is in degrees. Under angle_units=rad a 2.0 rad (~114 deg) angle passed the 2.0 < 90 check and assembled a plausible RNA for a non-physical geometry. Resolve to radians first, then bound the magnitude at pi/2, so the 90 deg limit holds whether the file was rad or deg. --- src/pybmodes/io/windio.py | 16 ++++++++++------ tests/test_windio_rna.py | 14 ++++++++++++++ 2 files changed, 24 insertions(+), 6 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index b4347dd..dcd7797 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -764,15 +764,19 @@ def _windio_rotor_angles( if units == "auto": biggest = max(abs(cone), abs(uptilt)) units = "deg" if biggest > 0.5 else "rad" - if max(abs(cone), abs(uptilt)) > 90.0: + if units == "deg": + cone_rad, uptilt_rad = float(np.radians(cone)), float(np.radians(uptilt)) + else: + cone_rad, uptilt_rad = cone, uptilt + # Bound the resolved radian value (not the raw input, whose scale depends + # on the unit), so a 90 deg limit holds whether the file was rad or deg. + if max(abs(cone_rad), abs(uptilt_rad)) > np.pi / 2: raise ValueError( f"rotor precone / shaft tilt (cone_angle={cone!r}, uptilt=" - f"{uptilt!r}) is non-physical in {units}; expected values well " - f"under 90 deg." + f"{uptilt!r}, units={units}) resolves to more than 90 deg, which " + f"is non-physical." ) - if units == "deg": - return float(np.radians(cone)), float(np.radians(uptilt)) - return cone, uptilt + return cone_rad, uptilt_rad def _positive_mass(value: Any, what: str) -> float: diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index 4346a48..a9bc680 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -293,6 +293,20 @@ def test_rna_rejects_nonphysical_angle(tmp_path: pathlib.Path) -> None: read_windio_rna(_write(o, tmp_path, "cone_bad.yaml")) +def test_rna_rejects_over_90deg_explicit_radians(tmp_path: pathlib.Path) -> None: + """The 90 deg physical bound is applied to the resolved radian value, so a + 2.0 rad (~114 deg) cone under angle_units='rad' is rejected rather than + slipping past a raw < 90 comparison (Codex review on #130).""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + o = _base_ontology() + o["components"]["hub"]["cone_angle"] = 2.0 # rad (~114 deg), non-physical + p = _write(o, tmp_path, "cone_2rad.yaml") + with pytest.raises(ValueError, match="non-physical"): + read_windio_rna(p, angle_units="rad") + + @pytest.mark.parametrize("n_bl", [1, 2]) def test_rna_rejects_one_or_two_bladed_rotor( tmp_path: pathlib.Path, n_bl: int From e38287b734e678a39ba0ee616d9aac0166bac3e9 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 6 Jul 2026 00:42:06 +0900 Subject: [PATCH 13/17] fix: project the hub radius through precone in the rotor lever The polar lever left the hub radius unconed in the rotor plane and dropped it from the axial offset. WISDEM and ElastoDyn place the blade root at the hub radius along the coned pitch axis, so the whole (hub_r + span) length projects through the cone (WISDEM rotorR = Rtip*cos(precone) with Rtip = Rhub + L). Fold hub_r into the along-axis distance before projecting: radial = (hub_r + span)*cos(cone) - x*sin(cone), axial = (hub_r + span)* sin(cone) + x*cos(cone). For a coned rotor with a real hub radius this adds the Rhub*sin(cone) shaft-axis offset that was missing. Unconed rotors (cone = 0) are unchanged. --- src/pybmodes/io/windio.py | 21 ++++++++++++++------- tests/test_windio_rna.py | 21 +++++++++++++++++++++ 2 files changed, 35 insertions(+), 7 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index dcd7797..991bafa 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -891,10 +891,11 @@ def _blade_span_mass_inertia( Returns the single-blade mass, the polar second moment about the rotor axis ``∫ (dm/ds) · (r² + t²) ds`` (in-plane radial - ``r = hub_r + z·cone_cos − x·cone_sin`` and tangential ``t = y`` so - prebend / sweep are included), the axial second moment - ``∫ (dm/ds) · a(s)² ds`` (axial offset ``a = z·cone_sin + x·cone_cos``), - and the axial first moment ``∫ (dm/ds) · a(s) ds``. + ``r = (hub_r + z)·cone_cos − x·cone_sin`` and tangential ``t = y`` so the + coned hub radius, prebend and sweep are all included), the axial second + moment ``∫ (dm/ds) · a(s)² ds`` (axial offset + ``a = (hub_r + z)·cone_sin + x·cone_cos``), and the axial first moment + ``∫ (dm/ds) · a(s) ds``. The caller builds the rotor's polar inertia from the second radial moment, the coned transverse (diametral) inertia from the radial and axial second moments, and places the rotor body at its coned centre of @@ -980,9 +981,15 @@ def _blade_span_mass_inertia( # ``r² + t²`` (so sweep and prebend are not dropped, issue #130); the # axial offset feeds the coned transverse term. x_off, y_off, span = coords[0], coords[1], coords[2] - radial = hub_r + span * cone_cos - x_off * cone_sin + # Distance from the hub centre along the (coned) pitch axis: the blade + # root sits at the hub radius and each section a further `span` outboard, + # so the whole (hub_r + span) length projects through the cone, matching + # WISDEM/ElastoDyn (rotorR = Rtip*cos(precone) = (Rhub + L)*cos(cone)) + # rather than leaving the hub radius unconed in the rotor plane. + axis_dist = hub_r + span + radial = axis_dist * cone_cos - x_off * cone_sin tangential = y_off - axial = span * cone_sin + x_off * cone_cos + axial = axis_dist * cone_sin + x_off * cone_cos polar_second_moment = _second_moment(mpl, radial, s) + _second_moment( mpl, tangential, s ) @@ -1035,7 +1042,7 @@ def read_windio_rna( Blade inertia (issue #130): the rotor carries the diametral inertia ``N_bl · ∫ (dm/ds) · r² ds`` from the blade mass spread along the span - (``r = hub_radius + span·cos(cone)``), as ``diag([I_polar, I_polar/2, + (``r = (hub_radius + span)·cos(cone)``), as ``diag([I_polar, I_polar/2, I_polar/2])`` about the hub (the same perpendicular-axis form as the hub tensor). Only each blade's own *sectional* spin inertia is left out (its parallel-axis / span contribution is captured). This intentionally diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index a9bc680..a8e11cd 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -347,6 +347,27 @@ def _mk(y_off: float): assert swept.izz > straight.izz +def test_rna_hub_radius_projected_through_cone(tmp_path: pathlib.Path) -> None: + """The hub radius is coned like the span (WISDEM Rtip*cos(precone)), so a + coned rotor with a hub radius carries an extra Rhub*sin(cone) axial offset + that raises the tower-top vertical CM through the shaft tilt, versus the + same rotor with a zero hub radius (Codex review on #130).""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + def _mk(hub_d: float): + o = _base_ontology() + o["components"]["hub"]["cone_angle"] = 0.15 + o["components"]["hub"]["diameter"] = hub_d + o["components"]["nacelle"]["drivetrain"]["uptilt"] = 0.10 + return read_windio_rna(_write(o, tmp_path, f"hr_{hub_d}.yaml")) + + no_hub = _mk(0.0) # hub_r = 0 + with_hub = _mk(6.0) # hub_r = 3 + assert with_hub.mass == pytest.approx(no_hub.mass) + assert with_hub.cm_axial > no_hub.cm_axial + + def test_rna_coned_cm_shift_with_uptilt(tmp_path: pathlib.Path) -> None: """A coned rotor's CM sits off the hub apex along the (tilted) shaft, so with uptilt the tower-top vertical CM shifts by m_blades·offset·sin(uptilt) From 6e27d40cf7ce01a27fef9a9f13dc65a6cb034d5e Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 6 Jul 2026 00:49:28 +0900 Subject: [PATCH 14/17] docs: sync the rotor-lever comment with the coned hub-radius formula A stale inline comment still described the polar lever as r = hub_r + span*cos(cone), from before the hub radius was coned. Update it to r = (hub_r + span)*cos(cone) so the comment matches the implementation and the surrounding docstrings. --- src/pybmodes/io/windio.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 991bafa..1ca78b6 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -1206,9 +1206,10 @@ def _nac_geom(key: str) -> float: # Rotor inertia from the spanwise blade mass (issue #130). In the shaft # frame the rotor tensor is diag([I_polar, I_diam, I_diam]): the polar - # term ``I_polar = N_bl·∫dm·r²`` (radial ``r = hub_r + span·cos(cone)``), - # and the transverse (diametral) term. For a coned rotor the blade mass - # sits off the hub plane by ``a = span·sin(cone)``, so the rotor CM shifts + # term ``I_polar = N_bl·∫dm·r²`` (radial ``r = (hub_r + span)·cos(cone)``, + # the coned pitch-axis distance from the hub centre), and the transverse + # (diametral) term. For a coned rotor the blade mass sits off the hub + # plane by ``a = (hub_r + span)·sin(cone)``, so the rotor CM shifts # along the shaft axis by ``offset = ∫dm·a / m_blades`` and the tensor is # formed about that CM: ``I_diam = I_polar/2 + N_bl·∫dm·a² − m_blades· # offset²`` (the axial second moment reduced to the CM by the parallel- From e0b1fcee5d63be320226b677cece9b69cc09f81f Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 6 Jul 2026 01:02:39 +0900 Subject: [PATCH 15/17] docs: correct the auto-RNA changelog note for the coned rotor The note described the pre-review lever (r = hub_radius + span*cos(cone)) and claimed the centre of mass was unchanged. The rotor is now assembled about its true (coned) centre of mass with the coned pitch-axis lever (hub_radius + span)*cos(cone) including prebend and sweep, and cm_axial shifts by the precone term for a non-zero cone_angle. Update the note to match: the tensor is diag([I_polar, I_diam, I_diam]) about the rotor CM, and both inertia and (for coned rotors) the CM change relative to 1.16.0; only total mass is unchanged. --- CHANGELOG.md | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9f998d8..a2810ed 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,14 +19,19 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). inertia (on IEA-22 it is ~2.5x the value the point-mass lump produced), so tower fore-aft / side-side frequencies from a rigid-RNA lump came out too high versus a rigid-rotor reference. The rotor is now assembled as a - rigid body, `diag([I_polar, I_polar/2, I_polar/2])` about the hub with - `I_polar = N_bl · ∫ (dm/ds) · r² ds` (`r = hub_radius + span·cos(cone)`), - using the hub diameter and cone angle when present. Only each blade's own - sectional spin inertia is still excluded. This changes the auto-RNA - inertia (and hence coupled tower frequencies) relative to 1.16.0; mass - and centre of mass are unchanged. This intentionally goes beyond the - ElastoDyn deck path's point-mass lumping, which the WindIO per-station - blade mass makes possible. + rigid body, `diag([I_polar, I_diam, I_diam])` about the rotor centre of + mass with `I_polar = N_bl · ∫ (dm/ds) · r² ds` and the in-plane lever + `r = (hub_radius + span)·cos(cone)` (the coned pitch-axis distance, with + prebend and sweep folded in), using the hub diameter and cone angle when + present. A coned rotor's mass sits off the hub plane, so the rotor is + placed at its true centre of mass before the parallel-axis shift; this + also moves the tower-top centre of mass slightly (by the precone term) + for a non-zero `cone_angle`. Only each blade's own sectional spin inertia + is still excluded. This changes the auto-RNA inertia (and hence coupled + tower frequencies), and the CM for coned rotors, relative to 1.16.0; the + total mass is unchanged. This intentionally goes beyond the ElastoDyn + deck path's point-mass lumping, which the WindIO per-station blade mass + makes possible. ## [1.16.0] — 2026-07-03 From bb3ddad7b0ddbc6066e697b7942925df7cc1dc6f Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 6 Jul 2026 05:53:22 +0900 Subject: [PATCH 16/17] fix: rotate the hub tensor by the correct WindIO hub-frame handedness The WindIO hub-aligned frame is the tower-top frame tilted about the lateral y-axis (x along the shaft, z up including tilt), so the shaft-tilt rotation carries sign_x on the sin (off-diagonal) terms, not the cos. The previous matrix put sign_x on the cos terms, which for an upwind rotor flipped the signs of a hub tensor Ixy / Iyz. A diagonal hub (every current reference, incl. IEA-22) and the symmetric rotor tensor are invariant to the choice, so results are unchanged there, but a full off-diagonal hub inertia is now oriented correctly. The outboard shaft axis for the coned-CM offset is built explicitly, since it is the shaft line (opposite sense to the hub-frame x basis for an upwind rotor). Also assert the >= 3-blade precondition at the axisymmetric split, which the blade-count guard upstream already enforces. --- src/pybmodes/io/windio.py | 42 +++++++++++++++++++++++---------------- tests/test_windio_rna.py | 19 ++++++++++++++++++ 2 files changed, 44 insertions(+), 17 deletions(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 1ca78b6..70a835d 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -1218,6 +1218,10 @@ def _nac_geom(key: str) -> float: # A flat rotor (``sin(cone) = 0``) leaves the offset and the axial term # zero, recovering the planar ``I_polar/2`` split about the apex. Lumping # the blades as a bare point mass at the apex would drop the whole tensor. + # n_blades is 0 or >= 3 here: 1 and 2 are rejected above, since the + # I_polar/2 transverse split holds only for an azimuthally symmetric rotor + # (>= 3 evenly spaced blades). n_blades == 0 gives a zero rotor tensor. + assert n_blades == 0 or n_blades >= 3 i_polar_rotor = n_blades * i_polar_each cm_axial_offset = (m_axial1_each / m_blade_each) if m_blade_each > 0.0 else 0.0 i_axial_cm = max( @@ -1235,25 +1239,29 @@ def _nac_geom(key: str) -> float: dtype=float, ) - # The hub and rotor inertia tensors are in the shaft frame; rotate them - # into the tower-top frame by the shaft tilt before they are summed as - # tower-top tensors (a tilted rotor otherwise loses the tensor izx - # product and misallocates the large polar term between ixx and izz). - # The shaft unit axis is ``[sign_x*cos(uptilt), 0, sin(uptilt)]``; R maps - # the shaft x-axis onto it and fixes the lateral y-axis. The nacelle - # tensor is already in the tower frame and is left alone. + # The hub and rotor inertia tensors are in the WindIO hub-aligned frame + # (x along the tilted shaft, y lateral, z up including tilt); rotate them + # into the tower-top frame before they are summed as tower-top tensors (a + # tilted rotor otherwise loses the tensor izx product and misallocates the + # large polar term between ixx and izz). The hub frame is the tower frame + # tilted about the lateral y-axis by the shaft tilt, so ``sign_x`` rides + # the sin (off-diagonal) terms, not the cos: for an upwind rotor the + # downwind-positive shaft tilts the opposite way. A diagonal hub or the + # symmetric rotor tensor is invariant to that choice, but a hub with + # off-diagonal Ixy / Iyz would be mis-signed with sign_x on the cos terms. + # The nacelle tensor is already in the tower frame and is left alone. c_t, s_t = float(np.cos(uptilt)), float(np.sin(uptilt)) - r_tilt = np.array( - [[sign_x * c_t, 0.0, -s_t], [0.0, 1.0, 0.0], [s_t, 0.0, sign_x * c_t]] + r_hub = np.array( + [[c_t, 0.0, -sign_x * s_t], [0.0, 1.0, 0.0], [sign_x * s_t, 0.0, c_t]] ) - i_hub = r_tilt @ i_hub @ r_tilt.T - i_blades = r_tilt @ i_blades @ r_tilt.T - - # The coned rotor CM sits ``cm_axial_offset`` beyond the apex along the - # shaft axis (r_tilt's first column, the outboard unit vector the overhang - # already runs along), so a precone displaces the blade mass further from - # the tower in the same sense for both rotor orientations. - shaft_axis = r_tilt[:, 0] + i_hub = r_hub @ i_hub @ r_hub.T + i_blades = r_hub @ i_blades @ r_hub.T + + # Outboard shaft direction (tower top toward the apex and beyond), along + # which the overhang runs and the coned blade CM is displaced. This is the + # shaft *line*, not the hub-frame x basis rotated above (they point in + # opposite senses for an upwind rotor), so build it explicitly. + shaft_axis = np.array([sign_x * c_t, 0.0, s_t]) rotor_cm = apex + cm_axial_offset * shaft_axis bodies = [ diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index a8e11cd..f242c8c 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -347,6 +347,25 @@ def _mk(y_off: float): assert swept.izz > straight.izz +def test_rna_upwind_hub_offdiagonal_inertia_sign(tmp_path: pathlib.Path) -> None: + """The WindIO hub frame has z up, so at zero tilt the hub tensor rotates by + the identity and its off-diagonal Iyz passes through unchanged, even for an + upwind rotor. A sign_x on the cos terms would flip it (Codex review on + #130). The hub sits on the tower axis (y = 0), so no other body adds iyz.""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + o = _base_ontology() + assert o["assembly"]["rotor_orientation"] == "Upwind" + o["components"]["nacelle"]["drivetrain"]["uptilt"] = 0.0 + # hub 6-vector [Ixx, Iyy, Izz, Ixy, Ixz, Iyz] with a nonzero Iyz + o["components"]["hub"]["elastic_properties_mb"]["system_inertia"] = [ + 2.0e6, 1.0e6, 1.0e6, 0.0, 0.0, 5.0e5, + ] + rna = read_windio_rna(_write(o, tmp_path, "hub_offdiag.yaml")) + assert rna.iyz == pytest.approx(5.0e5) + + def test_rna_hub_radius_projected_through_cone(tmp_path: pathlib.Path) -> None: """The hub radius is coned like the span (WISDEM Rtip*cos(precone)), so a coned rotor with a hub radius carries an extra Rhub*sin(cone) axial offset From 30fe99e049b347638af383bfeb26aed624742b5a Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Mon, 6 Jul 2026 06:02:55 +0900 Subject: [PATCH 17/17] fix: measure blade span from the root in the rotor lever reference_axis.z can be offset from zero (defined from a hub datum), but the rotor lever used the absolute z as the distance from the root, placing every section z[0] too far out and inflating the rotor inertia and coned CM offset. The mass integration already used the root-relative arc length, so only the lever was wrong. Subtract the root value: span = z - z[0]. A [20, 70] axis now gives the same rotor inertia and CM as [0, 50]. --- src/pybmodes/io/windio.py | 7 ++++++- tests/test_windio_rna.py | 24 ++++++++++++++++++++++++ 2 files changed, 30 insertions(+), 1 deletion(-) diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 70a835d..4017961 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -980,7 +980,12 @@ def _blade_span_mass_inertia( # unaffected. The polar lever is the in-plane distance from the shaft axis # ``r² + t²`` (so sweep and prebend are not dropped, issue #130); the # axial offset feeds the coned transverse term. - x_off, y_off, span = coords[0], coords[1], coords[2] + # Span measured from the blade root: reference_axis.z may be offset from + # zero (defined from a hub datum), so subtract the root value rather than + # treat the absolute z as the distance from the root, which would place + # every section too far out. + x_off, y_off = coords[0], coords[1] + span = coords[2] - coords[2][0] # Distance from the hub centre along the (coned) pitch axis: the blade # root sits at the hub radius and each section a further `span` outboard, # so the whole (hub_r + span) length projects through the cone, matching diff --git a/tests/test_windio_rna.py b/tests/test_windio_rna.py index f242c8c..31a970a 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -347,6 +347,30 @@ def _mk(y_off: float): assert swept.izz > straight.izz +def test_rna_blade_span_offset_from_root(tmp_path: pathlib.Path) -> None: + """reference_axis.z offset from zero (defined from a hub datum) is measured + relative to the root, so a [20, 70] axis gives the same rotor inertia and + CM as [0, 50] rather than placing every section 20 m too far out (Codex + review on #130).""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + def _mk(z0: float, z1: float): + o = _base_ontology() + o["components"]["hub"]["cone_angle"] = 0.15 + o["components"]["nacelle"]["drivetrain"]["uptilt"] = 0.10 + ra = o["components"]["blade"]["elastic_properties_mb"]["six_x_six"][ + "reference_axis" + ] + ra["z"] = {"grid": [0.0, 1.0], "values": [z0, z1]} + return read_windio_rna(_write(o, tmp_path, f"z_{z0}_{z1}.yaml")) + + base = _mk(0.0, 50.0) + offset = _mk(20.0, 70.0) + for attr in ("mass", "cm_axial", "ixx", "iyy", "izz"): + assert getattr(offset, attr) == pytest.approx(getattr(base, attr)) + + def test_rna_upwind_hub_offdiagonal_inertia_sign(tmp_path: pathlib.Path) -> None: """The WindIO hub frame has z up, so at zero tilt the hub tensor rotates by the identity and its off-diagonal Iyz passes through unchanged, even for an