From 489fb199db72ccda9a5c81b94fa316a7898fb0e4 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 30 Sep 2026 16:09:52 +0900 Subject: [PATCH 1/4] fix: lump HydroDyn filled-member ballast into deck-built floating platforms --- src/pybmodes/io/_hydrodyn_ballast.py | 362 +++++++++++++++++++++++++ src/pybmodes/models/_platform.py | 52 ++++ src/pybmodes/models/tower.py | 13 +- tests/test_hydrodyn_ballast.py | 388 +++++++++++++++++++++++++++ 4 files changed, 813 insertions(+), 2 deletions(-) create mode 100644 src/pybmodes/io/_hydrodyn_ballast.py create mode 100644 tests/test_hydrodyn_ballast.py diff --git a/src/pybmodes/io/_hydrodyn_ballast.py b/src/pybmodes/io/_hydrodyn_ballast.py new file mode 100644 index 0000000..7a5b758 --- /dev/null +++ b/src/pybmodes/io/_hydrodyn_ballast.py @@ -0,0 +1,362 @@ +# Copyright 2024-2026 Jae Hoon Seo +# Marine Structural Mechanics and Integrity Lab (SMI Lab), Inha University +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Filled-member (flooded ballast) mass properties from a HydroDyn ``.dat``. + +Some floating decks keep their water ballast out of ElastoDyn's +``PtfmMass`` and declare it instead as HydroDyn *filled member groups*. +The OC4 DeepCwind semi is the standard example, where ``PtfmMass`` is the +3.852e6 kg of steel and the 9.6e6 kg of ballast water lives in two fill +groups (upper and base columns). OpenFAST adds that ballast through +HydroDyn's Morison loads at run time, so a model assembled from the +ElastoDyn scalars alone is missing most of the platform mass, puts the +centre of gravity 5 m too high and loses its roll and pitch restoring. + +HydroDyn treats a closed fill group as body-fixed ballast. The fill mass, +centroid and inertia of each element are precomputed once +(``FloodedBallastPartSegmentCyl`` in ``Morison.f90``) and the internal +pressure datum follows the instantaneous high point of the ballast, so +the net fill load is the weight of a rigid mass at the fill centroid, +with no free-surface correction. This module reproduces those mass +properties so they can be lumped into the platform. + +Only the geometry needed for that is parsed, namely the member joints, +the cross-section property sets, the member table and the fill groups. +Both member-table layouts in circulation are accepted (with and without +the ``MSecGeom`` / ``MSpinOrient`` columns), as is the older single +``NPropSets`` property table. +""" + +from __future__ import annotations + +import pathlib +from dataclasses import dataclass, field + +import numpy as np + +from pybmodes.io.wamit_reader import ( + _is_fortran_default, + _parse_fortran_float, +) + +__all__ = ["FilledBallast", "read_filled_ballast"] + +#: Gauss-Legendre order for the along-axis integrals. The integrands are +#: polynomials of degree at most 4 in the axial coordinate (the inner +#: radius tapers linearly), so five points integrate them exactly. +_GAUSS_ORDER = 5 + + +@dataclass +class FilledBallast: + """Rigid-body mass properties of all filled members in a HydroDyn deck. + + Attributes + ---------- + mass : float + Total fill mass, kg. + cg : np.ndarray + Centre of mass ``(x, y, z)`` in the HydroDyn global frame (z up + from MSL), metres. + inertia_cg : np.ndarray + 3×3 inertia tensor about ``cg`` in kg·m², axes (x, y, z), i.e. + roll, pitch, yaw. Off-diagonal entries are the negated products + of inertia, ``I_ij = ∫ (|r|² δ_ij − r_i r_j) dm``. + member_mass : dict[int, float] + Fill mass of each filled member, keyed by ``MemberID``, for + diagnostics. + """ + + mass: float + cg: np.ndarray + inertia_cg: np.ndarray + member_mass: dict[int, float] = field(default_factory=dict) + + +@dataclass +class _Member: + member_id: int + joint1: int + joint2: int + prop1: int + prop2: int + sec_geom: int + + +def _tokens(line: str) -> list[str]: + return line.split("!", 1)[0].split() + + +def _count_after(lines: list[str], label: str) -> int | None: + """Index of the line whose second token is ``label``, or ``None``.""" + for i, raw in enumerate(lines): + parts = raw.split() + if len(parts) >= 2 and parts[1] == label: + return i + return None + + +def _table( + lines: list[str], label: str, path: pathlib.Path, +) -> tuple[list[str], list[list[str]]]: + """Return ``(header_tokens, rows)`` for the table counted by ``label``. + + HydroDyn tables are a count line, a column-name line, a units line + and then ``count`` data rows. A missing label means an empty table. + """ + i = _count_after(lines, label) + if i is None: + return [], [] + count_tok = lines[i].split()[0] + try: + n = int(count_tok) + except ValueError as err: + raise ValueError( + f"{label} in {path} is not an integer: {count_tok!r}" + ) from err + if n < 0: + raise ValueError(f"{label} in {path} is negative: {n}") + header = _tokens(lines[i + 1]) if i + 1 < len(lines) else [] + rows = [] + for j in range(i + 3, i + 3 + n): + if j >= len(lines): + raise ValueError( + f"{path} ends inside the {label} table " + f"(expected {n} rows)" + ) + rows.append(_tokens(lines[j])) + return header, rows + + +def _float(tok: str, what: str, path: pathlib.Path) -> float: + try: + return _parse_fortran_float(tok) + except ValueError as err: + raise ValueError(f"Malformed {what} in {path}: {tok!r}") from err + + +def _int(tok: str, what: str, path: pathlib.Path) -> int: + try: + return int(tok) + except ValueError as err: + raise ValueError(f"Malformed {what} in {path}: {tok!r}") from err + + +def _water_density(lines: list[str], path: pathlib.Path) -> float: + """``WtrDens`` from the deck, or the paired SeaState file default.""" + i = _count_after(lines, "WtrDens") + if i is not None: + tok = lines[i].split()[0] + if not _is_fortran_default(tok): + return _float(tok, "WtrDens", path) + return 1025.0 + + +def read_filled_ballast( + dat_path: str | pathlib.Path, + *, + water_density: float | None = None, +) -> FilledBallast | None: + """Mass properties of the filled members declared in a HydroDyn deck. + + Parameters + ---------- + dat_path + HydroDyn ``.dat`` file. + water_density + Density used for fill groups whose ``FillDens`` is ``DEFAULT``. + When omitted, the deck's own ``WtrDens`` is used, falling back to + 1025 kg/m³ (the SeaState default) for decks that delegate it. + + Returns + ------- + FilledBallast or None + ``None`` when the deck declares no fill groups or none of them + holds any fluid. + + Each filled member is a straight frustum whose outer diameter and + wall thickness vary linearly between its two property sets, as in + HydroDyn. The fill runs from the lower end up to ``FillFSLoc`` + (metres, relative to MSL), measured along the member axis exactly as + HydroDyn does. Rectangular members are rejected with a clear error + rather than silently skipped. + """ + path = pathlib.Path(dat_path) + if not path.is_file(): + raise FileNotFoundError(f"HydroDyn .dat not found at {path}") + lines = path.read_text(encoding="utf-8", errors="replace").splitlines() + + _, fill_rows = _table(lines, "NFillGroups", path) + if not fill_rows: + return None + + if water_density is None: + water_density = _water_density(lines, path) + if not np.isfinite(water_density) or water_density <= 0.0: + raise ValueError( + f"water_density must be finite and positive, got {water_density!r}" + ) + + _, joint_rows = _table(lines, "NJoints", path) + joints: dict[int, np.ndarray] = {} + for row in joint_rows: + if len(row) < 4: + raise ValueError(f"Short joint row in {path}: {row!r}") + jid = _int(row[0], "JointID", path) + joints[jid] = np.array( + [_float(t, f"joint {jid} coordinate", path) for t in row[1:4]] + ) + + props: dict[int, tuple[float, float]] = {} + label = "NPropSetsCyl" if _count_after(lines, "NPropSetsCyl") is not None \ + else "NPropSets" + for row in _table(lines, label, path)[1]: + if len(row) < 3: + raise ValueError(f"Short property-set row in {path}: {row!r}") + pid = _int(row[0], "PropSetID", path) + diam = _float(row[1], f"PropD of set {pid}", path) + thick = _float(row[2], f"PropThck of set {pid}", path) + if diam <= 0.0 or thick < 0.0 or 2.0 * thick > diam: + raise ValueError( + f"Property set {pid} in {path} has a non-physical section " + f"(PropD = {diam}, PropThck = {thick})" + ) + props[pid] = (diam, thick) + + header, member_rows = _table(lines, "NMembers", path) + cols = {name: k for k, name in enumerate(header)} + for need in ("MemberID", "MJointID1", "MJointID2", + "MPropSetID1", "MPropSetID2"): + if need not in cols: + raise ValueError( + f"Member table in {path} has no {need} column; " + f"header was {header!r}" + ) + members: dict[int, _Member] = {} + for row in member_rows: + if len(row) <= max(cols[c] for c in ("MemberID", "MJointID1", + "MJointID2", "MPropSetID1", + "MPropSetID2")): + raise ValueError(f"Short member row in {path}: {row!r}") + mid = _int(row[cols["MemberID"]], "MemberID", path) + sec_geom = 1 + if "MSecGeom" in cols: + sec_geom = _int(row[cols["MSecGeom"]], "MSecGeom", path) + members[mid] = _Member( + member_id=mid, + joint1=_int(row[cols["MJointID1"]], "MJointID1", path), + joint2=_int(row[cols["MJointID2"]], "MJointID2", path), + prop1=_int(row[cols["MPropSetID1"]], "MPropSetID1", path), + prop2=_int(row[cols["MPropSetID2"]], "MPropSetID2", path), + sec_geom=sec_geom, + ) + + xg, wg = np.polynomial.legendre.leggauss(_GAUSS_ORDER) + mass = 0.0 + first = np.zeros(3) + second = np.zeros((3, 3)) # ∫ (|r|² I − r rᵀ) dm about the origin + member_mass: dict[int, float] = {} + + for row in fill_rows: + n_mem = _int(row[0], "FillNumM", path) + if n_mem < 1 or len(row) < n_mem + 3: + raise ValueError(f"Malformed fill-group row in {path}: {row!r}") + ids = [_int(t, "FillMList entry", path) for t in row[1:1 + n_mem]] + fs_z = _float(row[1 + n_mem], "FillFSLoc", path) + dens_tok = row[2 + n_mem] + dens = water_density if _is_fortran_default(dens_tok) \ + else _float(dens_tok, "FillDens", path) + if dens < 0.0: + raise ValueError(f"Negative FillDens in {path}: {dens}") + + for mid in ids: + if mid not in members: + raise ValueError( + f"Fill group in {path} names member {mid}, which is not " + f"in the member table" + ) + mem = members[mid] + if mem.sec_geom != 1: + raise NotImplementedError( + f"Filled member {mid} in {path} is rectangular " + f"(MSecGeom = {mem.sec_geom}); only cylindrical filled " + f"members are supported" + ) + for jid in (mem.joint1, mem.joint2): + if jid not in joints: + raise ValueError( + f"Member {mid} in {path} uses undefined joint {jid}" + ) + for pid in (mem.prop1, mem.prop2): + if pid not in props: + raise ValueError( + f"Member {mid} in {path} uses undefined property " + f"set {pid}" + ) + + # Orient from the lower end, as HydroDyn does (Za <= Zb). + p_a, p_b = joints[mem.joint1], joints[mem.joint2] + (d_a, t_a), (d_b, t_b) = props[mem.prop1], props[mem.prop2] + if p_b[2] < p_a[2]: + p_a, p_b = p_b, p_a + d_a, t_a, d_b, t_b = d_b, t_b, d_a, t_a + axis = p_b - p_a + length = float(np.linalg.norm(axis)) + if length <= 0.0: + raise ValueError(f"Member {mid} in {path} has zero length") + e = axis / length + rin_a = 0.5 * d_a - t_a + rin_b = 0.5 * d_b - t_b + + if fs_z >= p_b[2]: + l_fill = length + elif p_a[2] >= fs_z: + l_fill = 0.0 + else: + # Only reached when p_b[2] > fs_z > p_a[2], so e[2] > 0. + l_fill = (fs_z - p_a[2]) / e[2] + if l_fill <= 0.0 or dens == 0.0: + member_mass[mid] = member_mass.get(mid, 0.0) + continue + + s = 0.5 * l_fill * (xg + 1.0) + w = 0.5 * l_fill * wg + rin = rin_a + (rin_b - rin_a) * s / length + dm = dens * np.pi * rin ** 2 * w # mass per Gauss slice + pos = p_a[None, :] + s[:, None] * e[None, :] + + m_mem = float(dm.sum()) + mass += m_mem + member_mass[mid] = member_mass.get(mid, 0.0) + m_mem + first += dm @ pos + # Point-mass part plus each slice's own disc inertia, + # (r²/4)(I + e eᵀ) per unit mass for a thin disc of radius r + # normal to e. The slices' axial extent is recovered exactly + # by the quadrature through the point-mass part. + r2 = np.einsum("ki,ki->k", pos, pos) + second += (dm * r2).sum() * np.eye(3) - np.einsum( + "k,ki,kj->ij", dm, pos, pos, + ) + second += float((dm * rin ** 2 / 4.0).sum()) * ( + np.eye(3) + np.outer(e, e) + ) + + if mass <= 0.0: + return None + cg = first / mass + inertia_cg = second - mass * (float(cg @ cg) * np.eye(3) - np.outer(cg, cg)) + return FilledBallast( + mass=mass, cg=cg, inertia_cg=inertia_cg, member_mass=member_mass, + ) diff --git a/src/pybmodes/models/_platform.py b/src/pybmodes/models/_platform.py index 877caca..e028cd5 100644 --- a/src/pybmodes/models/_platform.py +++ b/src/pybmodes/models/_platform.py @@ -142,6 +142,58 @@ def _platform_inertia_matrix(ptfm: dict[str, float]) -> np.ndarray: return i_mat +def _add_filled_ballast( + ptfm: dict[str, float], + i_mat: np.ndarray, + hydrodyn_dat_path: pathlib.Path | str, +) -> tuple[dict[str, float], np.ndarray]: + """Lump the HydroDyn filled-member ballast into the platform. + + Returns an updated copy of the ``Ptfm*`` scalars (``PtfmMass`` and + the ``PtfmCM*`` centre of mass) and the 6×6 inertia matrix about the + combined centre of mass. When the deck declares no fill, both come + back unchanged. + + OpenFAST applies filled-member ballast through HydroDyn rather than + ElastoDyn, so decks such as the OC4 DeepCwind semi carry only the + steel in ``PtfmMass``. Leaving the ballast out misses most of the + platform mass, raises its centre of gravity and, once the weight + restoring is formed from that mass, leaves roll and pitch unstable. + HydroDyn models a closed fill group as body-fixed ballast (see + :mod:`pybmodes.io._hydrodyn_ballast`), which is exactly a rigid mass + on the platform. The combined rotational block is a full 3×3 tensor; + for a symmetric platform its products of inertia vanish. + """ + import numpy as np + + from pybmodes.io._hydrodyn_ballast import read_filled_ballast + + ballast = read_filled_ballast(hydrodyn_dat_path) + if ballast is None: + return ptfm, i_mat + + m_s = float(ptfm["PtfmMass"]) + c_s = np.array([ptfm["PtfmCMxt"], ptfm["PtfmCMyt"], ptfm["PtfmCMzt"]]) + m_b = ballast.mass + m = m_s + m_b + c = (m_s * c_s + m_b * ballast.cg) / m + + def shift(mass: float, d: np.ndarray) -> np.ndarray: + return mass * (float(d @ d) * np.eye(3) - np.outer(d, d)) + + rot = ( + np.asarray(i_mat, dtype=float)[3:, 3:] + shift(m_s, c_s - c) + + ballast.inertia_cg + shift(m_b, ballast.cg - c) + ) + out = dict(ptfm) + out["PtfmMass"] = m + out["PtfmCMxt"], out["PtfmCMyt"], out["PtfmCMzt"] = (float(v) for v in c) + new_i = np.zeros((6, 6)) + new_i[0, 0] = new_i[1, 1] = new_i[2, 2] = m + new_i[3:, 3:] = 0.5 * (rot + rot.T) + return out, new_i + + #: Standard gravity (m/s²), the OpenFAST ``Gravity`` default. STANDARD_GRAVITY = 9.80665 diff --git a/src/pybmodes/models/tower.py b/src/pybmodes/models/tower.py index 63d41b7..43ae806 100644 --- a/src/pybmodes/models/tower.py +++ b/src/pybmodes/models/tower.py @@ -40,6 +40,7 @@ from pybmodes.io.sec_props import SectionProperties from pybmodes.models._pipeline import run_fem from pybmodes.models._platform import ( + _add_filled_ballast, _gravitational_restoring, _platform_inertia_matrix, _scan_platform_fields, @@ -1028,8 +1029,13 @@ def from_elastodyn_with_mooring( A_inf = wamit.A_inf C_hst = wamit.C_hst - M = ptfm["PtfmMass"] i_mat = _platform_inertia_matrix(ptfm) + # Ballast declared as HydroDyn filled members (the OC4 semi keeps + # all of its water ballast there) is platform mass that ElastoDyn + # never sees. Lump it in before the CM and weight terms are formed. + if hydrodyn_dat_path is not None: + ptfm, i_mat = _add_filled_ballast(ptfm, i_mat, hydrodyn_dat_path) + M = ptfm["PtfmMass"] # BModes file convention for these scalars (see the OC3 Hywind # sample BMI in ``src/pybmodes/_examples/sample_inputs/ @@ -1396,8 +1402,11 @@ def _tower_grid_od_wt( blade = read_elastodyn_blade(bp) rna_tip = _tower_top_assembly_mass(main, blade) ptfm = _scan_platform_fields(ed_path) - M = ptfm["PtfmMass"] i_mat = _platform_inertia_matrix(ptfm) + # Filled-member ballast, as in from_elastodyn_with_mooring. + if hydrodyn_dat is not None: + ptfm, i_mat = _add_filled_ballast(ptfm, i_mat, hydrodyn_dat) + M = ptfm["PtfmMass"] cm_pform = -ptfm["PtfmCMzt"] cm_x, cm_y = ptfm["PtfmCMxt"], ptfm["PtfmCMyt"] ref_msl = -ptfm["PtfmRefzt"] diff --git a/tests/test_hydrodyn_ballast.py b/tests/test_hydrodyn_ballast.py new file mode 100644 index 0000000..a90e41a --- /dev/null +++ b/tests/test_hydrodyn_ballast.py @@ -0,0 +1,388 @@ +"""Filled-member ballast from HydroDyn decks. + +HydroDyn treats a closed fill group as body-fixed ballast. The fill mass, +centroid and inertia are precomputed per element and the internal +pressure datum follows the ballast's high point, so the net load is the +weight of a rigid mass (``FloodedBallastPartSegmentCyl`` and the fill +terms in ``Morison.f90``). The self-contained tests check the reader +against the closed-form mass properties of solid cylinders and frustums. + +The integration anchor is the OC4 DeepCwind semi, which keeps all of its +water ballast in two fill groups while ElastoDyn's ``PtfmMass`` holds only +the 3.852e6 kg of steel. Robertson et al. (2014), *Definition of the +Semisubmersible Floating System for Phase II of OC4*, NREL/TP-5000-60601, +give the ballasted platform as 1.3473e7 kg with its centre of mass +13.46 m below the still water line. +""" + +from __future__ import annotations + +import pathlib +import warnings + +import numpy as np +import pytest + +from pybmodes.io._hydrodyn_ballast import read_filled_ballast +from pybmodes.models._platform import _add_filled_ballast + +RHO = 1025.0 + +REPO_ROOT = pathlib.Path(__file__).resolve().parents[1] +OC4_DECKS = ( + REPO_ROOT / "external" / "OpenFAST_files" / "r-test" / "glue-codes" + / "openfast" / "5MW_OC4Semi_WSt_WavesWN" +) +OC4_PREFIX = "NRELOffshrBsline5MW_OC4DeepCwindSemi_" +OC4_ELASTODYN = OC4_DECKS / f"{OC4_PREFIX}ElastoDyn.dat" +OC4_MOORDYN = OC4_DECKS / f"{OC4_PREFIX}MoorDyn.dat" +OC4_HYDRODYN = OC4_DECKS / f"{OC4_PREFIX}HydroDyn.dat" + +# Robertson et al. (2014) NREL/TP-5000-60601, platform including ballast. +OC4_PLATFORM_MASS = 1.3473e7 +OC4_PLATFORM_CM_Z = -13.46 + + +def _deck( + tmp_path: pathlib.Path, + joints: list[tuple[int, float, float, float]], + props: list[tuple[int, float, float]], + members: list[tuple[int, int, int, int, int]], + fills: list[str], + *, + legacy: bool = False, + sec_geom: dict[int, int] | None = None, + wtr_dens: str | None = None, +) -> pathlib.Path: + """Write the HydroDyn sections the ballast reader looks at.""" + sec_geom = sec_geom or {} + out = ["------- HydroDyn Input File (synthetic) -------"] + if wtr_dens is not None: + out.append(f"{wtr_dens} WtrDens - Water density (kg/m^3)") + out.append("---------------------- MEMBER JOINTS ----------------------") + out.append(f"{len(joints)} NJoints - Number of joints (-)") + out.append("JointID Jointxi Jointyi Jointzi JointAxID JointOvrlp") + out.append(" (-) (m) (m) (m) (-) (switch)") + out += [f"{j} {x} {y} {z} 1 0" for j, x, y, z in joints] + if legacy: + out.append("------- MEMBER CROSS-SECTION PROPERTIES -------") + out.append(f"{len(props)} NPropSets - Number of member property sets") + else: + out.append("------- CYLINDRICAL MEMBER CROSS-SECTION PROPERTIES -------") + out.append(f"{len(props)} NPropSetsCyl - Number of cylindrical sets") + out.append("PropSetID PropD PropThck") + out.append(" (-) (m) (m)") + out += [f"{p} {d} {t}" for p, d, t in props] + if not legacy: + out.append("------- RECTANGULAR MEMBER CROSS-SECTION PROPERTIES -------") + out.append("0 NPropSetsRec - Number of rectangular sets") + out.append("PropSetID PropA PropB PropThck") + out.append(" (-) (m) (m) (m)") + out.append("-------------------- MEMBERS --------------------") + out.append(f"{len(members)} NMembers - Number of members (-)") + if legacy: + out.append("MemberID MJointID1 MJointID2 MPropSetID1 MPropSetID2 " + "MDivSize MCoefMod PropPot") + out.append(" (-) (-) (-) (-) (-) (m) (switch) (flag)") + out += [f"{m} {a} {b} {p} {q} 1.0 1 TRUE" + for m, a, b, p, q in members] + else: + out.append("MemberID MJointID1 MJointID2 MPropSetID1 MPropSetID2 " + "MSecGeom MSpinOrient MDivSize MCoefMod MHstLMod PropPot") + out.append(" (-) (-) (-) (-) (-) (switch) (deg) (m) (switch) " + "(switch) (flag)") + out += [f"{m} {a} {b} {p} {q} {sec_geom.get(m, 1)} 0 1.0 1 1 TRUE" + for m, a, b, p, q in members] + out.append("---------------------- FILLED MEMBERS ----------------------") + out.append(f"{len(fills)} NFillGroups - Number of filled member groups") + out.append("FillNumM FillMList FillFSLoc FillDens") + out.append(" (-) (-) (m) (kg/m^3)") + out += fills + out.append("---------------------- MARINE GROWTH ----------------------") + out.append("0 NMGDepths - Number of marine-growth depths specified") + path = tmp_path / "hd.dat" + path.write_text("\n".join(out) + "\n", encoding="utf-8") + return path + + +def _vertical_column(tmp_path, *, fs: float, **kw) -> pathlib.Path: + # Outer D = 10 m, wall 0.5 m, so inner radius 4.5 m; z from -20 to +10. + return _deck( + tmp_path, + joints=[(1, 3.0, -2.0, -20.0), (2, 3.0, -2.0, 10.0)], + props=[(1, 10.0, 0.5)], + members=[(1, 1, 2, 1, 1)], + fills=[f"1 1 {fs} {RHO}"], + **kw, + ) + + +class TestClosedForms: + def test_partially_filled_vertical_cylinder(self, tmp_path): + b = read_filled_ballast(_vertical_column(tmp_path, fs=-8.0)) + r, h = 4.5, 12.0 + m = RHO * np.pi * r ** 2 * h + assert b.mass == pytest.approx(m, rel=1e-12) + np.testing.assert_allclose(b.cg, [3.0, -2.0, -14.0], atol=1e-9) + # Solid cylinder about its own centroid (axis along z). + i_perp = m * (3.0 * r ** 2 + h ** 2) / 12.0 + np.testing.assert_allclose( + b.inertia_cg, np.diag([i_perp, i_perp, m * r ** 2 / 2.0]), + rtol=1e-10, atol=1e-6 * i_perp, + ) + + def test_member_ends_are_order_independent(self, tmp_path): + a = read_filled_ballast(_vertical_column(tmp_path, fs=-8.0)) + path = _deck( + tmp_path, + joints=[(1, 3.0, -2.0, -20.0), (2, 3.0, -2.0, 10.0)], + props=[(1, 10.0, 0.5)], + members=[(1, 2, 1, 1, 1)], + fills=[f"1 1 -8.0 {RHO}"], + ) + b = read_filled_ballast(path) + assert b.mass == pytest.approx(a.mass, rel=1e-12) + np.testing.assert_allclose(b.inertia_cg, a.inertia_cg, rtol=1e-10) + + def test_free_surface_above_member_fills_it_completely(self, tmp_path): + b = read_filled_ballast(_vertical_column(tmp_path, fs=50.0)) + assert b.mass == pytest.approx(RHO * np.pi * 4.5 ** 2 * 30.0) + + def test_free_surface_below_member_holds_nothing(self, tmp_path): + assert read_filled_ballast(_vertical_column(tmp_path, fs=-25.0)) is None + + def test_tapered_frustum_volume(self, tmp_path): + # Inner radius tapers 5.0 -> 2.0 m over a fully filled 9 m length. + path = _deck( + tmp_path, + joints=[(1, 0.0, 0.0, -30.0), (2, 0.0, 0.0, -21.0)], + props=[(1, 10.2, 0.1), (2, 4.2, 0.1)], + members=[(1, 1, 2, 1, 2)], + fills=[f"1 1 0.0 {RHO}"], + ) + b = read_filled_ballast(path) + r1, r2, h = 5.0, 2.0, 9.0 + vol = np.pi * h / 3.0 * (r1 ** 2 + r1 * r2 + r2 ** 2) + z_c = h * (r1 ** 2 + 2 * r1 * r2 + 3 * r2 ** 2) / ( + 4 * (r1 ** 2 + r1 * r2 + r2 ** 2)) + assert b.mass == pytest.approx(RHO * vol, rel=1e-12) + assert b.cg[2] == pytest.approx(-30.0 + z_c, rel=1e-12) + # Axial moment of a frustum, (3/10) m (r1^5 - r2^5)/(r1^3 - r2^3). + izz = 0.3 * RHO * vol * (r1 ** 5 - r2 ** 5) / (r1 ** 3 - r2 ** 3) + assert b.inertia_cg[2, 2] == pytest.approx(izz, rel=1e-12) + + def test_horizontal_member_is_all_or_nothing(self, tmp_path): + def deck(fs): + return _deck( + tmp_path, + joints=[(1, -5.0, 0.0, -10.0), (2, 5.0, 0.0, -10.0)], + props=[(1, 2.0, 0.0)], + members=[(1, 1, 2, 1, 1)], + fills=[f"1 1 {fs} {RHO}"], + ) + full = read_filled_ballast(deck(-9.0)) + assert full.mass == pytest.approx(RHO * np.pi * 10.0) + m, r, L = full.mass, 1.0, 10.0 + np.testing.assert_allclose( + np.diag(full.inertia_cg), + [m * r ** 2 / 2, m * (3 * r ** 2 + L ** 2) / 12, + m * (3 * r ** 2 + L ** 2) / 12], + rtol=1e-10, + ) + # HydroDyn tests FillFSLoc >= Zb first, so a surface level with + # the member still counts it as full. + assert read_filled_ballast(deck(-10.0)).mass == pytest.approx(full.mass) + assert read_filled_ballast(deck(-10.5)) is None + + def test_inclined_member_fills_along_its_axis(self, tmp_path): + # 45 degree member from z = -10 to z = 0; free surface at -4 gives + # an along-axis fill length of 6·sqrt(2). + path = _deck( + tmp_path, + joints=[(1, 0.0, 0.0, -10.0), (2, 10.0, 0.0, 0.0)], + props=[(1, 2.0, 0.0)], + members=[(1, 1, 2, 1, 1)], + fills=[f"1 1 -4.0 {RHO}"], + ) + b = read_filled_ballast(path) + assert b.mass == pytest.approx(RHO * np.pi * 6.0 * np.sqrt(2.0)) + np.testing.assert_allclose(b.cg, [3.0, 0.0, -7.0], atol=1e-12) + + def test_two_groups_combine_by_parallel_axis(self, tmp_path): + path = _deck( + tmp_path, + joints=[(1, 10.0, 0.0, -20.0), (2, 10.0, 0.0, 0.0), + (3, -10.0, 0.0, -20.0), (4, -10.0, 0.0, 0.0)], + props=[(1, 2.0, 0.0)], + members=[(1, 1, 2, 1, 1), (2, 3, 4, 1, 1)], + fills=[f"1 1 -10.0 {RHO}", f"1 2 -10.0 {RHO}"], + ) + b = read_filled_ballast(path) + m1 = RHO * np.pi * 10.0 + assert b.mass == pytest.approx(2 * m1) + np.testing.assert_allclose(b.cg, [0.0, 0.0, -15.0], atol=1e-12) + i_own_perp = m1 * (3.0 + 100.0) / 12.0 + assert b.inertia_cg[1, 1] == pytest.approx( + 2 * (i_own_perp + m1 * 100.0), rel=1e-12) + assert b.inertia_cg[0, 0] == pytest.approx(2 * i_own_perp, rel=1e-12) + assert b.inertia_cg[2, 2] == pytest.approx( + 2 * (m1 / 2.0 + m1 * 100.0), rel=1e-12) + assert set(b.member_mass) == {1, 2} + + +class TestDeckHandling: + def test_no_fill_groups_returns_none(self, tmp_path): + path = _deck( + tmp_path, joints=[(1, 0, 0, -1), (2, 0, 0, 1)], + props=[(1, 1.0, 0.1)], members=[(1, 1, 2, 1, 1)], fills=[], + ) + assert read_filled_ballast(path) is None + + def test_legacy_layout(self, tmp_path): + new = read_filled_ballast(_vertical_column(tmp_path, fs=-8.0)) + old = read_filled_ballast( + _vertical_column(tmp_path, fs=-8.0, legacy=True)) + assert old.mass == pytest.approx(new.mass, rel=1e-12) + + @pytest.mark.parametrize("token", ["DEFAULT", '"default"']) + def test_default_fill_density_is_the_water_density(self, tmp_path, token): + path = _deck( + tmp_path, joints=[(1, 0, 0, -10), (2, 0, 0, 0)], + props=[(1, 2.0, 0.0)], members=[(1, 1, 2, 1, 1)], + fills=[f"1 1 0.0 {token}"], wtr_dens="1000", + ) + b = read_filled_ballast(path) + assert b.mass == pytest.approx(1000.0 * np.pi * 10.0) + b2 = read_filled_ballast(path, water_density=1025.0) + assert b2.mass == pytest.approx(1025.0 * np.pi * 10.0) + + def test_rectangular_member_is_rejected(self, tmp_path): + path = _vertical_column(tmp_path, fs=-8.0, sec_geom={1: 2}) + with pytest.raises(NotImplementedError, match="rectangular"): + read_filled_ballast(path) + + def test_unknown_member_is_rejected(self, tmp_path): + path = _deck( + tmp_path, joints=[(1, 0, 0, -10), (2, 0, 0, 0)], + props=[(1, 2.0, 0.0)], members=[(1, 1, 2, 1, 1)], + fills=[f"1 7 0.0 {RHO}"], + ) + with pytest.raises(ValueError, match="member 7"): + read_filled_ballast(path) + + def test_non_physical_section_is_rejected(self, tmp_path): + path = _deck( + tmp_path, joints=[(1, 0, 0, -10), (2, 0, 0, 0)], + props=[(1, 2.0, 1.5)], members=[(1, 1, 2, 1, 1)], + fills=[f"1 1 0.0 {RHO}"], + ) + with pytest.raises(ValueError, match="non-physical"): + read_filled_ballast(path) + + @pytest.mark.parametrize("bad", [0.0, -1.0, float("nan"), float("inf")]) + def test_bad_water_density_is_rejected(self, tmp_path, bad): + path = _vertical_column(tmp_path, fs=-8.0) + with pytest.raises(ValueError, match="water_density"): + read_filled_ballast(path, water_density=bad) + + def test_missing_file(self, tmp_path): + with pytest.raises(FileNotFoundError): + read_filled_ballast(tmp_path / "absent.dat") + + +class TestMergeIntoPlatform: + def _ptfm(self): + return { + "PtfmMass": 4.0e6, "PtfmCMxt": 0.0, "PtfmCMyt": 0.0, + "PtfmCMzt": -8.0, "PtfmRIner": 2.0e9, "PtfmPIner": 2.0e9, + "PtfmYIner": 3.0e9, + } + + def _i_mat(self, p): + return np.diag([p["PtfmMass"]] * 3 + + [p["PtfmRIner"], p["PtfmPIner"], p["PtfmYIner"]]) + + def test_combined_mass_cm_and_parallel_axis(self, tmp_path): + path = _vertical_column(tmp_path, fs=-8.0) + b = read_filled_ballast(path) + p = self._ptfm() + out, i_mat = _add_filled_ballast(p, self._i_mat(p), path) + + m = p["PtfmMass"] + b.mass + c_s = np.array([0.0, 0.0, -8.0]) + c = (p["PtfmMass"] * c_s + b.mass * b.cg) / m + assert out["PtfmMass"] == pytest.approx(m) + np.testing.assert_allclose( + [out["PtfmCMxt"], out["PtfmCMyt"], out["PtfmCMzt"]], c) + np.testing.assert_allclose(np.diag(i_mat)[:3], [m] * 3) + + def shift(mass, d): + return mass * (d @ d * np.eye(3) - np.outer(d, d)) + expect = (np.diag([2.0e9, 2.0e9, 3.0e9]) + shift(p["PtfmMass"], c_s - c) + + b.inertia_cg + shift(b.mass, b.cg - c)) + np.testing.assert_allclose(i_mat[3:, 3:], expect, rtol=1e-12) + # An off-axis column gives a genuine roll-yaw product of inertia. + assert i_mat[3, 5] != 0.0 + np.testing.assert_allclose(i_mat, i_mat.T) + # The input dict is not mutated. + assert p["PtfmMass"] == 4.0e6 + + def test_no_fill_leaves_platform_unchanged(self, tmp_path): + path = _deck( + tmp_path, joints=[(1, 0, 0, -1), (2, 0, 0, 1)], + props=[(1, 1.0, 0.1)], members=[(1, 1, 2, 1, 1)], fills=[], + ) + p = self._ptfm() + i0 = self._i_mat(p) + out, i_mat = _add_filled_ballast(p, i0, path) + assert out is p and i_mat is i0 + + +_needs_oc4 = pytest.mark.skipif( + not (OC4_ELASTODYN.is_file() and OC4_MOORDYN.is_file() + and OC4_HYDRODYN.is_file()), + reason=f"OC4 DeepCwind r-test decks not present under {OC4_DECKS}", +) + + +@pytest.mark.integration +@_needs_oc4 +class TestOC4DeepCwind: + def _tower(self): + from pybmodes.models import Tower + + return Tower.from_elastodyn_with_mooring( + OC4_ELASTODYN, OC4_MOORDYN, OC4_HYDRODYN, + ) + + def test_ballast_matches_the_published_platform(self): + ps = self._tower()._bmi.support + assert ps.mass_pform == pytest.approx(OC4_PLATFORM_MASS, rel=1e-3) + assert -ps.cm_pform == pytest.approx(OC4_PLATFORM_CM_Z, abs=0.01) + + def test_roll_and_pitch_restoring_is_positive(self): + k = self._tower()._bmi.support.hydro_K + assert k[3, 3] > 0.0 and k[4, 4] > 0.0 + + def test_all_six_rigid_body_modes_are_returned(self): + with warnings.catch_warnings(): + warnings.simplefilter("error", UserWarning) + res = self._tower().run(n_modes=8, check_model=False) + assert set(res.mode_labels[:6]) == { + "surge", "sway", "heave", "roll", "pitch", "yaw", + } + + def test_heave_matches_the_single_dof_estimate(self): + """Heave barely couples to the tower, so it should sit close to + sqrt(C33 / (m_total + A33)) / 2π with the ballasted mass.""" + tower = self._tower() + ps = tower._bmi.support + res = tower.run(n_modes=8, check_model=False) + f = dict(zip(res.mode_labels, res.frequencies))["heave"] + m_tower = float(np.trapezoid( + tower._sp.mass_den, tower._sp.span_loc)) * ( + tower._bmi.radius + ps.draft) + m_total = ps.mass_pform + m_tower + tower._bmi.tip_mass.mass + k33 = ps.hydro_K[2, 2] + ps.mooring_K[2, 2] + f_est = np.sqrt(k33 / (m_total + ps.hydro_M[2, 2])) / (2 * np.pi) + assert f == pytest.approx(f_est, rel=0.02) From 30d1503d9407713592ab39159e066c4077e66a34 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 30 Sep 2026 16:23:27 +0900 Subject: [PATCH 2/4] chore: regenerate floating samples with weight restoring and fill ballast The OC4 DeepCwind sample now carries its ballasted platform and positive roll and pitch restoring; the IEA-15, IEA-22 and UPSCALE semis pick up the weight term in hydro_K. A default-run test keeps every bundled floater statically stable with all six rigid-body modes labelled. --- CHANGELOG.md | 32 ++++++++++++++++--- VALIDATION.md | 4 +++ docs/conventions.rst | 6 ++++ .../01_nrel5mw_land/README.md | 8 ++--- .../02_nrel5mw_oc3monopile/README.md | 8 ++--- .../03_iea34_land/README.md | 8 ++--- .../04_iea10_monopile/README.md | 8 ++--- .../05_iea15_monopile/README.md | 8 ++--- .../06_iea22_monopile/README.md | 8 ++--- .../07_nrel5mw_oc3hywind_spar/README.md | 6 ++-- .../08_nrel5mw_oc4semi_tower.bmi | 20 ++++++------ .../08_nrel5mw_oc4semi/README.md | 8 ++--- .../09_iea15_umainesemi_tower.bmi | 4 +-- .../09_iea15_umainesemi/README.md | 8 ++--- .../10_iea22_semi/10_iea22_semi_tower.bmi | 12 +++---- .../10_iea22_semi/README.md | 6 ++-- .../11_upscale25_centraltower_tower.bmi | 14 ++++---- .../11_upscale25_centraltower/README.md | 8 ++--- src/pybmodes/models/_pipeline.py | 7 ++-- tests/test_floating_samples_spectra.py | 20 ++++++++++++ tests/test_hydrodyn_ballast.py | 4 +-- 21 files changed, 130 insertions(+), 77 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 83b7465..432498d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -81,11 +81,33 @@ The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). `Tower.run` on a free-free floating model now warns (`UserWarning`) when `hydro_K + mooring_K` has a negative direction, naming the dominant DOFs: the platform is statically unstable there and the - general solver leaves those modes out of the result. The bundled - `08_nrel5mw_oc4semi` sample triggers it — the OC4 DeepCwind deck keeps - its ballast in HydroDyn filled members rather than `PtfmMass`, which - the platform model does not include, so its roll / pitch restoring is - negative and the modes labelled roll / pitch there are tower bending. + general solver leaves those modes out of the result. +- **Ballast held in HydroDyn filled members was missing from floating + models.** Some decks keep their water ballast out of ElastoDyn's + `PtfmMass` and declare it as HydroDyn filled-member groups, which + OpenFAST applies through HydroDyn at run time. The OC4 DeepCwind semi + is the standard case, with 3.852e6 kg of steel in `PtfmMass` and + 9.62e6 kg of water in two fill groups. `Tower.from_elastodyn_with_mooring` + and the deck-backed tier of `Tower.from_windio_floating` now read the + joints, cross-sections, members and fill groups of the HydroDyn deck + and lump the fill into the platform as the rigid, body-fixed ballast + HydroDyn models it as: mass, centre of mass and full 3×3 inertia about + the combined centre of mass. On OC4 this gives 1.3476e7 kg with the + centre of mass 13.460 m below MSL, against 1.3473e7 kg and 13.46 m in + Robertson et al. (2014), NREL/TP-5000-60601. Roll / pitch restoring + goes from −3.8e8 to +9.85e8 N·m/rad and all six rigid-body modes come + back (surge / sway 0.0088, heave 0.0579, roll / pitch 0.0390, yaw + 0.0125 Hz); before, the modes labelled roll and pitch were tower + bending. Decks without fill groups are unaffected. Only cylindrical + filled members are supported; a rectangular one raises + `NotImplementedError` rather than being skipped. +- **Bundled floating samples regenerated.** `08_nrel5mw_oc4semi` now + carries the ballasted platform above (its first tower FA mode moves + from 0.4507 to 0.4284 Hz), and `09_iea15_umainesemi`, + `10_iea22_semi` and `11_upscale25_centraltower` now include the weight + restoring in `hydro_K` (roll / pitch +3.4e8, +1.3e7 and −1.4e9 + N·m/rad; all stay positive). Every sample README picks up the + corrected dense-solver frequencies. - **`Tower.add_point_mass` on a deck-built floater now updates the weight restoring.** The weight term is formed when the model is constructed, so a lump added afterwards reached the FEM mass matrix but diff --git a/VALIDATION.md b/VALIDATION.md index 2d36ce7..2f10b11 100644 --- a/VALIDATION.md +++ b/VALIDATION.md @@ -129,6 +129,7 @@ by title + publisher: - Jonkman, Butterfield, Musial & Scott (2009). *Definition of a 5-MW Reference Wind Turbine for Offshore System Development*. NREL/TP-500-38060. [PDF](https://docs.nrel.gov/docs/fy09osti/38060.pdf). - Jonkman & Musial (2010). *Offshore Code Comparison Collaboration (OC3) for IEA Wind Task 23*. NREL/TP-5000-48191. [PDF](https://docs.nrel.gov/docs/fy11osti/48191.pdf). - Jonkman (2010). *Definition of the Floating System for Phase IV of OC3*. NREL/TP-500-47535. [PDF](https://docs.nrel.gov/docs/fy10osti/47535.pdf). +- Robertson, Jonkman, Masciola, Song, Goupee, Coulling & Luan (2014). *Definition of the Semisubmersible Floating System for Phase II of OC4*. NREL/TP-5000-60601. [PDF](https://docs.nrel.gov/docs/fy14osti/60601.pdf). - Bortolotti, Tarrés, Dykes, Merz, Sethuraman, Verelst & Zahle (2019). *IEA Wind TCP Task 37: Systems Engineering in Wind Energy — WP2.1 Reference Wind Turbines*. NREL/TP-5000-73492. [PDF](https://docs.nrel.gov/docs/fy19osti/73492.pdf). - Gaertner et al. (2020). *IEA Wind TCP Task 37: Definition of the IEA 15-Megawatt Offshore Reference Wind Turbine*. NREL/TP-5000-75698. [PDF](https://docs.nrel.gov/docs/fy20osti/75698.pdf). - Allen et al. (2020). *Definition of the UMaine VolturnUS-S Reference Platform Developed for the IEA Wind 15-Megawatt Offshore Reference Wind Turbine*. NREL/TP-5000-76773. [PDF](https://docs.nrel.gov/docs/fy20osti/76773.pdf). @@ -230,6 +231,9 @@ is behavioural / contract-style. | Floating weight restoring $-g \sum m_i (z_i - z_{ref})$ in roll / pitch (platform + tower + RNA) | closed form; classical $C_{44} = \rho g V (z_B - z_G) + \rho g I_{wp}$ (Faltinsen 1990, ch. 2) | 6×6 weight matrix, reference-point shift, mass scaling, point masses | `rel = 1e-9` | exact | [`tests/test_floating_weight_restoring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_floating_weight_restoring.py) | no | | Negative platform restoring is reported, not silently dropped | construction (OC3 Hywind sample with the weight term removed) | `UserWarning` naming roll / pitch | always | (within tol) | [`tests/test_floating_weight_restoring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_floating_weight_restoring.py) | no | | **`Tower.from_elastodyn_with_mooring` — OC3 Hywind rigid-body frequencies** (WAMIT `.hst` + weight of platform, tower and RNA) | Jonkman (2010) NREL/TP-500-47535, FAST full-system natural frequencies: surge/sway 0.0080, heave 0.0324, roll 0.0342, pitch 0.0343 Hz | six labelled rigid-body modes | surge/sway/roll/pitch 3 %, heave 2 % | surge/sway 0.8 %, heave 0.0 %, roll 1.1 %, pitch 1.4 % (roll/pitch were missing before the weight term) | [`tests/test_floating_weight_restoring.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_floating_weight_restoring.py) | yes | +| HydroDyn filled-member ballast mass properties (body-fixed, as HydroDyn precomputes them) | closed form: solid cylinder, frustum volume / centroid / axial moment, parallel axis | mass, centre of mass, 3×3 inertia; partial / full / empty, horizontal, inclined, tapered members; both member-table layouts | `rel = 1e-12` | exact | [`tests/test_hydrodyn_ballast.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_hydrodyn_ballast.py) | no | +| **`Tower.from_elastodyn_with_mooring` — OC4 DeepCwind ballasted platform** (steel `PtfmMass` + two HydroDyn fill groups) | Robertson et al. (2014) NREL/TP-5000-60601: platform incl. ballast 1.3473e7 kg, CM 13.46 m below SWL | platform mass, CM depth; positive roll / pitch restoring; six labelled rigid-body modes; heave vs single-DOF estimate | mass `rel = 1e-3`, CM `abs = 0.01 m`, heave `rel = 2%` | mass +0.02 %, CM 0.000 m | [`tests/test_hydrodyn_ballast.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_hydrodyn_ballast.py) | yes | +| Bundled floating samples are statically stable | construction (samples 07-11) | six labelled rigid-body modes, no negative-restoring warning | always | (within tol) | [`tests/test_floating_samples_spectra.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_floating_samples_spectra.py) | no | | `MudlineFoundation.from_soil_properties` — Yu & Amdahl (2023) Table 9 DTU 10 MW flexible pile, E\_SO = 30 MPa | Yu and Amdahl (2023) Marine Structures 92, 103482, Table 9 row 2 | K\_hh, K\_hr, K\_rr via Shadlou and Bhattacharya (2016) formulas | `rel = 10%` | within 3.0 % on all three springs | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | | `MudlineFoundation.from_soil_properties` — Yu & Amdahl (2023) Table 9 DTU 10 MW rigid pile, E\_SO = 30 MPa | Yu and Amdahl (2023) Table 9 row 1 | K\_hh, K\_hr, K\_rr via Shadlou and Bhattacharya (2016) rigid-pile formulas | `rel = 5%` | within 2.9 % on all three springs | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | | `MudlineFoundation.as_mooring_K` — cross-coupling sign matches Jonkman 2010 OC3 convention | construction | K\[0,4\] < 0, K\[1,3\] > 0, K\[0,4\]·K\[1,3\] < 0 | always | (within tol) | [`tests/test_foundation.py`](https://github.com/SMI-Lab-Inha/pyBModes/blob/master/tests/test_foundation.py) | no | diff --git a/docs/conventions.rst b/docs/conventions.rst index af9345b..9ed759b 100644 --- a/docs/conventions.rst +++ b/docs/conventions.rst @@ -299,6 +299,12 @@ Common pitfalls a ballast-stabilised spar is negative in roll / pitch (OC3 Hywind: ``−5.0e9 N·m/rad``), and ``Tower.run`` warns that those rigid-body modes are unstable rather than returning them. +* **Leaving out ballast declared in HydroDyn.** Water ballast given as + HydroDyn filled-member groups is not in ElastoDyn's ``PtfmMass`` + (OC4 DeepCwind keeps all of it there). The deck constructors lump it + into ``mass_pform``, ``cm_pform`` and ``i_matrix`` as the rigid ballast + HydroDyn treats it as; a hand-built ``PlatformSupport`` has to include + it itself. * **Confusing OpenFAST's two frames.** ``PtfmCMzt`` (CM) and ``PtfmRefzt`` (reference) are both in the tower-base *t*-frame; diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/01_nrel5mw_land/README.md b/src/pybmodes/_examples/sample_inputs/reference_turbines/01_nrel5mw_land/README.md index 2c55b29..4cea766 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/01_nrel5mw_land/README.md +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/01_nrel5mw_land/README.md @@ -33,13 +33,13 @@ This sample is a cantilever (`hub_conn = 1`) clamped at TowerBsHt with the RNA l ### Tower -- 1st FA tower-bending: **0.3332 Hz** +- 1st FA tower-bending: **0.3304 Hz** ### Blade (spinning at deck `RotSpeed = 12.1 rpm`) -- 1st flap: **0.7343 Hz** -- 1st edge: **1.1083 Hz** -- 2nd flap: **2.0177 Hz** +- 1st flap: **0.7342 Hz** +- 1st edge: **1.1080 Hz** +- 2nd flap: **2.0175 Hz** ## Comparison with published values diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/02_nrel5mw_oc3monopile/README.md b/src/pybmodes/_examples/sample_inputs/reference_turbines/02_nrel5mw_oc3monopile/README.md index 4913a17..39b61dd 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/02_nrel5mw_oc3monopile/README.md +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/02_nrel5mw_oc3monopile/README.md @@ -33,13 +33,13 @@ This sample is a combined pile + tower cantilever (`hub_conn = 1`) clamped at th ### Tower -- 1st FA tower-bending: **0.2856 Hz** +- 1st FA tower-bending: **0.2869 Hz** ### Blade (spinning at deck `RotSpeed = 12.1 rpm`) -- 1st flap: **0.7343 Hz** -- 1st edge: **1.1083 Hz** -- 2nd flap: **2.0177 Hz** +- 1st flap: **0.7342 Hz** +- 1st edge: **1.1080 Hz** +- 2nd flap: **2.0175 Hz** ## Comparison with published values diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/03_iea34_land/README.md b/src/pybmodes/_examples/sample_inputs/reference_turbines/03_iea34_land/README.md index 56d4477..b904cb6 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/03_iea34_land/README.md +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/03_iea34_land/README.md @@ -33,13 +33,13 @@ This sample is a cantilever (`hub_conn = 1`) clamped at TowerBsHt with the RNA l ### Tower -- 1st FA tower-bending: **0.4172 Hz** +- 1st FA tower-bending: **0.4179 Hz** ### Blade (spinning at deck `RotSpeed = 5 rpm`) -- 1st flap: **0.7958 Hz** -- 1st edge: **1.0294 Hz** -- 2nd flap: **2.3593 Hz** +- 1st flap: **0.7960 Hz** +- 1st edge: **1.0296 Hz** +- 2nd flap: **2.3592 Hz** ## Comparison with published values diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/04_iea10_monopile/README.md b/src/pybmodes/_examples/sample_inputs/reference_turbines/04_iea10_monopile/README.md index fd6b25f..2fc5687 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/04_iea10_monopile/README.md +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/04_iea10_monopile/README.md @@ -33,13 +33,13 @@ This sample is a combined pile + tower cantilever (`hub_conn = 1`) clamped at th ### Tower -- 1st FA tower-bending: **0.3008 Hz** +- 1st FA tower-bending: **0.3000 Hz** ### Blade (spinning at deck `RotSpeed = 5 rpm`) -- 1st flap: **0.4044 Hz** -- 1st edge: **0.6496 Hz** -- 2nd flap: **1.0430 Hz** +- 1st flap: **0.4050 Hz** +- 1st edge: **0.6494 Hz** +- 2nd flap: **1.0433 Hz** ## Comparison with published values diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/05_iea15_monopile/README.md b/src/pybmodes/_examples/sample_inputs/reference_turbines/05_iea15_monopile/README.md index ec97393..6ff5815 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/05_iea15_monopile/README.md +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/05_iea15_monopile/README.md @@ -33,13 +33,13 @@ This sample is a combined pile + tower cantilever (`hub_conn = 1`) clamped at th ### Tower -- 1st FA tower-bending: **0.1906 Hz** +- 1st FA tower-bending: **0.1883 Hz** ### Blade (spinning at deck `RotSpeed = 7.55 rpm`) -- 1st flap: **0.5583 Hz** -- 1st edge: **0.7278 Hz** -- 2nd flap: **1.5915 Hz** +- 1st flap: **0.5585 Hz** +- 1st edge: **0.7281 Hz** +- 2nd flap: **1.5916 Hz** ## Comparison with published values diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/06_iea22_monopile/README.md b/src/pybmodes/_examples/sample_inputs/reference_turbines/06_iea22_monopile/README.md index 552902f..bab2531 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/06_iea22_monopile/README.md +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/06_iea22_monopile/README.md @@ -33,13 +33,13 @@ This sample is a combined pile + tower cantilever (`hub_conn = 1`) clamped at th ### Tower -- 1st FA tower-bending: **0.1819 Hz** +- 1st FA tower-bending: **0.1701 Hz** ### Blade (spinning at deck `RotSpeed = 3 rpm`) -- 1st flap: **0.3998 Hz** -- 1st edge: **0.5381 Hz** -- 2nd flap: **1.1032 Hz** +- 1st flap: **0.3987 Hz** +- 1st edge: **0.5388 Hz** +- 2nd flap: **1.1031 Hz** ## Comparison with published values diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/07_nrel5mw_oc3hywind_spar/README.md b/src/pybmodes/_examples/sample_inputs/reference_turbines/07_nrel5mw_oc3hywind_spar/README.md index d7ccda5..bcf5a30 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/07_nrel5mw_oc3hywind_spar/README.md +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/07_nrel5mw_oc3hywind_spar/README.md @@ -37,9 +37,9 @@ This sample is a free-free flexible-tower (`hub_conn = 2`) with a full PlatformS ### Blade (spinning at deck `RotSpeed = 12.1 rpm`) -- 1st flap: **0.7343 Hz** -- 1st edge: **1.1083 Hz** -- 2nd flap: **2.0177 Hz** +- 1st flap: **0.7342 Hz** +- 1st edge: **1.1080 Hz** +- 2nd flap: **2.0175 Hz** ## Comparison with published values diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/08_nrel5mw_oc4semi/08_nrel5mw_oc4semi_tower.bmi b/src/pybmodes/_examples/sample_inputs/reference_turbines/08_nrel5mw_oc4semi/08_nrel5mw_oc4semi_tower.bmi index be93077..7e87403 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/08_nrel5mw_oc4semi/08_nrel5mw_oc4semi_tower.bmi +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/08_nrel5mw_oc4semi/08_nrel5mw_oc4semi_tower.bmi @@ -50,12 +50,12 @@ Distance of element boundary nodes from tower base (norm. by tower length), el_l --------- Properties of tower support subsystem --------------------------------------------- 1 tow_support: 1 = floating-platform with optional tension wires (-) -10.000000 draft : depth of tower base from MSL (m) -8.658800 cm_pform : distance of platform c.m. below MSL (m) -3.852180e+06 mass_pform : platform mass (kg) +13.460321 0.000002 0.000000 cm_pform_xyz : platform c.m. below MSL + horizontal (x downwind, y lateral) offsets from the tower axis (m) +1.347603e+07 mass_pform : platform mass (kg) Platform mass inertia 3X3 matrix (i_matrix_pform): [OC4 DeepCwind semi (Robertson 2014) via HydroDyn+MoorDyn] - 2.561930e+09 0.000000e+00 0.000000e+00 - 0.000000e+00 2.561930e+09 0.000000e+00 - 0.000000e+00 0.000000e+00 4.242650e+09 + 7.099696e+09 0.000000e+00 6.165441e+01 + 0.000000e+00 7.099696e+09 0.000000e+00 + 6.165441e+01 0.000000e+00 1.280537e+10 -0.000000 ref_msl : distance of platform reference point below MSL (m) Hydrodynamic 6X6 mass matrix at platform ref point (hydro_M) [OC4 DeepCwind semi (Robertson 2014) via HydroDyn+MoorDyn]: 8.740059e+06 0.000000e+00 -3.871260e-01 0.000000e+00 -1.075816e+08 0.000000e+00 @@ -68,16 +68,16 @@ Hydrodynamic 6X6 stiffness matrix at platform ref point (hydro_K) [OC4 DeepCwin 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 3.820308e+06 0.000000e+00 -1.541053e+02 0.000000e+00 - 0.000000e+00 0.000000e+00 0.000000e+00 -3.807176e+08 0.000000e+00 0.000000e+00 - 0.000000e+00 0.000000e+00 -1.541053e+02 0.000000e+00 -3.807153e+08 0.000000e+00 + 0.000000e+00 0.000000e+00 0.000000e+00 9.852074e+08 0.000000e+00 0.000000e+00 + 0.000000e+00 0.000000e+00 -1.541053e+02 0.000000e+00 9.852097e+08 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 Mooring-system 6X6 stiffness matrix (mooring_K): [OC4 DeepCwind semi (Robertson 2014) via HydroDyn+MoorDyn] 7.083634e+04 -0.000000e+00 1.049534e-01 -0.000000e+00 -1.095301e+05 -0.000000e+00 -0.000000e+00 7.083592e+04 -0.000000e+00 1.095165e+05 -0.000000e+00 -1.105049e+01 - 1.049534e-01 -0.000000e+00 1.914043e+04 -0.000000e+00 -2.706354e+02 -5.820766e-10 + 1.049534e-01 -0.000000e+00 1.914043e+04 -0.000000e+00 -2.706354e+02 -0.000000e+00 -0.000000e+00 1.095165e+05 -0.000000e+00 8.770363e+07 -0.000000e+00 5.655611e+04 - -1.095301e+05 -0.000000e+00 -2.706354e+02 -0.000000e+00 8.770406e+07 -2.328306e-09 - -0.000000e+00 -1.105049e+01 -5.820766e-10 5.655611e+04 -2.328306e-09 1.180255e+08 + -1.095301e+05 -0.000000e+00 -2.706354e+02 -0.000000e+00 8.770406e+07 2.328306e-09 + -0.000000e+00 -1.105049e+01 -0.000000e+00 5.655611e+04 2.328306e-09 1.180255e+08 Distributed (hydrodynamic) added-mass per unit length along the flexible tower: 0 n_secs_m_distr: number of distributed-added-mass stations (-) diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/08_nrel5mw_oc4semi/README.md b/src/pybmodes/_examples/sample_inputs/reference_turbines/08_nrel5mw_oc4semi/README.md index 63936e2..56e6697 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/08_nrel5mw_oc4semi/README.md +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/08_nrel5mw_oc4semi/README.md @@ -33,13 +33,13 @@ This sample is a free-free flexible-tower (`hub_conn = 2`) with a full PlatformS ### Tower -- 1st FA tower-bending: **0.4507 Hz** +- 1st FA tower-bending: **0.4284 Hz** ### Blade (spinning at deck `RotSpeed = 9 rpm`) -- 1st flap: **0.7118 Hz** -- 1st edge: **1.1043 Hz** -- 2nd flap: **1.9913 Hz** +- 1st flap: **0.7117 Hz** +- 1st edge: **1.1042 Hz** +- 2nd flap: **1.9914 Hz** ## Comparison with published values diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/09_iea15_umainesemi/09_iea15_umainesemi_tower.bmi b/src/pybmodes/_examples/sample_inputs/reference_turbines/09_iea15_umainesemi/09_iea15_umainesemi_tower.bmi index 4a45a70..570c086 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/09_iea15_umainesemi/09_iea15_umainesemi_tower.bmi +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/09_iea15_umainesemi/09_iea15_umainesemi_tower.bmi @@ -68,8 +68,8 @@ Hydrodynamic 6X6 stiffness matrix at platform ref point (hydro_K) [UMaine Voltu 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 4.453443e+06 0.000000e+00 -4.033086e+03 0.000000e+00 - 0.000000e+00 0.000000e+00 0.000000e+00 2.193713e+09 0.000000e+00 -9.130073e+04 - 0.000000e+00 0.000000e+00 -4.033086e+03 0.000000e+00 2.193473e+09 0.000000e+00 + 0.000000e+00 0.000000e+00 0.000000e+00 2.538614e+09 0.000000e+00 -9.130073e+04 + 0.000000e+00 0.000000e+00 -4.033086e+03 0.000000e+00 2.538374e+09 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 -9.130073e+04 0.000000e+00 0.000000e+00 Mooring-system 6X6 stiffness matrix (mooring_K): [UMaine VolturnUS-S semi (Allen 2020) via HydroDyn+MoorDyn] 7.189209e+04 -0.000000e+00 -3.420710e-01 -0.000000e+00 1.146452e+06 -0.000000e+00 diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/09_iea15_umainesemi/README.md b/src/pybmodes/_examples/sample_inputs/reference_turbines/09_iea15_umainesemi/README.md index 221b94d..316e921 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/09_iea15_umainesemi/README.md +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/09_iea15_umainesemi/README.md @@ -33,13 +33,13 @@ This sample is a free-free flexible-tower (`hub_conn = 2`) with a full PlatformS ### Tower -- 1st FA tower-bending: **0.5255 Hz** +- 1st FA tower-bending: **0.5256 Hz** ### Blade (spinning at deck `RotSpeed = 7.55 rpm`) -- 1st flap: **0.5583 Hz** -- 1st edge: **0.7278 Hz** -- 2nd flap: **1.5915 Hz** +- 1st flap: **0.5585 Hz** +- 1st edge: **0.7281 Hz** +- 2nd flap: **1.5916 Hz** ## Comparison with published values diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/10_iea22_semi/10_iea22_semi_tower.bmi b/src/pybmodes/_examples/sample_inputs/reference_turbines/10_iea22_semi/10_iea22_semi_tower.bmi index b21af9d..f704925 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/10_iea22_semi/10_iea22_semi_tower.bmi +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/10_iea22_semi/10_iea22_semi_tower.bmi @@ -68,16 +68,16 @@ Hydrodynamic 6X6 stiffness matrix at platform ref point (hydro_K) [IEA-22 semi 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 4.837443e+06 0.000000e+00 -3.755001e+02 0.000000e+00 - 0.000000e+00 0.000000e+00 0.000000e+00 4.182874e+09 0.000000e+00 4.600509e+04 - 0.000000e+00 0.000000e+00 -3.755001e+02 0.000000e+00 4.182872e+09 0.000000e+00 + 0.000000e+00 0.000000e+00 0.000000e+00 4.196176e+09 0.000000e+00 4.600509e+04 + 0.000000e+00 0.000000e+00 -3.755001e+02 0.000000e+00 4.196174e+09 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 4.600509e+04 0.000000e+00 0.000000e+00 Mooring-system 6X6 stiffness matrix (mooring_K): [IEA-22 semi (Bortolotti 2024) via HydroDyn+MoorDyn] - 1.287829e+05 -0.000000e+00 6.722716e-02 -0.000000e+00 1.740192e+06 -0.000000e+00 + 1.287829e+05 -0.000000e+00 6.722721e-02 -0.000000e+00 1.740192e+06 -0.000000e+00 -0.000000e+00 1.287828e+05 -0.000000e+00 -1.740195e+06 -0.000000e+00 -3.233520e+00 - 6.722716e-02 -0.000000e+00 9.519838e+04 -0.000000e+00 8.539450e+02 -0.000000e+00 - -0.000000e+00 -1.740195e+06 -0.000000e+00 4.882500e+08 -0.000000e+00 2.892178e+04 + 6.722721e-02 -0.000000e+00 9.519838e+04 -4.656613e-09 8.539450e+02 -4.656613e-09 + -0.000000e+00 -1.740195e+06 -4.656613e-09 4.882500e+08 -0.000000e+00 2.892178e+04 1.740192e+06 -0.000000e+00 8.539450e+02 -0.000000e+00 4.882518e+08 -0.000000e+00 - -0.000000e+00 -3.233520e+00 -0.000000e+00 2.892178e+04 -0.000000e+00 4.879339e+08 + -0.000000e+00 -3.233520e+00 -4.656613e-09 2.892178e+04 -0.000000e+00 4.879339e+08 Distributed (hydrodynamic) added-mass per unit length along the flexible tower: 0 n_secs_m_distr: number of distributed-added-mass stations (-) diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/10_iea22_semi/README.md b/src/pybmodes/_examples/sample_inputs/reference_turbines/10_iea22_semi/README.md index 285b2d7..155c4fc 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/10_iea22_semi/README.md +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/10_iea22_semi/README.md @@ -37,9 +37,9 @@ This sample is a free-free flexible-tower (`hub_conn = 2`) with a full PlatformS ### Blade (spinning at deck `RotSpeed = 3 rpm`) -- 1st flap: **0.3998 Hz** -- 1st edge: **0.5381 Hz** -- 2nd flap: **1.1032 Hz** +- 1st flap: **0.3987 Hz** +- 1st edge: **0.5388 Hz** +- 2nd flap: **1.1031 Hz** ## Comparison with published values diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/11_upscale25_centraltower/11_upscale25_centraltower_tower.bmi b/src/pybmodes/_examples/sample_inputs/reference_turbines/11_upscale25_centraltower/11_upscale25_centraltower_tower.bmi index 68c6439..70ccbb9 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/11_upscale25_centraltower/11_upscale25_centraltower_tower.bmi +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/11_upscale25_centraltower/11_upscale25_centraltower_tower.bmi @@ -68,16 +68,16 @@ Hydrodynamic 6X6 stiffness matrix at platform ref point (hydro_K) [UPSCALE 25MW 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 6.667346e+06 0.000000e+00 -1.385558e+04 0.000000e+00 - 0.000000e+00 0.000000e+00 0.000000e+00 4.870685e+09 0.000000e+00 0.000000e+00 - 0.000000e+00 0.000000e+00 -1.385558e+04 0.000000e+00 4.871378e+09 0.000000e+00 + 0.000000e+00 0.000000e+00 0.000000e+00 3.466707e+09 0.000000e+00 0.000000e+00 + 0.000000e+00 0.000000e+00 -1.385558e+04 0.000000e+00 3.467399e+09 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 Mooring-system 6X6 stiffness matrix (mooring_K): [UPSCALE 25MW (Sandua-Fernández 2023) via HydroDyn+MoorDyn] - 1.168959e+05 -0.000000e+00 2.653586e-02 -2.910383e-10 2.410390e+06 -0.000000e+00 - -0.000000e+00 1.168959e+05 -0.000000e+00 -2.410386e+06 -0.000000e+00 4.525265e-01 - 2.653586e-02 -0.000000e+00 1.004196e+05 -4.656613e-09 1.700183e+03 -0.000000e+00 - -2.910383e-10 -2.410386e+06 -4.656613e-09 4.943385e+08 -0.000000e+00 -7.852156e+04 + 1.168959e+05 -0.000000e+00 2.653586e-02 -0.000000e+00 2.410390e+06 -0.000000e+00 + -0.000000e+00 1.168959e+05 -0.000000e+00 -2.410386e+06 -0.000000e+00 4.525264e-01 + 2.653586e-02 -0.000000e+00 1.004196e+05 4.656613e-09 1.700183e+03 -0.000000e+00 + -0.000000e+00 -2.410386e+06 4.656613e-09 4.943385e+08 -0.000000e+00 -7.852156e+04 2.410390e+06 -0.000000e+00 1.700183e+03 -0.000000e+00 4.943404e+08 -0.000000e+00 - -0.000000e+00 4.525265e-01 -0.000000e+00 -7.852156e+04 -0.000000e+00 4.823907e+08 + -0.000000e+00 4.525264e-01 -0.000000e+00 -7.852156e+04 -0.000000e+00 4.823907e+08 Distributed (hydrodynamic) added-mass per unit length along the flexible tower: 0 n_secs_m_distr: number of distributed-added-mass stations (-) diff --git a/src/pybmodes/_examples/sample_inputs/reference_turbines/11_upscale25_centraltower/README.md b/src/pybmodes/_examples/sample_inputs/reference_turbines/11_upscale25_centraltower/README.md index 8885d18..87edb69 100644 --- a/src/pybmodes/_examples/sample_inputs/reference_turbines/11_upscale25_centraltower/README.md +++ b/src/pybmodes/_examples/sample_inputs/reference_turbines/11_upscale25_centraltower/README.md @@ -33,13 +33,13 @@ This sample is a free-free flexible-tower (`hub_conn = 2`) with a full PlatformS ### Tower -- 1st FA tower-bending: **0.4372 Hz** +- 1st FA tower-bending: **0.4367 Hz** ### Blade (spinning at deck `RotSpeed = 6 rpm`) -- 1st flap: **0.3697 Hz** -- 1st edge: **0.4949 Hz** -- 2nd flap: **1.0694 Hz** +- 1st flap: **0.3706 Hz** +- 1st edge: **0.4955 Hz** +- 2nd flap: **1.0693 Hz** ## Comparison with published values diff --git a/src/pybmodes/models/_pipeline.py b/src/pybmodes/models/_pipeline.py index 58cbc8f..ef23ca2 100644 --- a/src/pybmodes/models/_pipeline.py +++ b/src/pybmodes/models/_pipeline.py @@ -253,9 +253,10 @@ def _warn_on_negative_platform_restoring(plat: PlatformSupport) -> None: f"eigenvalues, which the general (non-symmetric) solver path " f"leaves out of the returned modes. A WAMIT " f".hst excludes the body weight (the -m g z_G term), and ballast " - f"carried outside PtfmMass — e.g. HydroDyn filled members — is not " - f"in the model either; check that the restoring includes every " - f"weight and buoyancy contribution.", + f"carried outside the platform mass is easy to miss in a hand-built " + f"model (the HydroDyn-deck constructors add filled-member ballast " + f"themselves); check that the restoring includes every weight and " + f"buoyancy contribution.", UserWarning, stacklevel=4, # past run_fem and Tower.run, to the caller ) diff --git a/tests/test_floating_samples_spectra.py b/tests/test_floating_samples_spectra.py index a56b7dd..d639a5b 100644 --- a/tests/test_floating_samples_spectra.py +++ b/tests/test_floating_samples_spectra.py @@ -34,6 +34,7 @@ from __future__ import annotations import pathlib +import warnings import numpy as np import pytest @@ -133,3 +134,22 @@ def test_oc3hywind_sample_matches_bmodes_jj() -> None: f"(allowed 0.5%).\n pyBmodes : {np.array2string(f[:9], precision=5)}" f"\n BModes JJ: {np.array2string(_OC3_BMODES_JJ, precision=5)}" ) + + +@pytest.mark.parametrize("sample_id", _FLOATING) +def test_floating_sample_is_statically_stable(sample_id: str) -> None: + """Every bundled floater has positive roll / pitch restoring and returns + all six rigid-body modes without the negative-stiffness warning. + + The OC4 DeepCwind sample used to fail this: its water ballast sits in + HydroDyn filled members rather than ElastoDyn's PtfmMass, and until the + deck constructors lumped that ballast into the platform the sample + carried 3.85e6 kg of steel instead of 1.35e7 kg of platform, read + C44 = C55 = -3.8e8 N·m/rad, and labelled tower-bending modes as roll + and pitch.""" + with warnings.catch_warnings(): + warnings.simplefilter("error", UserWarning) + res = Tower(_tower_bmi(sample_id)).run(n_modes=9, check_model=False) + assert set(res.mode_labels[:6]) == { + "surge", "sway", "heave", "roll", "pitch", "yaw", + }, (sample_id, res.mode_labels[:6]) diff --git a/tests/test_hydrodyn_ballast.py b/tests/test_hydrodyn_ballast.py index a90e41a..0d07d7e 100644 --- a/tests/test_hydrodyn_ballast.py +++ b/tests/test_hydrodyn_ballast.py @@ -23,6 +23,7 @@ import numpy as np import pytest +from pybmodes._numeric import trapezoid from pybmodes.io._hydrodyn_ballast import read_filled_ballast from pybmodes.models._platform import _add_filled_ballast @@ -379,8 +380,7 @@ def test_heave_matches_the_single_dof_estimate(self): ps = tower._bmi.support res = tower.run(n_modes=8, check_model=False) f = dict(zip(res.mode_labels, res.frequencies))["heave"] - m_tower = float(np.trapezoid( - tower._sp.mass_den, tower._sp.span_loc)) * ( + m_tower = trapezoid(tower._sp.mass_den, tower._sp.span_loc) * ( tower._bmi.radius + ps.draft) m_total = ps.mass_pform + m_tower + tower._bmi.tip_mass.mass k33 = ps.hydro_K[2, 2] + ps.mooring_K[2, 2] From a7b4ed5d1f2746e532f4c48c1a003e9e2aebf11b Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 30 Sep 2026 19:26:56 +0900 Subject: [PATCH 3/4] fix: resolve DEFAULT fill density with the WindIO path's rho --- src/pybmodes/models/_platform.py | 9 ++++++++- src/pybmodes/models/tower.py | 8 ++++++-- tests/test_hydrodyn_ballast.py | 13 +++++++++++++ 3 files changed, 27 insertions(+), 3 deletions(-) diff --git a/src/pybmodes/models/_platform.py b/src/pybmodes/models/_platform.py index e028cd5..650ac68 100644 --- a/src/pybmodes/models/_platform.py +++ b/src/pybmodes/models/_platform.py @@ -146,6 +146,8 @@ def _add_filled_ballast( ptfm: dict[str, float], i_mat: np.ndarray, hydrodyn_dat_path: pathlib.Path | str, + *, + water_density: float | None = None, ) -> tuple[dict[str, float], np.ndarray]: """Lump the HydroDyn filled-member ballast into the platform. @@ -163,12 +165,17 @@ def _add_filled_ballast( :mod:`pybmodes.io._hydrodyn_ballast`), which is exactly a rigid mass on the platform. The combined rotational block is a full 3×3 tensor; for a symmetric platform its products of inertia vanish. + + ``water_density`` resolves ``FillDens = DEFAULT``; when omitted the + deck's own ``WtrDens`` (or 1025 kg/m³) is used. """ import numpy as np from pybmodes.io._hydrodyn_ballast import read_filled_ballast - ballast = read_filled_ballast(hydrodyn_dat_path) + ballast = read_filled_ballast( + hydrodyn_dat_path, water_density=water_density, + ) if ballast is None: return ptfm, i_mat diff --git a/src/pybmodes/models/tower.py b/src/pybmodes/models/tower.py index 43ae806..a517f35 100644 --- a/src/pybmodes/models/tower.py +++ b/src/pybmodes/models/tower.py @@ -1403,9 +1403,13 @@ def _tower_grid_od_wt( rna_tip = _tower_top_assembly_mass(main, blade) ptfm = _scan_platform_fields(ed_path) i_mat = _platform_inertia_matrix(ptfm) - # Filled-member ballast, as in from_elastodyn_with_mooring. + # Filled-member ballast, as in from_elastodyn_with_mooring, with + # DEFAULT fill density resolved to the same rho the hydrostatics + # and mooring on this path use. if hydrodyn_dat is not None: - ptfm, i_mat = _add_filled_ballast(ptfm, i_mat, hydrodyn_dat) + ptfm, i_mat = _add_filled_ballast( + ptfm, i_mat, hydrodyn_dat, water_density=rho, + ) M = ptfm["PtfmMass"] cm_pform = -ptfm["PtfmCMzt"] cm_x, cm_y = ptfm["PtfmCMxt"], ptfm["PtfmCMyt"] diff --git a/tests/test_hydrodyn_ballast.py b/tests/test_hydrodyn_ballast.py index 0d07d7e..e307126 100644 --- a/tests/test_hydrodyn_ballast.py +++ b/tests/test_hydrodyn_ballast.py @@ -328,6 +328,19 @@ def shift(mass, d): # The input dict is not mutated. assert p["PtfmMass"] == 4.0e6 + def test_water_density_resolves_default_fill(self, tmp_path): + """from_windio_floating passes its rho so the ballast and the + hydrostatics agree on the water density.""" + path = _deck( + tmp_path, joints=[(1, 0, 0, -10), (2, 0, 0, 0)], + props=[(1, 2.0, 0.0)], members=[(1, 1, 2, 1, 1)], + fills=["1 1 0.0 DEFAULT"], + ) + p = self._ptfm() + out, _ = _add_filled_ballast( + p, self._i_mat(p), path, water_density=1000.0) + assert out["PtfmMass"] == pytest.approx(4.0e6 + 1000.0 * np.pi * 10.0) + def test_no_fill_leaves_platform_unchanged(self, tmp_path): path = _deck( tmp_path, joints=[(1, 0, 0, -1), (2, 0, 0, 1)], From ab0f561e448c011d3101f7431aa31cda88e165bb Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 30 Sep 2026 19:36:19 +0900 Subject: [PATCH 4/4] fix: reject non-finite HydroDyn ballast scalars explicitly --- src/pybmodes/io/_hydrodyn_ballast.py | 11 ++++++++++- tests/test_hydrodyn_ballast.py | 21 +++++++++++++++++++++ 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/src/pybmodes/io/_hydrodyn_ballast.py b/src/pybmodes/io/_hydrodyn_ballast.py index 7a5b758..342c2f0 100644 --- a/src/pybmodes/io/_hydrodyn_ballast.py +++ b/src/pybmodes/io/_hydrodyn_ballast.py @@ -141,10 +141,19 @@ def _table( def _float(tok: str, what: str, path: pathlib.Path) -> float: + """Parse one deck scalar, rejecting anything non-finite. + + A NaN would slip past every later range check (comparisons with NaN + are false) and reach the platform mass matrix, so finiteness is + checked here rather than left to the shared parser. + """ try: - return _parse_fortran_float(tok) + value = _parse_fortran_float(tok) except ValueError as err: raise ValueError(f"Malformed {what} in {path}: {tok!r}") from err + if not np.isfinite(value): + raise ValueError(f"Non-finite {what} in {path}: {tok!r}") + return value def _int(tok: str, what: str, path: pathlib.Path) -> int: diff --git a/tests/test_hydrodyn_ballast.py b/tests/test_hydrodyn_ballast.py index e307126..89238da 100644 --- a/tests/test_hydrodyn_ballast.py +++ b/tests/test_hydrodyn_ballast.py @@ -280,6 +280,27 @@ def test_non_physical_section_is_rejected(self, tmp_path): with pytest.raises(ValueError, match="non-physical"): read_filled_ballast(path) + @pytest.mark.parametrize("token", ["nan", "NaN", "inf", "-Infinity"]) + @pytest.mark.parametrize("where", ["FillDens", "FillFSLoc", "joint", "PropD"]) + def test_non_finite_deck_scalars_are_rejected(self, tmp_path, token, where): + joints = [(1, 0, 0, -10), (2, 0, 0, 0)] + props = [(1, 2.0, 0.0)] + fs, dens = "0.0", f"{RHO}" + if where == "FillDens": + dens = token + elif where == "FillFSLoc": + fs = token + elif where == "joint": + joints = [(1, 0, token, -10), (2, 0, 0, 0)] + else: + props = [(1, token, 0.0)] + path = _deck( + tmp_path, joints=joints, props=props, + members=[(1, 1, 2, 1, 1)], fills=[f"1 1 {fs} {dens}"], + ) + with pytest.raises(ValueError, match="(?i)non-finite|malformed"): + read_filled_ballast(path) + @pytest.mark.parametrize("bad", [0.0, -1.0, float("nan"), float("inf")]) def test_bad_water_density_is_rejected(self, tmp_path, bad): path = _vertical_column(tmp_path, fs=-8.0)