diff --git a/CHANGELOG.md b/CHANGELOG.md index 64537ac..a2810ed 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,7 +8,30 @@ 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_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 diff --git a/src/pybmodes/io/windio.py b/src/pybmodes/io/windio.py index 21741cc..4017961 100644 --- a/src/pybmodes/io/windio.py +++ b/src/pybmodes/io/windio.py @@ -67,6 +67,50 @@ 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 _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. @@ -688,6 +732,53 @@ def _finite_float(value: Any, what: str) -> float: return f +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). + ``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") + 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 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}, units={units}) resolves to more than 90 deg, which " + f"is non-physical." + ) + return cone_rad, uptilt_rad + + def _positive_mass(value: Any, what: str) -> float: """Coerce ``value`` to a positive, finite mass (kg), else raise.""" m = _finite_float(value, what) @@ -783,8 +874,11 @@ 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, cone_sin: float, +) -> tuple[float, 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 @@ -794,6 +888,19 @@ 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, 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 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 + 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) @@ -864,7 +971,36 @@ 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) + # 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. + # 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 + # 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 = axis_dist * 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 def read_windio_rna( @@ -873,6 +1009,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``. @@ -892,11 +1029,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 +1043,22 @@ 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. + + ``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). """ @@ -972,7 +1121,7 @@ def _nac_geom(key: str) -> float: return _finite_float(val, f"nacelle.drivetrain.{key}") overhang = _nac_geom("overhang") - uptilt = _nac_geom("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 --- @@ -989,6 +1138,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_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)) # --- Blades: assembly.number_of_blades x integrated span mass --- n_blades = _require_key(assembly, "number_of_blades", "assembly") @@ -997,6 +1161,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( @@ -1004,6 +1183,9 @@ def _nac_geom(key: str) -> float: f"{assembly.get('rotor_orientation')!r}." ) 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}" @@ -1017,9 +1199,42 @@ 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, + 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 is diag([I_polar, I_diam, I_diam]): the polar + # 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- + # 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. + # 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( + 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 # x (downwind-positive frame); the vertical hub position is # distance_tt_hub directly (= Twr2Shft + overhang*sin(uptilt)). @@ -1029,10 +1244,35 @@ def _nac_geom(key: str) -> float: dtype=float, ) + # 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_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_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 = [ (m_nac, r_nac, i_nac), (m_hub, apex, i_hub), - (m_blades, apex.copy(), np.zeros((3, 3))), + (m_blades, rotor_cm, i_blades), ] m_total = float(sum(m for m, _, _ in bodies)) 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 88ebf61..31a970a 100644 --- a/tests/test_windio_rna.py +++ b/tests/test_windio_rna.py @@ -96,15 +96,350 @@ 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] = 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) +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 + # 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. + 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_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_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_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 + + 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_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_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.""" + pytest.importorskip("yaml") + from pybmodes.io.windio import read_windio_rna + + o = _base_ontology() + o["components"]["hub"]["cone_angle"] = 200.0 + with pytest.raises(ValueError, match="non-physical"): + 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 +) -> 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 + 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_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 + 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 + 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) + 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 @@ -463,9 +798,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 ( @@ -481,7 +822,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