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.

Author: SMI-Lab-Inha

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 = [

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