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.

Author: SMI-Lab-Inha

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

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)