API reference

Public API

Op3: OptumGX-OpenSeesPy-OpenFAST integration framework for offshore wind.

Public API:

from op3 import build_foundation, compose_tower_model
foundation = build_foundation(mode='fixed')
model = compose_tower_model(
    rotor='nrel_5mw_baseline',
    tower='nrel_5mw_tower',
    foundation=foundation,
)
freqs = model.eigen(n_modes=3)

class op3.FoundationMode(*values)

Bases: str, Enum

The five Op^3 foundation representations (LEGACY, frozen at v1.0).

The first four are the original v1.0 fidelity ladder (Modes A-D). distributed_bnwf_nonlinear (v1.1+) is the blueprint Q1(a) primary: a physically-distributed skirt column with per-depth PySimple1/ TzSimple1 backbones.

Deprecation notice (v1.1+): these values represent SSI fidelity levels, not foundation types. New code should use the type/SSI split under op3.foundations.types and op3.ssi. The legacy build_foundation(mode=...) path remains functional throughout v1.x.

DISSIPATION_WEIGHTED = 'dissipation_weighted'

DISTRIBUTED_BNWF = 'distributed_bnwf'

DISTRIBUTED_BNWF_NONLINEAR = 'distributed_bnwf_nonlinear'

FIXED = 'fixed'

STIFFNESS_6X6 = 'stiffness_6x6'

op3.build_foundation(mode: str, *, spring_profile: str | Path | None = None, stiffness_matrix: str | Path | ndarray | None = None, ogx_dissipation: str | Path | None = None, ogx_capacity: str | Path | None = None, scour_depth: float = 0.0, mode_d_alpha: float = 1.0, mode_d_beta: float = 0.05, diameter_m: float = 8.0, skirt_length_m: float | None = None, skirt_thickness_m: float = 0.025, physical: bool = False, _suppress_deprecation_warning: bool = False) → Foundation

Construct a Foundation handle ready for the composer (LEGACY API).

Deprecated (v1.1+): use op3.foundations.types and op3.ssi for new code. Example:

>>> from op3.foundations.types import Monopile
>>> from op3.ssi import Stiffness6x6
>>> mono = Monopile.from_oc3_spec(soil_profile=...)
>>> mono.with_ssi(Stiffness6x6(K=mono.head_stiffness_6x6()))

Parameters:

• mode (str) – One of ‘fixed’, ‘stiffness_6x6’, ‘distributed_bnwf’, ‘dissipation_weighted’, ‘distributed_bnwf_nonlinear’. Case-insensitive.

• spring_profile (str or Path, optional) – Path to a CSV with columns (depth_m, k_ini_kN_per_m, p_ult_kN_per_m, spring_type). Required for Modes C and D.

• stiffness_matrix (str, Path, or ndarray, optional) – A 6x6 stiffness matrix, either a path to a CSV or a NumPy array. Required for Mode B.

• ogx_dissipation (str or Path, optional) – Path to a CSV with columns (depth_m, w_z, D_total_kJ). Required for Mode D.

• ogx_capacity (str or Path, optional) – Path to the OptumGX power-law capacity CSV with columns (param_name, param_value). Used by Mode D for the ultimate resistance scaling.

• scour_depth (float, default 0.0) – Scour depth in meters. Affects Modes B, C, D.

• _suppress_deprecation_warning (bool) – Internal flag used by op3.foundations.types adapters to avoid spurious deprecation noise on the back-compat bridge.

op3.compose_tower_model(rotor: str, tower: str, foundation: Foundation, damping_ratio: float = 0.01) → TowerModel

Build a TowerModel from the given rotor, tower, and foundation.

Parameters:

• rotor (str) –

  Name of a rotor template. Valid values:

  ’nrel_5mw_baseline’, ‘iea_15mw_rwt’, ‘ref_4mw_owt’, ‘nrel_1.72_103’, ‘nrel_2.8_127’, ‘vestas_v27’

• tower (str) –

  Name of a tower template. Valid values:

  ’nrel_5mw_tower’, ‘iea_15mw_tower’, ‘site_a_rt1_tower’, ‘iea_land_onshore_tower’

• foundation (Foundation) – A Foundation handle from build_foundation().

• damping_ratio (float) – Rayleigh structural damping ratio (fraction of critical). Default 0.01 (1%) matches most NREL reference decks.

Returns:

Not yet built. Call model.eigen() or model.extract_6x6_stiffness() which will trigger .build() internally.

Return type:

TowerModel

op3.cross_compare(rotor: str, tower: str, scour_levels: Sequence[float], *, spring_profile: str | None = None, stiffness_matrix: str | None = None, ogx_dissipation: str | None = None, modes: Sequence[str] | None = None) → DataFrame

Run a cross-comparison of foundation modes and scour levels.

Parameters:

• rotor (str) – Rotor template name (see composer.compose_tower_model).

• tower (str) – Tower template name.

• scour_levels (sequence of float) – Scour depths in meters.

• spring_profile (str, optional) – Path to the distributed-BNWF spring CSV. Required for Modes C and D.

• stiffness_matrix (str, optional) – Path to the 6x6 stiffness matrix CSV. Required for Mode B.

• ogx_dissipation (str, optional) – Path to the dissipation profile CSV. Required for Mode D.

• modes (sequence of str, optional) – Which modes to include. Default: all four.

Returns:

One row per (mode, scour_m) with columns:

mode, scour_m, f1_Hz, f2_Hz, f3_Hz, integrated_K_kN_per_m, wall_clock_s

Return type:

pandas.DataFrame

op3.load_site_config(path: str | Path) → dict[str, Any]

Load and lightly validate a site configuration YAML.

Parameters:

path (str or Path) – Path to the YAML file. Relative paths are resolved from the repository root.

Returns:

Parsed YAML with a _source_path key added for provenance.

Return type:

dict

Raises:

• FileNotFoundError – If the YAML file does not exist.

• ValueError – If the YAML is missing a required top-level section.

Foundation modes

Op^3 foundations package.

Two coexisting APIs:

1. Legacy API (frozen at v1.0) — re-exported from op3.foundations._legacy. Emits DeprecationWarning when build_foundation() is called. Will be removed in v2.0.

2. New type/SSI API (v1.1+) — op3.foundations.types provides concrete foundation topologies (Monopile, Tripod, Jacket, SuctionBucket); op3.ssi provides SSI fidelity strategies (Fixed, Stiffness6x6, BNWFLumped, BNWFPhysical, CraigBampton). Each model under op3.models is a validated-dossier instance that combines the two.

Import cheat sheet

>>> # Legacy (still works; DeprecationWarning)
>>> from op3.foundations import build_foundation, FoundationMode
>>>
>>> # New types
>>> from op3.foundations.types import Monopile
>>> from op3.ssi import Stiffness6x6
>>>
>>> # Foundation protocol (for duck-typing / type hints)
>>> from op3.foundations.base import FoundationProtocol

The package import surface is intentionally thin: everything below is re-exported from either _legacy (with a deprecation trail) or base (the new protocol). Concrete types are NOT re-exported here; import them explicitly from op3.foundations.types.

class op3.foundations.BaseFoundation

Bases: ABC

Optional base class with a default back-compat bridge.

Concrete types may subclass BaseFoundation for the as_legacy_foundation() default (which delegates to head_stiffness_6x6()), or implement FoundationProtocol directly via duck typing.

Subclasses MUST:

• set type_name and foundation_type class attributes

• implement head_stiffness_6x6()

Subclasses SHOULD:

• provide a classmethod .from_dossier(path) that loads site.yaml, geometry.yaml, soil.yaml from a op3.models subdirectory and returns an instance

• carry a ssi attribute set by with_ssi() so the same topology can be re-used across fidelity levels

as_legacy_foundation() → Foundation

Default back-compat bridge: wrap head_stiffness_6x6() in a legacy Foundation(mode=STIFFNESS_6X6) object.

The _suppress_deprecation_warning kwarg on build_foundation() is intentionally NOT used here because Foundation(...) is constructed directly (no factory call). The legacy enum itself is silent until called.

foundation_type: FoundationType = 'monopile'

abstractmethod head_stiffness_6x6() → ndarray

Compute or return the 6x6 interface stiffness.

Concrete types typically delegate to self.ssi.compute_head_stiffness(self).

ssi: SSIProtocol | None = None

type_name: str = '<abstract>'

with_ssi(ssi: SSIProtocol) → BaseFoundation

Attach an SSI strategy and return self for chaining.

The strategy is stored on the instance; repeated calls replace the previous strategy. The same type instance can therefore be evaluated under multiple fidelity levels by calling with_ssi() then head_stiffness_6x6() in sequence.

class op3.foundations.Foundation(mode: FoundationMode, spring_table: DataFrame | None = None, stiffness_matrix: ndarray | None = None, dissipation_weights: DataFrame | None = None, mode_d_alpha: float = 1.0, mode_d_beta: float = 0.05, scour_depth: float = 0.0, diameter_m: float = 8.0, skirt_length_m: float | None = None, skirt_thickness_m: float = 0.025, physical: bool = False, base_H_stiffness_fraction: float | None = None, base_V_to_H_ratio: float | None = None, shaft_t_to_p_ratio: float | None = None, shaft_kz_to_kx_ratio: float | None = None, missing_pult_fallback_factor: float | None = None, source: str = 'analytical', diagnostics: dict = <factory>)

Bases: object

Opaque handle returned by build_foundation.

The composer calls attach_to_opensees(base_node) at model-build time. The foundation is therefore described declaratively here and only touches OpenSees tags when the composer passes it a base node.

attach_to_opensees(base_node: int) → dict

Instantiate the foundation as OpenSees elements and return a diagnostics dict.

This is called by the composer at model-build time. The concrete OpenSees commands live in the opensees_foundations submodule because they require an active OpenSees model.

Parameters:

base_node (int) – The OpenSees node tag at the tower base that the foundation attaches to.

Returns:

Diagnostic info: number of springs, integrated stiffness, energy balance check, etc.

Return type:

dict

base_H_stiffness_fraction: float | None = None

base_V_to_H_ratio: float | None = None

diagnostics: dict

diameter_m: float = 8.0

dissipation_weights: DataFrame | None = None

missing_pult_fallback_factor: float | None = None

mode: FoundationMode

mode_d_alpha: float = 1.0

mode_d_beta: float = 0.05

physical: bool = False

scour_depth: float = 0.0

shaft_kz_to_kx_ratio: float | None = None

shaft_t_to_p_ratio: float | None = None

skirt_length_m: float | None = None

skirt_thickness_m: float = 0.025

source: str = 'analytical'

spring_table: DataFrame | None = None

stiffness_matrix: ndarray | None = None

class op3.foundations.FoundationMode(*values)

Bases: str, Enum

The five Op^3 foundation representations (LEGACY, frozen at v1.0).

The first four are the original v1.0 fidelity ladder (Modes A-D). distributed_bnwf_nonlinear (v1.1+) is the blueprint Q1(a) primary: a physically-distributed skirt column with per-depth PySimple1/ TzSimple1 backbones.

Deprecation notice (v1.1+): these values represent SSI fidelity levels, not foundation types. New code should use the type/SSI split under op3.foundations.types and op3.ssi. The legacy build_foundation(mode=...) path remains functional throughout v1.x.

DISSIPATION_WEIGHTED = 'dissipation_weighted'

DISTRIBUTED_BNWF = 'distributed_bnwf'

DISTRIBUTED_BNWF_NONLINEAR = 'distributed_bnwf_nonlinear'

FIXED = 'fixed'

STIFFNESS_6X6 = 'stiffness_6x6'

class op3.foundations.FoundationProtocol(*args, **kwargs)

Bases: Protocol

Contract for a foundation TYPE.

Implementations own geometry + mass + the OpenSees topology they need to build, and compose with an SSIProtocol strategy for soil-structure interaction. The protocol is deliberately small so new types (e.g. gravity-based structures) can be added without touching the core pipeline.

as_legacy_foundation() → Foundation

Back-compat bridge: return a legacy Foundation (mode=STIFFNESS_6X6, stiffness_matrix=head_stiffness_6x6()) so the v1.0 op3.composer.compose_tower_model() pipeline works unchanged.

Future types with non-trivial topology (skirt, braces, etc.) will grow a second bridge, build_opensees(ops, base_node), that instantiates the real topology directly. For v1.1 the condensed-6x6 bridge is the only supported handover.

foundation_type: FoundationType

High-level classification (see FoundationType).

head_stiffness_6x6() → ndarray

Return the 6x6 stiffness at the tower-base (foundation-top) interface, in SI units (N/m, N·m/rad, mixed off-diagonals).

For fixed-base and elastic-head-spring SSI strategies this is analytical. For physical BNWF with Craig-Bampton reduction it may trigger an OpenSees model build and matrix condensation.

type_name: str

Canonical short name, e.g. "monopile" or "tripod". Used for logging, dossier lookup, and registry keys.

class op3.foundations.FoundationType(*values)

Bases: str, Enum

High-level classification of supported foundation topologies.

These are DESCRIPTIVE labels attached to each concrete type; they do NOT replace the legacy FoundationMode enum (which describes SSI fidelity, not topology). A single FoundationType value can be instantiated with any compatible SSIProtocol.

GBS = 'gbs'

JACKET = 'jacket'

MONOPILE = 'monopile'

SUCTION_BUCKET = 'suction_bucket'

TRIPOD = 'tripod'

op3.foundations.apply_scour_relief(spring_table: DataFrame, scour_depth: float) → DataFrame

Apply a stress-relief factor to a spring table for a given scour depth.

Rules:

• Nodes above the scoured mudline (z < scour_depth) have zero stiffness and zero capacity (the soil is gone).

• Nodes below but near the scour front have a smoothly tapered factor relief(z) = sqrt((z - scour) / z) to account for stress relief in the remaining soil column.

• Nodes far below the scour front are unchanged.

This is the stress-correction procedure from Chapter 6 of the dissertation and is mathematically identical to the factor documented in Appendix A Section A.3.

op3.foundations.build_foundation(mode: str, *, spring_profile: str | Path | None = None, stiffness_matrix: str | Path | ndarray | None = None, ogx_dissipation: str | Path | None = None, ogx_capacity: str | Path | None = None, scour_depth: float = 0.0, mode_d_alpha: float = 1.0, mode_d_beta: float = 0.05, diameter_m: float = 8.0, skirt_length_m: float | None = None, skirt_thickness_m: float = 0.025, physical: bool = False, _suppress_deprecation_warning: bool = False) → Foundation

Construct a Foundation handle ready for the composer (LEGACY API).

Deprecated (v1.1+): use op3.foundations.types and op3.ssi for new code. Example:

>>> from op3.foundations.types import Monopile
>>> from op3.ssi import Stiffness6x6
>>> mono = Monopile.from_oc3_spec(soil_profile=...)
>>> mono.with_ssi(Stiffness6x6(K=mono.head_stiffness_6x6()))

Parameters:

• mode (str) – One of ‘fixed’, ‘stiffness_6x6’, ‘distributed_bnwf’, ‘dissipation_weighted’, ‘distributed_bnwf_nonlinear’. Case-insensitive.

• spring_profile (str or Path, optional) – Path to a CSV with columns (depth_m, k_ini_kN_per_m, p_ult_kN_per_m, spring_type). Required for Modes C and D.

• stiffness_matrix (str, Path, or ndarray, optional) – A 6x6 stiffness matrix, either a path to a CSV or a NumPy array. Required for Mode B.

• ogx_dissipation (str or Path, optional) – Path to a CSV with columns (depth_m, w_z, D_total_kJ). Required for Mode D.

• ogx_capacity (str or Path, optional) – Path to the OptumGX power-law capacity CSV with columns (param_name, param_value). Used by Mode D for the ultimate resistance scaling.

• scour_depth (float, default 0.0) – Scour depth in meters. Affects Modes B, C, D.

• _suppress_deprecation_warning (bool) – Internal flag used by op3.foundations.types adapters to avoid spurious deprecation noise on the back-compat bridge.

op3.foundations.foundation_from_pisa(*, diameter_m: float, embed_length_m: float, soil_profile: list, n_segments: int = 50) → Foundation

Construct a STIFFNESS_6X6 Foundation whose K matrix is derived from the PISA framework (Burd 2020 / Byrne 2020). This is the canonical Op^3 entry point for monopile foundations when site-specific p-y curves are not available.

Parameters:

• diameter_m – Pile geometry.

• embed_length_m – Pile geometry.

• soil_profile (list[op3.standards.pisa.SoilState]) – Layered soil definition (depth, G, su or phi, soil_type).

• n_segments – Vertical discretisation for the PISA integration (default 50).

Composer

Op^3 tower + foundation composer.

The composer assembles an OpenSeesPy model from three orthogonal choices:

(rotor template) x (tower template) x (foundation module)

The rotor and tower templates live in op3.opensees_foundations.templates and map to published reference designs (NREL 5MW, IEA 15MW, Reference 4 MW OWT, etc.). The foundation module is one of the four Op^3 foundation modes from op3.foundations.

Example

>>> from op3 import build_foundation, compose_tower_model
>>> f = build_foundation(mode='distributed_bnwf',
...                      spring_profile='data/fem_results/opensees_spring_stiffness.csv',
...                      scour_depth=1.0)
>>> model = compose_tower_model(
...     rotor='nrel_5mw_baseline',
...     tower='nrel_5mw_tower',
...     foundation=f,
...     damping_ratio=0.01,
... )
>>> freqs = model.eigen(n_modes=6)
>>> print(freqs)

class op3.composer.TowerModel(rotor_name: str, tower_name: str, foundation: Foundation, _built: bool = False, _eigen_freqs: ndarray | None = None, _k_ssi_6x6: ndarray | None = None)

Bases: object

Handle returned by the composer.

Wraps an OpenSees domain and exposes the standard three analyses (eigenvalue, pushover, transient) plus a 6x6 static-condensation stiffness extraction for OpenFAST SubDyn export.

build() → None

Instantiate the OpenSees model. Idempotent.

eigen(n_modes: int = 6) → ndarray

Run an eigenvalue analysis and return the first n natural frequencies in Hz.

extract_6x6_stiffness() → ndarray

Static condensation of the foundation DOFs onto the tower base node. Returns a symmetric 6x6 stiffness matrix suitable for OpenFAST SubDyn.

foundation: Foundation

pushover(target_disp_m: float = 1.0, n_steps: int = 50) → dict

Run a static lateral pushover at the hub node.

Returns a dict with displacement_m and reaction_kN arrays.

rotor_name: str

tower_name: str

transient(duration_s: float = 10.0, dt_s: float = 0.01, damping_ratio: float = 0.01) → dict

Run a free-vibration transient with hub-node initial perturbation.

Returns a dict with time_s and hub_disp_m arrays for the hub node.

op3.composer.compose_tower_model(rotor: str, tower: str, foundation: Foundation, damping_ratio: float = 0.01) → TowerModel

Build a TowerModel from the given rotor, tower, and foundation.

Parameters:

• rotor (str) –

  Name of a rotor template. Valid values:

  ’nrel_5mw_baseline’, ‘iea_15mw_rwt’, ‘ref_4mw_owt’, ‘nrel_1.72_103’, ‘nrel_2.8_127’, ‘vestas_v27’

• tower (str) –

  Name of a tower template. Valid values:

  ’nrel_5mw_tower’, ‘iea_15mw_tower’, ‘site_a_rt1_tower’, ‘iea_land_onshore_tower’

• foundation (Foundation) – A Foundation handle from build_foundation().

• damping_ratio (float) – Rayleigh structural damping ratio (fraction of critical). Default 0.01 (1%) matches most NREL reference decks.

Returns:

Not yet built. Call model.eigen() or model.extract_6x6_stiffness() which will trigger .build() internally.

Return type:

TowerModel

Industry standards (Mode B)

DNVGL-ST-0126 Section 5.5 — Equivalent linear soil stiffness for offshore wind turbine support structures.

This module implements the closed-form expressions from DNVGL-ST-0126 (Edition 2016, amended 2021). The standard provides diagonal 6x6 stiffness matrices for the three main offshore wind foundation types: monopiles, jackets, and suction buckets.

Reference:

DNVGL-ST-0126 (2021). “Support structures for wind turbines”. DNV. https://www.dnv.com/oilgas/download/dnvgl-st-0126-support-structures-for-wind-turbines.html

The expressions implemented here are the equivalent linear values suitable for first-mode eigenvalue analysis. They are NOT intended for ultimate limit state design, fatigue assessment, or any nonlinear analysis. Users requiring those should run a full nonlinear soil-pile interaction analysis (e.g. PISA, p-y curves, or 3D finite element).

All formulas are documented inline with the page or section number of the standard. Where the standard refers to a coefficient table (e.g. for shape factors), the values are reproduced verbatim with attribution.

op3.standards.dnv_st_0126.dnv_jacket_stiffness(leg_spacing_m: float, n_legs: int = 4, pile_diameter_m: float = 1.5, pile_embedment_m: float = 30.0, soil_type: str = 'dense_sand', G: float | None = None, nu: float | None = None) → ndarray

Equivalent 6x6 stiffness for a jacket on n piles per DNVGL-ST-0126.

Combines individual pile stiffnesses with the jacket leg geometry to give an effective foundation stiffness at the tower base. The rotational stiffness of a multi-leg jacket is dominated by the pile spacing rather than the individual pile rotational stiffness.

Parameters:

• leg_spacing_m (float) – Distance from jacket centerline to each leg, in meters.

• n_legs (int, default 4) – Number of jacket legs (typically 3 or 4).

• pile_diameter_m (float, default 1.5) – Per-pile diameter in meters.

• pile_embedment_m (float, default 30.0) – Per-pile embedded length in meters.

• soil_type (as in dnv_monopile_stiffness.)

• G (as in dnv_monopile_stiffness.)

• nu (as in dnv_monopile_stiffness.)

• Reference

• ---------

• (2021) (DNVGL-ST-0126)

• structures" (Section 5.5.3 "Jacket support)

op3.standards.dnv_st_0126.dnv_monopile_stiffness(diameter_m: float, embedment_m: float, soil_type: str = 'dense_sand', G: float | None = None, nu: float | None = None) → ndarray

Equivalent linear 6x6 stiffness for a monopile per DNVGL-ST-0126 §5.5.

Uses the cylindrical pile in elastic half-space approximation (Randolph 1981, adopted by DNV) with depth-correction factor for embedment > 1 diameter.

Parameters:

• diameter_m (float) – Pile outer diameter in meters.

• embedment_m (float) – Embedded length in meters (mudline to pile tip).

• soil_type (str, default 'dense_sand') – One of the keys in DEFAULT_SHEAR_MODULUS.

• G (float, optional) – Soil shear modulus in Pa. Overrides soil_type default.

• nu (float, optional) – Poisson’s ratio. Overrides soil_type default.

Returns:

• numpy.ndarray – 6x6 diagonal stiffness matrix in SI units (N/m for translation, N*m/rad for rotation).

• Reference

• ———

• DNVGL-ST-0126 (2021), Section 5.5.2 “Linear soil-pile stiffness”

• Randolph, M.F. (1981). “The response of flexible piles to lateral – loading”. Geotechnique 31(2), 247-259.

op3.standards.dnv_st_0126.dnv_suction_bucket_stiffness(diameter_m: float, skirt_length_m: float, n_buckets: int = 1, spacing_m: float | None = None, soil_type: str = 'soft_clay', G: float | None = None, nu: float | None = None) → ndarray

Equivalent 6x6 stiffness for one or more suction bucket caissons.

For a single bucket (n_buckets = 1), the formula is the embedded cylindrical foundation per DNVGL-ST-0126 Annex F (which references Gazetas 1991). For multi-bucket configurations (tripod n=3, quadrupod n=4), the buckets are treated as independent feet at spacing_m from the centerline and combined like a jacket.

Parameters:

• diameter_m (float) – Bucket outer diameter in meters.

• skirt_length_m (float) – Skirt embedment depth in meters.

• n_buckets (int, default 1) – 1 for monopod, 3 for tripod, 4 for quadrupod.

• spacing_m (float, optional) – Center-to-leg distance in meters. Required for n_buckets > 1.

• soil_type (as before.)

• G (as before.)

• nu (as before.)

Returns:

• numpy.ndarray – 6x6 diagonal stiffness matrix.

• Reference

• ———

• DNVGL-ST-0126 (2021), Annex F “Suction bucket foundations”

• Gazetas, G. (1991). “Formulas and charts for impedances of – surface and embedded foundations”. J. Geotech. Eng. 117(9), 1363-1381.

ISO 19901-4 Annex E — Stiffness coefficients for shallow and pile foundations of offshore structures.

Reference:

ISO 19901-4:2016 “Petroleum and natural gas industries — Specific requirements for offshore structures — Part 4: Geotechnical and foundation design considerations”. International Organization for Standardization.

ISO 19901-4 is the international consensus standard for offshore geotechnical foundation design and is the reference adopted by DNVGL-ST-0126 §5.5 for offshore wind turbines. The expressions in this module are reproduced from Annex E of ISO 19901-4.

The expressions are similar to (and in many cases identical to) the Gazetas (1991) closed-form impedance formulas, with corrections for embedment, layered soil profiles, and dynamic effects. For Op^3 Mode B (linear stiffness), only the static expressions are used.

op3.standards.iso_19901_4.iso_pile_stiffness(diameter_m: float, embedment_m: float, soil_type: str = 'dense_sand', G: float | None = None, nu: float | None = None, pile_type: str = 'flexible') → ndarray

6x6 stiffness for a single pile per ISO 19901-4 Annex E.

Distinguishes ‘rigid’ (slender, kinematic rotation) and ‘flexible’ (long, elastic bending) pile responses.

Reference

ISO 19901-4:2016 Annex E “Pile foundations”.

op3.standards.iso_19901_4.iso_shallow_foundation_stiffness(radius_m: float, embedment_m: float = 0.0, shape: str = 'circular', soil_type: str = 'dense_sand', G: float | None = None, nu: float | None = None) → ndarray

6x6 stiffness for a shallow circular or square foundation per ISO 19901-4 Annex E Table E.1.

Parameters:

• radius_m (float) – For ‘circular’: the radius. For ‘square’: half-side length.

• embedment_m (float, default 0.0) – Embedment depth (0 = surface foundation).

• shape (str, default 'circular') – ‘circular’ or ‘square’. Square uses an equivalent radius.

• soil_type (see dnv_st_0126.)

• G (see dnv_st_0126.)

• nu (see dnv_st_0126.)

• Reference

• ---------

• 19901-4 (ISO) – foundations”.

API RP 2GEO — Geotechnical and Foundation Design Considerations.

Reference:

API RP 2GEO (2011, addendum 2014). “Geotechnical and Foundation Design Considerations”. American Petroleum Institute, 1st Edition.

API RP 2GEO is the American consensus standard for offshore geotechnical foundation design. Section 8 provides closed-form impedance functions for shallow and pile foundations that are essentially the Gazetas (1991) formulas with API-specific shape correction factors.

This module also provides the full Gazetas 6x6 coupled stiffness matrix including translational-rotational coupling terms, which is the most accurate analytical 6x6 stiffness available for shallow foundations and is the recommended Mode B input when no site-specific PISA or FE analysis is available.

op3.standards.api_rp_2geo.api_pile_stiffness(diameter_m: float, embedment_m: float, soil_type: str = 'dense_sand', G: float | None = None, nu: float | None = None) → ndarray

6x6 pile stiffness per API RP 2GEO Section 8.3.

Reference

API RP 2GEO (2011) Section 8.3 “Linear soil stiffness for foundation

impedance analysis”.

op3.standards.api_rp_2geo.gazetas_full_6x6(radius_m: float, embedment_m: float = 0.0, soil_type: str = 'dense_sand', G: float | None = None, nu: float | None = None) → ndarray

Full 6x6 coupled stiffness matrix per Gazetas (1991).

Includes translational-rotational coupling terms K_xz, K_yz that are zero in the diagonal models. The coupling terms become significant for embedded foundations and are the right representation when used as a SubDyn 6x6 in OpenFAST.

Reference

Gazetas, G. (1991). “Formulas and charts for impedances of surface

and embedded foundations”. Journal of Geotechnical Engineering 117(9), 1363-1381.

Carbon Trust Offshore Wind Accelerator (OWA) suction bucket foundation guidance.

Reference:

Carbon Trust OWA (2019). “Bearing Capacity Report — Suction Bucket Jackets”. The Carbon Trust, OWA Joint Industry Project.

Houlsby, G. T., & Byrne, B. W. (2005). “Design procedures for installation of suction caissons in clay and other materials”. Proc. ICE — Geotechnical Engineering, 158(2), 75-82.

Houlsby, G. T., & Byrne, B. W. (2005). “Design procedures for installation of suction caissons in sand”. Proc. ICE — Geotechnical Engineering, 158(3), 135-144.

The OWA Bearing Capacity Report is the most authoritative open guidance on suction bucket foundation design and is the standard that the offshore wind industry adopted for the early commercial suction bucket projects (Borkum Riffgrund, Aberdeen Bay, Aberdeen Offshore Wind Farm). Op^3 Mode B can use these expressions for preliminary suction bucket stiffness when no site-specific PISA or 3D FE analysis is available.

op3.standards.owa_bearing.houlsby_byrne_caisson_stiffness(diameter_m: float, skirt_length_m: float, soil_type: str = 'soft_clay', G: float | None = None, nu: float | None = None) → ndarray

Direct call to the single-caisson Houlsby & Byrne (2005) formulas.

Convenience alias for owa_suction_bucket_stiffness(n_buckets=1, …).

op3.standards.owa_bearing.owa_suction_bucket_stiffness(diameter_m: float, skirt_length_m: float, n_buckets: int = 1, spacing_m: float | None = None, soil_type: str = 'soft_clay', G: float | None = None, nu: float | None = None) → ndarray

6x6 stiffness for a suction bucket per OWA Bearing Capacity Report.

The OWA expressions are based on the Houlsby & Byrne (2005) caisson impedance functions extended to multi-bucket configurations (tripod, quadrupod) by treating each bucket as an independent foundation and combining via parallel axis theorem.

Parameters:

• diameter_m (float) – Bucket outer diameter.

• skirt_length_m (float) – Skirt embedment depth.

• n_buckets (int, default 1) – Number of buckets (1, 3, or 4).

• spacing_m (float, optional) – Center-to-leg distance for multi-bucket configurations.

• soil_type (as in dnv_st_0126.)

• G (as in dnv_st_0126.)

• nu (as in dnv_st_0126.)

• Reference

• ---------

• (2019) (Carbon Trust OWA) – for design”

• stiffness (Section 4.3 "Equivalent linear) – for design”

• Houlsby

• T. (G.)

• Byrne (&)

• (2005) (B. W.)

• 5 (Section)

PISA framework (Burd 2020 / Byrne 2020)

PISA monopile soil reaction framework (Phase 3 / Task 3.1).

Implements the PISA (Pile Soil Analysis) one-dimensional macro-element formulation for laterally loaded monopiles in sand and clay. PISA replaces the legacy API RP 2GEO p-y curves with FOUR distributed reaction components per pile depth segment:

1. p(z, v) distributed lateral load [N/m] f(local lateral disp v)

2. m(z, psi) distributed moment [N] f(local rotation psi)

3. H_b(v_b) pile-base shear [N] f(base lateral disp)

4. M_b(psi_b) pile-base moment [Nm] f(base rotation)

The non-dimensional shape functions are 4-parameter conic curves:

y / y_u = ((x_u - x) / (x_u - x_y)) - sqrt(((x_u - x) / (x_u - x_y))^2 - 4 n x / x_u * ((x_u - x_y) / x_u))

with parameters (k, n, x_u, y_u) calibrated for each component and each soil type. This module exposes:

pisa_lateral_pl(z, v, D, L, soil) -> p [N/m] pisa_moment_pl (z, psi, D, L, soil) -> m [N] pisa_base_shear (v_b, D, soil) -> H_b [N] pisa_base_moment(psi_b, D, soil) -> M_b [Nm] pisa_pile_stiffness_6x6(D, L, soil_profile) -> 6x6 K via small-strain

References

Burd, H. J. et al. (2020). “PISA design model for monopiles for

offshore wind turbines: application to a stiff glacial clay till”. Geotechnique 70(11), 1030-1047.

Byrne, B. W. et al. (2020). “PISA design model for monopiles for

offshore wind turbines: application to a marine sand”. Geotechnique 70(11), 1048-1066.

McAdam, R. A. et al. (2020). “Monotonic laterally loaded pile testing

in a dense marine sand at Dunkirk”. Geotechnique 70(11), 986-998.

Calibration regime

The published PISA coefficients are calibrated for:

• L/D = 2..10 (rigid to semi-rigid monopiles)

• D = 2..10 m

The parameters in this module are reproduced from Burd 2020 Table 6 (clay) and Byrne 2020 Table 7 (sand). Other soils require independent calibration.

class op3.standards.pisa.SoilState(depth_m: float, G_Pa: float, su_or_phi: float, soil_type: Literal['sand', 'clay'])

Bases: object

Local soil properties at a given depth.

G_Pa: float

depth_m: float

soil_type: Literal['sand', 'clay']

su_or_phi: float

op3.standards.pisa.conic(x: float, k: float, n: float, x_u: float, y_u: float) → float

PISA 4-parameter conic function.

y / y_u = c1 - sqrt(c1^2 - 4 n (x / x_u) c2)

where:

• c1 = 1 + n (1 - x / x_u)

• c2 = 1 - n

For n -> 0 the curve is bilinear (k for x < x_y, then plateau at y_u). For n -> 1 it is asymptotically elastic-perfectly plastic.

op3.standards.pisa.conic_initial_slope(k: float, y_u: float, x_u: float) → float

Initial slope at x = 0 of the normalised conic = k * y_u / x_u.

op3.standards.pisa.effective_head_stiffness(K: ndarray, h_load_m: float) → float

Secant stiffness H / v_G for a horizontal load H applied at height h_load_m above the ground-line reference of the 6x6 pile-head stiffness matrix.

For a rigid-pile 2x2 {translation, rotation} block

[ Kxx Kxrx ] [v ] [ H ] [ Kxrx Krxrx] [psi] = [ H h ]

the ground-level translation is

v = H (Krxrx - h Kxrx) / det

where det = Kxx Krxrx - Kxrx^2. Therefore

k_eff = H / v = det / (Krxrx - h Kxrx).

This is the apples-to-apples comparator for the McAdam 2020 and Byrne 2020 field-test k_Hinit values, which are defined as the secant slope of H vs v_G at small displacement.

op3.standards.pisa.pisa_base_moment(psi_b: float, D: float, soil: SoilState, L: float = 0.0) → float

Pile-base moment M_b [Nm] as a function of base rotation.

op3.standards.pisa.pisa_base_shear(v_b: float, D: float, soil: SoilState, L: float = 0.0) → float

Pile-base shear H_b [N] as a function of base lateral displacement.

op3.standards.pisa.pisa_coeffs(component: str, soil_type: str, z_over_D: float = 0.0, L_over_D: float = 0.0) → dict

Evaluate the depth-function-adjusted conic parameters for a specific soil reaction component at a given depth.

For distributed components lateral_p and moment_m the variable is z/D. For base components base_shear and base_moment the variable is L/D (the full pile slenderness) and the z/D argument is ignored.

op3.standards.pisa.pisa_lateral_pl(z: float, v: float, D: float, L: float, soil: SoilState) → float

Distributed lateral reaction p [N/m] at depth z, displacement v.

op3.standards.pisa.pisa_moment_pl(z: float, psi: float, D: float, L: float, soil: SoilState) → float

Distributed moment reaction m [N] at depth z, rotation psi.

op3.standards.pisa.pisa_pile_stiffness_6x6(diameter_m: float, embed_length_m: float, soil_profile: list[SoilState], n_segments: int = 50) → ndarray

Initial small-strain 6x6 K matrix at the pile head (mudline) by integrating the PISA distributed reactions along the embedded length and adding the base contributions.

The cross-coupling K_xrx (lateral-rocking) is computed via the second moment of the distributed lateral stiffness:

K_xx = int_0^L k_p(z) dz + K_b_shear K_xrx = int_0^L k_p(z) z dz (-) K_rxrx = int_0^L k_p(z) z^2 dz + int_0^L k_m(z) dz + K_b_moment

Vertical and torsional terms use elementary skin-friction estimates since PISA does not address them; this is consistent with how monopile designers extend PISA in DNV-ST-0126 commentary.

Cyclic degradation (Hardin-Drnevich / Vucetic-Dobry)

Cyclic soil stiffness degradation (Phase 3 / Task 3.2).

Implements the Hardin-Drnevich (1972) modified-hyperbolic backbone curve and the Vucetic-Dobry (1991) plasticity-index-dependent reference shear strain, plus convenience wrappers that take a layered soil profile and an applied (or estimated) cyclic shear strain and return a degraded profile suitable for re-running the PISA Mode B 6x6 stiffness assembly under design-storm conditions.

Backbone curve

The modified hyperbolic Hardin-Drnevich:

G / G_max = 1 / (1 + (gamma / gamma_ref) ** a)

with curvature exponent a (default 1.0 for the original 1972 form). For a = 0.92 the curve matches the Darendeli (2001) formulation widely used in seismic site response analysis.

Reference strain

Vucetic & Dobry (1991) Figure 5: gamma_ref increases with plasticity index (PI). For sands (PI = 0): gamma_ref ~ 1e-4. For high-plasticity clays (PI = 200): gamma_ref ~ 5e-3. Linear-in-log fit between.

References

Hardin, B. O., & Drnevich, V. P. (1972). “Shear modulus and damping

in soils: design equations and curves”. J. Soil Mech. Found. Div. ASCE, 98(7), 667-692.

Vucetic, M., & Dobry, R. (1991). “Effect of soil plasticity on cyclic

response”. J. Geotech. Eng., 117(1), 89-107.

Darendeli, M. B. (2001). “Development of a new family of normalized

modulus reduction and material damping curves”. PhD dissertation, Univ. of Texas at Austin.

op3.standards.cyclic_degradation.cyclic_stiffness_6x6(*, diameter_m: float, embed_length_m: float, soil_profile: list[SoilState], cyclic_strain: float, PI_percent: float | None = None, a: float = 1.0, n_segments: int = 50) → ndarray

Convenience wrapper: degrade the profile via Hardin-Drnevich, then call pisa_pile_stiffness_6x6 on the degraded profile.

op3.standards.cyclic_degradation.damping_ratio(gamma: float, gamma_ref: float, D_min: float = 0.01, D_max: float = 0.2) → float

Hyperbolic damping ratio companion to the Hardin-Drnevich modulus reduction. Linearly interpolates between D_min (small strain) and D_max (large strain) using the same gamma_ref.

op3.standards.cyclic_degradation.degrade_profile(soil_profile: list[SoilState], cyclic_strain: float, PI_percent: float | None = None, a: float = 1.0) → list[SoilState]

Apply Hardin-Drnevich modulus reduction layer-by-layer to a soil profile and return a NEW profile with degraded G values.

The original profile is not mutated. The new profile is suitable for re-running pisa_pile_stiffness_6x6 to obtain a cyclic / strain-softened 6x6 K matrix.

Parameters:

• soil_profile – List of SoilState (same convention used by op3.standards.pisa).

• cyclic_strain – The shear strain at which the modulus is to be evaluated. Engineering strain, dimensionless.

• PI_percent – Plasticity index for clay layers. Ignored for sand layers.

• a – Hardin-Drnevich curvature exponent (default 1.0).

Returns:

Degraded profile with G_Pa replaced by G_max * (G/G_max).

Return type:

list[SoilState]

op3.standards.cyclic_degradation.gamma_ref_for(soil: SoilState, PI_percent: float | None = None) → float

Pick a reference strain for a given SoilState. For sand (the Op^3 SoilState convention is su_or_phi = friction angle), use the sand baseline 1e-4. For clay use the Vucetic-Dobry curve with the user-supplied PI; if PI is not provided, default to PI = 30%.

op3.standards.cyclic_degradation.hardin_drnevich(gamma: float, gamma_ref: float, a: float = 1.0) → float

Modulus reduction G / G_max for shear strain gamma.

Parameters:

• gamma – Cyclic shear strain (engineering strain, dimensionless).

• gamma_ref – Reference shear strain at which G / G_max = 0.5.

• a – Curvature exponent. Default 1.0 (Hardin-Drnevich 1972). Use 0.92 for the Darendeli 2001 form.

Returns:

G / G_max in [0, 1].

Return type:

float

op3.standards.cyclic_degradation.hardin_drnevich_array(gamma: ndarray, gamma_ref: float, a: float = 1.0) → ndarray

op3.standards.cyclic_degradation.vucetic_dobry_gamma_ref(PI_percent: float) → float

Reference shear strain as a function of plasticity index.

Linear-in-log interpolation through the Vucetic-Dobry (1991) Figure 5 curve. PI_percent = 0 recovers the sand value (gamma_ref ~ 1e-4) and PI_percent = 200 saturates at ~5e-3.

HSsmall constitutive wrapper

HSsmall constitutive wrapper (Phase 3 / Task 3.3).

Bridges the OptumGX Hardening-Soil-with-small-strain-stiffness output to the Op^3 PISA / cyclic-degradation pipeline.

The Hardening Soil model with small-strain stiffness (Benz 2007; Schanz et al. 1999) parameterises the soil with:

E_50_ref reference secant stiffness at 50% strength E_oed_ref oedometric reference stiffness E_ur_ref unload-reload reference stiffness G_0_ref small-strain shear modulus at reference pressure gamma_07 shear strain at G/G_0 = 0.7 m power-law stress exponent c, phi, psi strength parameters p_ref reference pressure (usually 100 kPa)

For Op^3 the relevant outputs are:

G_0(z) = G_0_ref * ((c cot(phi) + sigma_3’(z)) / (c cot(phi) + p_ref)) ** m

which gives a depth-dependent small-strain shear modulus that feeds the PISA framework directly. For sand (c = 0):

G_0(z) = G_0_ref * (sigma_3’(z) / p_ref) ** m

This module exposes:

• HSsmallParams : dataclass for HSsmall material parameters

• hssmall_G_at_depth() : evaluate G_0(z) for one HSsmall material

• load_hssmall_profile() : read an OptumGX-exported CSV -> list[SoilState]

• hssmall_to_pisa() : convert an HSsmall material list to a PISA-ready soil profile

Reference

Benz, T. (2007). “Small-strain stiffness of soils and its numerical

consequences”. PhD dissertation, Univ. of Stuttgart.

Schanz, T., Vermeer, P. A., & Bonnier, P. G. (1999). “The hardening

soil model: formulation and verification”. Beyond 2000 in Computational Geotechnics, 281-296.

class op3.standards.hssmall.HSsmallParams(layer_name: str, z_top_m: float, z_bot_m: float, soil_type: str, G0_ref_Pa: float, p_ref_Pa: float = 100000.0, m_exp: float = 0.5, phi_deg: float = 0.0, c_Pa: float = 0.0, su_Pa: float = 0.0, gamma_07: float = 0.0001, PI_percent: float | None = None, notes: str = '')

Bases: object

One HSsmall material layer.

G0_ref_Pa: float

PI_percent: float | None = None

c_Pa: float = 0.0

gamma_07: float = 0.0001

layer_name: str

m_exp: float = 0.5

notes: str = ''

p_ref_Pa: float = 100000.0

phi_deg: float = 0.0

soil_type: str

su_Pa: float = 0.0

z_bot_m: float

z_top_m: float

op3.standards.hssmall.hssmall_G_at_depth(params: HSsmallParams, z_m: float, gamma_eff_kN_per_m3: float = 10.0, K0: float = 0.5) → float

Stress-dependent small-strain shear modulus G_0(z) per the HSsmall power law.

For sand (c = 0): G_0(z) = G_0_ref * (sigma_3’/p_ref) ** m For clay (phi = 0, c > 0): use cohesion-shifted form.

op3.standards.hssmall.hssmall_to_pisa(layers: list[HSsmallParams], n_points_per_layer: int = 3) → list[SoilState]

Sample each HSsmall layer at n_points_per_layer depths and return a flat list of SoilState records that PISA can integrate.

op3.standards.hssmall.load_hssmall_profile(csv_path: str | Path) → list[HSsmallParams]

Read an OptumGX HSsmall layer export.

Expected columns (case-sensitive aliases supported):

layer_name, z_top_m, z_bot_m, soil_type, G0_ref_Pa, [p_ref_Pa], [m_exp], [phi_deg], [c_Pa], [su_Pa], [gamma_07], [PI_percent]

Anything in [brackets] is optional and falls back to the HSsmallParams dataclass default.

OpenSees foundation builder

Op^3 OpenSeesPy model builder.

Implements the four foundation modes as OpenSeesPy element patterns that all share a common tower-base interface. The public entry points are called by the composer and cross-compare utilities; a user should not need to call anything in this module directly.

Tag allocation convention (avoids collisions across rebuilds):

1..999 reserved for tower element nodes 1000..1099 tower base node (1000) + tower element tags (1001..1099) 1100..1199 foundation interface nodes (6x6 Mode, BNWF anchor) 1200..1999 BNWF spring soil nodes (Modes C, D) 2000..2999 BNWF zero-length elements 3000..3099 rotor-nacelle-assembly mass node 9000..9999 temporary / analysis-only tags

op3.opensees_foundations.builder.attach_foundation(foundation: Foundation, base_node: int) → dict

Attach a foundation to the OpenSees domain at the given base node.

op3.opensees_foundations.builder.build_opensees_model(tower_model: TowerModel) → None

Build the OpenSeesPy domain for a TowerModel.

This function is called by TowerModel.build(). It wipes the current OpenSees domain, instantiates the tower stick model, attaches the foundation according to the chosen mode, and places the RNA mass at the hub node.

op3.opensees_foundations.builder.run_cyclic_analysis(tower_model: TowerModel, n_cycles: int = 10, amplitude_m: float = 0.5, n_steps_per_half: int = 25, load_node: int | None = None) → dict

Run cyclic lateral push-pull at the hub and track degradation.

Applies symmetric cyclic displacement at the hub node (lateral DOF 1) and records permanent (residual) rotation at the tower base and the secant stiffness after each full cycle.

Parameters:

• tower_model (TowerModel) – A built TowerModel.

• n_cycles (int) – Number of full push-pull cycles.

• amplitude_m (float) – Peak lateral displacement at the hub per half-cycle [m].

• n_steps_per_half (int) – Analysis steps per half-cycle (push or pull).

• load_node (int or None) – Hub node tag. Auto-detected if None.

Returns:

• cycle_number : list[int]

• permanent_rotation_deg : list[float] – residual base rotation after returning to zero hub displacement.

• stiffness_kN_m_per_deg : list[float] – secant stiffness per cycle.

• peak_force_kN : list[float] – peak lateral reaction per cycle.

Return type:

dict

op3.opensees_foundations.builder.run_eigen_analysis(tower_model: TowerModel, n_modes: int = 6) → np.ndarray

Run eigenvalue analysis and return natural frequencies in Hz.

Tries the default -frequency (Arpack) solver first for speed. Arpack state does not always survive a second eigen call in the same Python process; when it returns NaN or zero frequencies we retry with -fullGenLapack which is deterministic but slower.

op3.opensees_foundations.builder.run_pushover_analysis(tower_model: TowerModel, target_disp_m: float = 1.0, n_steps: int = 50, load_node: int | None = None) → dict

Run a static lateral pushover at the hub node and return the force-displacement curve.

Returns a dict with keys ‘displacement_m’ and ‘reaction_kN’.

op3.opensees_foundations.builder.run_pushover_moment_rotation(tower_model: TowerModel, target_rotation_rad: float = 0.05, n_steps: int = 100) → dict

Run a moment-controlled pushover at the tower base.

Applies an incrementally increasing rotation at the tower base node (1000) and records the resisting moment. This is the standard test for foundation moment capacity and stiffness degradation under monotonic loading.

Parameters:

• tower_model (TowerModel) – A built TowerModel (tower_model.build() must have been called).

• target_rotation_rad (float) – Maximum imposed rotation at the base [rad]. Default 0.05 (~2.9 deg).

• n_steps (int) – Number of load increments.

Returns:

rotation_deg : list[float] — base rotation in degrees. moment_MNm : list[float] — resisting moment in MN*m. base_node : int — node tag where rotation was imposed.

Return type:

dict

op3.opensees_foundations.builder.run_static_condensation(tower_model: TowerModel) → np.ndarray

Extract the 6x6 stiffness matrix at the tower base / foundation interface via analytic Winkler integration.

The earlier version of this function used a finite-difference probe with SP-imposed unit displacements at the base node; that approach failed on Mode C (distributed BNWF) models because the springs absorb the imposed displacement and the nodal reactions read back as ~zero. The analytic Winkler integral is the correct closed-form condensation for any elastic Winkler foundation and matches the same expression used internally by pisa_pile_stiffness_6x6.

For Mode A (fixed) a diagonal dummy matrix with 1e20 is returned (effectively rigid). For Mode B the stored stiffness_matrix is returned directly. For Modes C / D the spring table is integrated via the Winkler formula:

K_xx = sum k(z) dz K_xrx = sum k(z) z dz K_rxrx = sum k(z) z^2 dz

op3.opensees_foundations.builder.run_transient_analysis(tower_model: TowerModel, duration_s: float = 10.0, dt_s: float = 0.01, damping_ratio: float = 0.01) → dict

Run a free vibration transient with an initial hub-node displacement.

Returns a dict with ‘time_s’ and ‘hub_disp_m’ arrays.

Tower template loader: parse OpenFAST ElastoDyn tower input files and return distributed tower properties suitable for assembling an OpenSeesPy stick model with byte-identical NREL section properties.

This is Phase 1 / Task 1.1 of the Op^3 Track C industry-grade plan: replace the hand-tuned TOWER_TEMPLATES dict in op3/opensees_foundations/builder.py with values read directly from the canonical NREL *_ElastoDyn*Tower.dat files. The expected outcome is that examples 1, 2, 3, 7, 8 frequency errors fall from 10-30% to <2% on the fixed-base configuration.

Usage

>>> from op3.opensees_foundations.tower_loader import load_elastodyn_tower
>>> tpl = load_elastodyn_tower(
...     ed_main="nrel_reference/openfast_rtest/5MW_OC3Mnpl_DLL_WTurb_WavesIrr/"
...             "NRELOffshrBsline5MW_OC3Monopile_ElastoDyn.dat",
... )
>>> tpl.tower_height_m, len(tpl.ht_fract)
(87.6, 11)

Schema

Returned TowerTemplate exposes:

• tower_height_m : TowerHt from ElastoDyn main file

• tower_base_z_m : TowerBsHt (offset from MSL or ground)

• ht_fract : array of fractional heights along tower (0..1)

• mass_density_kg_m : array of TMassDen at each station

• ei_fa_Nm2 : array of TwFAStif (fore-aft EI)

• ei_ss_Nm2 : array of TwSSStif (side-side EI)

• damping_ratio_fa : (TwrFADmp1, TwrFADmp2) percent of critical

• damping_ratio_ss : (TwrSSDmp1, TwrSSDmp2)

• adj_tower_mass : AdjTwMa scalar

• adj_fa_stiff : AdjFASt scalar

• adj_ss_stiff : AdjSSSt scalar

• source_files : provenance for V&V records

class op3.opensees_foundations.tower_loader.RNAProperties(hub_mass_kg: float, nac_mass_kg: float, blade_mass_kg: float, n_blades: int, hub_iner_kgm2: float, nac_yiner_kgm2: float, nac_cm_xn_m: float, nac_cm_zn_m: float, twr2shft_m: float, tip_rad_m: float, hub_rad_m: float, source_files: list[str] = <factory>)

Bases: object

Mass + inertia of the rotor-nacelle assembly parsed from ElastoDyn.

blade_mass_kg: float

hub_iner_kgm2: float

hub_mass_kg: float

hub_rad_m: float

n_blades: int

nac_cm_xn_m: float

nac_cm_zn_m: float

nac_mass_kg: float

nac_yiner_kgm2: float

source_files: list[str]

tip_rad_m: float

property total_rna_mass_kg: float

twr2shft_m: float

class op3.opensees_foundations.tower_loader.TowerTemplate(tower_height_m: float, tower_base_z_m: float, ht_fract: ndarray, mass_density_kg_m: ndarray, ei_fa_Nm2: ndarray, ei_ss_Nm2: ndarray, damping_ratio_fa: tuple[float, float]=(1.0, 1.0), damping_ratio_ss: tuple[float, float]=(1.0, 1.0), adj_tower_mass: float = 1.0, adj_fa_stiff: float = 1.0, adj_ss_stiff: float = 1.0, source_files: list[str] = <factory>)

Bases: object

Distributed tower properties parsed from an ElastoDyn deck.

adj_fa_stiff: float = 1.0

adj_ss_stiff: float = 1.0

adj_tower_mass: float = 1.0

damping_ratio_fa: tuple[float, float] = (1.0, 1.0)

damping_ratio_ss: tuple[float, float] = (1.0, 1.0)

ei_fa_Nm2: ndarray

ei_ss_Nm2: ndarray

ht_fract: ndarray

mass_density_kg_m: ndarray

property n_stations: int

section_at(ht_fract: float) → dict

Linearly interpolated section properties at a given height fraction.

source_files: list[str]

station_elevations_m() → ndarray

Absolute z-coordinates of stations (base + ht_fract * height).

tower_base_z_m: float

tower_height_m: float

op3.opensees_foundations.tower_loader.discretise(tpl: TowerTemplate, n_segments: int = 20) → list[dict]

Subdivide the tower into n_segments equal-length elements and return per-element section properties (interpolated at the element midpoint). The resulting list is what opensees_foundations.builder._build_tower_stick consumes.

Each entry has keys: z_bot, z_top, length, mass_kg_m, EI_fa_Nm2, EI_ss_Nm2, EI_avg_Nm2.

op3.opensees_foundations.tower_loader.evaluate_mode_shape(coeffs: ndarray, eta: ndarray) → ndarray

Sample psi(eta) given the [c2..c6] polynomial coefficients.

op3.opensees_foundations.tower_loader.load_elastodyn_rna(ed_main: str | Path) → RNAProperties

op3.opensees_foundations.tower_loader.load_elastodyn_tower(ed_main: str | Path, ed_tower: str | Path | None = None) → TowerTemplate

Parse an OpenFAST ElastoDyn deck and return a TowerTemplate.

Parameters:

• ed_main – Path to the ElastoDyn primary input file (the one referenced by EDFile in the top-level .fst). Provides TowerHt and TowerBsHt plus the path to the tower file via TwrFile.

• ed_tower – Optional explicit override for the tower properties file. If omitted, the loader resolves it from TwrFile in ed_main.

Notes

Some NREL decks (the OC3 monopile in particular) embed the distributed properties directly in the main ElastoDyn file rather than in a separate *_Tower.dat. We handle both layouts.

op3.opensees_foundations.tower_loader.published_mode_shape(twr_text: str, mode: str = 'TwFAM1Sh') → ndarray | None

Return the 5 polynomial coefficients [c2, c3, c4, c5, c6] for one of the NREL ElastoDyn tower mode shapes (TwFAM1Sh, TwFAM2Sh, TwSSM1Sh, TwSSM2Sh).

Mode shape is psi(eta) = sum_{i=2..6} c_i * eta^i where eta in [0,1].

OpenFAST coupling

Op^3 -> OpenFAST SoilDyn exporter (Phase 4 / SoilDyn bridge).

OpenFAST v5 introduced the SoilDyn module with three calculation options:

CalcOption = 1 Stiffness / Damping matrices (6x6 K + 6x6 D) CalcOption = 2 P-Y curves (currently unavailable) CalcOption = 3 Coupled REDWIN DLL (binary plug-in)

Op^3 already produces a calibrated 6x6 head stiffness via the PISA, DNV, ISO, API, OWA, and Mode-D pipelines. This module writes that matrix into the canonical SoilDyn input file format so that any Op^3 foundation can be plugged DIRECTLY into a SoilDyn-enabled OpenFAST deck without manual conversion.

Reference

Bergua, Robertson, Jonkman, Platt (2021). “Specification Document for

OC6 Phase II: Verification of an Advanced Soil-Structure Interaction Model for Offshore Wind Turbines”. NREL/TP-5000-79989. https://doi.org/10.2172/1811648

op3.openfast_coupling.soildyn_export.redwin_backbones_from_spring_table(spring_table, *, diameter_m: float, skirt_length_m: float | None = None, n_points: int = 25, max_disp_m: float = 0.2, max_rot_rad: float = 0.015) → dict

Build REDWIN Model-2 H-u and M-theta backbones from an Op^3 spring table.

This is a convenience adapter that approximates the blueprint’s Q1(b) macro-element calibration without running an OptumG2 VHM Limit-Analysis sweep. It integrates the per-depth lateral springs (PySimple1 conic backbone y = p*u / (1 + |u/y_r|)) analytically over the skirt length to produce head H-u and M-theta curves.

The resulting backbones are conservative (rigid-pile assumption) and should be replaced with OptumG2-calibrated curves once available. Until then the approximation is good to ~10-20% on ultimate capacity and ~30% on initial stiffness, which is within the REDWIN calibration envelope reported by Skau et al. 2018.

Returns:

K_seabed : 6x6 stiffness from Winkler integration. pushover_H_u : (n_points, 2) lateral backbone. pushover_M_theta: (n_points, 2) rocking backbone.

Return type:

dict

op3.openfast_coupling.soildyn_export.write_soildyn_from_pisa(out_path: str | Path, *, diameter_m: float, embed_length_m: float, soil_profile: list, location_xyz: tuple[float, float, float] = (0.0, 0.0, 0.0), n_segments: int = 50) → Path

Convenience: build a PISA K matrix and write it as a SoilDyn .dat in one call. Equivalent to calling pisa_pile_stiffness_6x6() then write_soildyn_input().

op3.openfast_coupling.soildyn_export.write_soildyn_input(out_path: str | Path, K: ndarray, *, location_xyz: tuple[float, float, float] = (0.0, 0.0, 0.0), damping: ndarray | None = None, provenance: str = 'Op^3 PISA') → Path

Write a SoilDyn input file using the Op^3 6x6 stiffness matrix.

Parameters:

• out_path – Destination file path.

• K – 6x6 stiffness matrix in SI units (N/m for translational, Nm/rad for rotational, off-diagonal in mixed units per OpenFAST convention).

• location_xyz – Coupling point in OpenFAST global coordinates (typically the tower base or the SubDyn interface node).

• damping – Optional 6x6 damping matrix. If None, a zero matrix is written.

• provenance – String describing where K came from (e.g. “Op^3 PISA Burd 2020” or “Op^3 Mode D alpha=2.0”).

op3.openfast_coupling.soildyn_export.write_soildyn_multipoint(out_path: str | Path, points: list[dict], *, provenance: str = 'Op^3 multi-point') → Path

Write a multi-point SoilDyn input file using CalcOption=3 layout (one independent 6x6 K per coupling point). Each points entry must be a dict with keys:

• location: (x, y, z) tuple

• K : 6x6 ndarray

• damping : optional 6x6 ndarray (defaults to zero)

• label : optional string label

This is the format used by the OC6 Phase II REDWIN DLL test case and is the natural target for the Op^3 Mode D dissipation-weighted custom DLL: each tripod leg gets its own K matrix derived from a different soil profile (e.g. one leg in dense sand, two legs in clay) and the dissipation weighting can be applied per-point.

Note: stock OpenFAST v5.0.0 SoilDyn requires the corresponding SubDyn mesh nodes to exist within 0.1 m of each location.

op3.openfast_coupling.soildyn_export.write_soildyn_redwin(out_path: str | Path, *, points: list[dict], dll_name: str = 'REDWIN.dll', model: int = 2, props_dir: str | Path | None = None, ldisp_dir: str | Path | None = None, provenance: str = 'Op^3 REDWIN shim', acknowledge_dll_missing: bool = False) → dict

Write a SoilDyn CalcOption=3 input deck wired to a REDWIN DLL.

CRITICAL PLACEHOLDER NOTICE (2026-04-20): the Op³ stack does not currently ship a REDWIN-compatible DLL. Running OpenFAST against the deck produced by this function WILL FAIL at DLL-load time. This function exists to publish the file-interface contract (deck + PropsFile + LDispFile) so that an NGI-supplied or Python- shim DLL can be dropped in later. Call sites MUST pass acknowledge_dll_missing=True to confirm they understand this is a research artefact, not a runnable coupling.

Each points entry must be a dict with:

• location : (x, y, z) tuple in OpenFAST global coordinates.

• label : short string label (e.g. "bucket_1_scour_0m").

• K_seabed : 6x6 elastic stiffness at the seabed.

• pushover_H_u : Nx2 backbone of lateral H-u curve.

• pushover_M_theta: Nx2 backbone of rocking M-theta curve.

Returns:

out_path : SoilDyn deck path (str). props_paths : list of REDWIN PropsFile paths. ldisp_paths : list of LDispFile stub paths. runnable_by_openfast: always False until a DLL is wired in. dll_status : "missing" or "configured" once DLL found.

Return type:

dict

Fatigue

Fatigue damage-equivalent load (DEL) computation module for Op^3.

Implements rainflow counting and DEL per DNV-RP-C203 / Hayman (2012).

DEL formula:

DEL = (sum(n_i * S_i^m) / N_eq) ^ (1/m)

where:

n_i = cycle count for stress range S_i (half-cycles count as 0.5) S_i = stress range (peak-to-valley amplitude, NOT half-amplitude) m = inverse Woehler slope (3-5 for steel, 8-12 for composites) N_eq = equivalent number of cycles (typically 1 Hz * T_lifetime)

References

• Hayman, G. (2012). MLife Theory Manual. NREL/TP-XXXX.

• DNV-RP-C203 (2021). Fatigue design of offshore steel structures.

• ASTM E1049-85. Standard practices for cycle counting in fatigue analysis.

op3.fatigue.compute_del(signal: ndarray, m: float, n_eq: float | None = None, dt: float | None = None, f_eq: float = 1.0) → float

Compute the damage-equivalent load for a time series.

Parameters:

• signal (array_like) – Load time series (e.g. TwrBsMyt in kN-m). 1-D array.

• m (float) –

  Inverse S-N slope (Woehler exponent). Typical values:

  • 3 for welded steel (DNV-RP-C203 detail categories)

  • 4 for base-metal steel

  • 10 for glass-fibre composites

• n_eq (float, optional) – Equivalent cycle count. If None, computed as f_eq * T where T = len(signal) * dt.

• dt (float, optional) – Time step in seconds. Required if n_eq is None.

• f_eq (float) – Equivalent frequency in Hz (default 1.0).

Returns:

del_value – Damage-equivalent load in the same units as signal.

Return type:

float

Raises:

ValueError – If neither n_eq nor dt is provided.

op3.fatigue.compute_del_multi_slope(signal: ndarray, m_values: list[float], n_eq: float | None = None, dt: float | None = None, f_eq: float = 1.0) → dict[float, float]

Compute DEL for multiple Woehler exponents.

Useful for comparing steel (m=3-5) vs composite (m=10) fatigue.

Returns a dict mapping m -> DEL.

op3.fatigue.rainflow_count(signal: ndarray) → Tuple[ndarray, ndarray]

Rainflow cycle counting using the rainflow package.

Falls back to a built-in 4-point algorithm (ASTM E1049-85) if the package is unavailable.

Parameters:

signal (array_like) – 1-D load time series.

Returns:

• ranges (ndarray) – Stress ranges (peak-to-valley, always positive).

• counts (ndarray) – Cycle counts (0.5 for half-cycles, 1.0 for full cycles).

Visualization

Op3 Structural Visualization via opsvis.

Provides publication-quality figures of the OpenSeesPy model:

• Model geometry (nodes, elements, supports)

• Mode shapes from eigenvalue analysis

• Pushover deformed shape

• Section force diagrams (moment, shear, axial)

Usage:

from op3 import build_foundation, compose_tower_model
from op3.visualization import plot_all

fnd = build_foundation(mode="fixed")
model = compose_tower_model(rotor="nrel_5mw_baseline",
                            tower="nrel_5mw_tower",
                            foundation=fnd)
freqs = model.eigen(n_modes=3)
plot_all(model, freqs, output_dir="figures/")

All functions accept an optional output_dir to save PNG files. If not provided, figures are displayed interactively (if backend allows).

op3.visualization.plot_all(model=None, freqs: list[float] | None = None, pushover_result: dict | None = None, mr_result: dict | None = None, output_dir: str | Path | None = None, n_modes: int = 3) → dict[str, str | None]

Generate all available visualization figures.

Returns a dict of {name: filepath_or_None}.

op3.visualization.plot_deformed(scale_factor: float = 0.0, output_dir: str | Path | None = None, az_el: tuple[float, float] = (-60.0, 30.0), title: str = 'Deformed Shape') → str | None

Plot deformed shape after a static analysis (pushover, etc.).

scale_factor=0 means auto-scale.

op3.visualization.plot_mode_shapes(n_modes: int = 3, output_dir: str | Path | None = None, az_el: tuple[float, float] = (-60.0, 30.0), freqs: list[float] | None = None) → list[str | None]

Plot mode shapes from eigenvalue analysis.

Must be called after model.eigen(n_modes >= n_modes).

op3.visualization.plot_model(output_dir: str | Path | None = None, az_el: tuple[float, float] = (-60.0, 30.0), title: str = 'Op3 Structural Model') → str | None

Plot the OpenSeesPy model geometry (nodes, elements, supports).

Must be called after model.eigen() or model.build() so the OpenSees domain is populated.

op3.visualization.plot_moment_rotation(mr_result: dict, output_dir: str | Path | None = None, title: str = 'Moment-Rotation at Foundation Head', ref_My: float | None = None, ref_theta_y: float | None = None) → str | None

Plot M-theta from moment-rotation analysis.

op3.visualization.plot_pushover_curve(pushover_result: dict, output_dir: str | Path | None = None, title: str = 'Pushover Curve') → str | None

Plot force-displacement from pushover analysis.

op3.visualization.plot_section_forces(force_type: str = 'Mz', output_dir: str | Path | None = None, az_el: tuple[float, float] = (-60.0, 30.0), title: str | None = None) → str | None

Plot section force diagram after static analysis.

force_type: ‘N’ (axial), ‘Vy’/’Vz’ (shear), ‘My’/’Mz’ (moment), ‘T’ (torsion)

Op3 OptumGX Visualization via PyVista.

3D visualization of finite-element limit analysis results:

• Suction bucket mesh (skirt + lid + soil)

• Contact pressure contours on plate elements

• Bearing capacity factor Np(z) depth profile

• Plastic dissipation field (collapse mechanism)

• Velocity vector glyphs at collapse

• Spring profile visualization (k and p_ult vs depth)

Usage:

from op3.viz_optumgx import plot_bucket_pressure, plot_spring_profile
plot_spring_profile("data/fem_results/spring_profile_op3.csv",
                    output_dir="validation/figures/")

op3.viz_optumgx.create_bucket_mesh(D: float = 8.0, L: float = 9.3, n_theta: int = 24, n_z: int = 20) → PolyData

Create a 3D suction bucket mesh (lid + skirt wall).

Returns a PyVista PolyData with the bucket surface mesh.

op3.viz_optumgx.create_soil_block(D: float = 8.0, L: float = 9.3, domain_factor: float = 3.0) → PolyData

Create a soil block around the bucket for context.

op3.viz_optumgx.plot_bucket_pressure(plate_df=None, D: float = 8.0, L: float = 9.3, output_dir: str | Path | None = None, title: str = 'Contact Pressure at Collapse') → str | None

Plot contact pressure on bucket surface from OptumGX plate data.

plate_df: DataFrame with columns Xc, Yc, Zc, sig_net (or pressure). If None, generates a synthetic demonstration.

op3.viz_optumgx.plot_collapse_mechanism(D: float = 8.0, L: float = 9.3, output_dir: str | Path | None = None, title: str = 'Collapse Mechanism (Plastic Dissipation)') → str | None

Plot the plastic dissipation field around the bucket.

Uses a synthetic dissipation field for demonstration. Real data would come from OptumGX solid element extraction.

op3.viz_optumgx.plot_np_profile(pult_csv: str | Path, output_dir: str | Path | None = None, title: str = 'Bearing Capacity Factor Np(z) from OptumGX') → str | None

Plot the normalized lateral bearing capacity factor vs depth.

op3.viz_optumgx.plot_spring_profile(spring_csv: str | Path, output_dir: str | Path | None = None, title: str = 'OptumGX-Derived Spring Profile') → str | None

Plot k(z) and p_ult(z) spring profile from Op3 CSV.

Op3 OpenFAST Postprocessing via welib + pCrunch.

Advanced visualization and analysis of OpenFAST outputs:

• Time series with channel overlay

• Power spectral density (PSD)

• Power curve from DLC sweep

• Fatigue damage-equivalent loads (DEL)

• DLC batch statistics (mean, std, min, max, DEL per wind speed)

• Campbell diagram (from linearization files)

Usage:

from op3.viz_openfast import (
    plot_time_series, plot_psd, plot_power_curve,
    plot_del_bar, compute_dlc_statistics,
)
stats = compute_dlc_statistics("validation/dlc11_partial/")
plot_power_curve(stats, output_dir="validation/figures/")

op3.viz_openfast.compute_dlc_statistics(dlc_dir: str | Path, channels: list[str] = None) → dict

Compute statistics from a DLC sweep directory.

Returns dict with ‘wind_speeds’, ‘stats’ (DataFrame per channel).

op3.viz_openfast.plot_del_bar(outb_paths: list[str | Path], channels: list[str] = None, m_values: list[float] = None, output_dir: str | Path | None = None, title: str = 'Damage-Equivalent Loads') → str | None

Compute and plot DEL for multiple channels across runs.

op3.viz_openfast.plot_power_curve(stats: dict, output_dir: str | Path | None = None, title: str = 'Power Curve (DLC 1.1)') → str | None

Plot power curve from DLC statistics.

op3.viz_openfast.plot_psd(outb_path: str | Path, channel: str = None, output_dir: str | Path | None = None, title: str = None) → str | None

Plot PSD of a channel using Welch’s method.

op3.viz_openfast.plot_time_series(outb_path: str | Path, channels: list[str] = None, output_dir: str | Path | None = None, title: str = 'OpenFAST Time Series') → str | None

Plot selected channels from a single OpenFAST output file.

op3.viz_openfast.run_pcrunch_batch(outb_files: list[str | Path], wind_speeds: list[float], m_values: list[float] = None) → dict

Run pCrunch batch statistics on a set of OpenFAST outputs.

Returns dict with summary_stats and DELs DataFrames.

Op3 Tier 1 Visualization: Defense-Quality Figures.

Four high-impact figures that tell the thesis story:

1. VHM failure envelope with scour degradation

2. Cross-pipeline composite (OptumGX + OpenSeesPy + OpenFAST)

3. Scour progression sweep (frequency + mode shape evolution)

4. Mode C vs Mode D dissipation comparison

Usage:

python -m op3.viz_tier1 –output validation/figures/tier1/

op3.viz_tier1.fig_cross_pipeline(output_dir: Path) → str

3-panel figure showing the full Op3 pipeline.

1. OptumGX: spring profile with bucket sketch

2. OpenSeesPy: eigenvalue mode shape

3. OpenFAST: time series with frequency highlight

op3.viz_tier1.fig_mode_cd_comparison(output_dir: Path) → str

2-panel figure comparing Mode C (elastic) vs Mode D (dissipation).

Left: spring profile with color-coded w(z) weighting Right: frequency comparison at multiple alpha values

op3.viz_tier1.fig_scour_sweep(output_dir: Path) → str

4-panel figure showing mode shape evolution with scour.

Uses Op3 eigenvalue at 4 scour depths: S/D = 0, 0.1, 0.3, 0.5.

op3.viz_tier1.fig_vhm_envelope(output_dir: Path) → str

3-panel VHM failure envelope showing scour shrinkage.

Panel (a): V-H envelope at S/D = 0, 0.25, 0.5, 1.0 Panel (b): Capacity reduction vs S/D (V and H separately) Panel (c): Factor of Safety trajectory with scour

op3.viz_tier1.main()

Op3 Tier 2 Visualization: Journal-Paper Quality Figures.

5. Foundation sketch + depth profile overlay (geotechnical standard)

6. Rainflow cycle matrix heatmap

7. Campbell diagram (from .lin files)

8. M-theta backbone with published reference points

Usage:

python -m op3.viz_tier2 –output validation/figures/tier2/

op3.viz_tier2.fig_campbell_diagram(output_dir: Path) → str

Campbell diagram showing natural frequencies vs rotor speed.

Since Op3 linearization runs are not available, constructs a parametric Campbell from eigenvalue analyses at different RPM.

op3.viz_tier2.fig_foundation_profile(output_dir: Path) → str

Publication-quality foundation cross-section with depth profiles.

Left panel: bucket cross-section with soil layers and CPT-style G(z) Right panels: k(z) stiffness and p_ult(z) capacity as smooth curves

op3.viz_tier2.fig_moment_rotation(output_dir: Path) → str

Moment-rotation backbone curve with published reference points.

Overlays: DJ Kim 2014 centrifuge, Houlsby 2005 field trial, Barari 2021 FE prediction.

op3.viz_tier2.fig_rainflow_heatmap(output_dir: Path) → str

2D heatmap of rainflow cycles (range vs mean) from OpenFAST output.

Shows where fatigue damage concentrates in the load spectrum.

op3.viz_tier2.main()

Op3 Tier 3 Visualization: Interactive Web Dashboard Components.

9. Interactive 3D foundation model (Plotly)

10. Live sensor overlay with Bayesian prediction band

These generate standalone HTML files or Plotly figures that can be embedded in the op3_viz Dash app or Quarto reports.

Usage:

python -m op3.viz_tier3 –output validation/figures/tier3/

op3.viz_tier3.fig_interactive_3d(output_dir: Path) → str

Interactive 3D model with bucket mesh, springs, and tower.

Generates a standalone HTML file with rotation, zoom, and click-to-inspect functionality.

op3.viz_tier3.fig_sensor_overlay(output_dir: Path) → str

Frequency prediction band with field measurement overlay.

Shows Op3 prior/posterior distribution vs field-measured f1, demonstrating the digital twin concept.

op3.viz_tier3.main()

Uncertainty quantification (Phase 5)

Soil parameter UQ propagation (Phase 5 / Task 5.1).

Monte Carlo propagation of geotechnical input uncertainty through the PISA / cyclic / HSsmall pipeline. Samples are drawn from a layered soil prior and the resulting 6x6 head-stiffness distribution is summarised by mean, std, and 5/50/95 percentiles for each diagonal term.

The motivating use case: published soil investigation reports give G_max as a band (e.g. 60-120 MPa) rather than a point value. This module turns that band into a defensible distribution on the foundation stiffness that downstream Op^3 stages (eigenvalue, DLC simulation, fatigue) can consume.

Reference

Phoon, K. K., & Kulhawy, F. H. (1999). “Characterization of

geotechnical variability”. Canadian Geotechnical Journal, 36(4), 612-624. – baseline COVs for soil parameters.

class op3.uq.propagation.SoilPrior(depth_m: float, G_mean_Pa: float, G_cov: float = 0.3, soil_type: str = 'sand', su_or_phi_mean: float = 35.0, su_or_phi_cov: float = 0.1)

Bases: object

Probabilistic specification for one soil layer. The prior is parameterised by mean and coefficient of variation (COV = std / mean) for each fluctuating quantity. Lognormal sampling is used so that all draws stay strictly positive.

G_cov: float = 0.3

G_mean_Pa: float

depth_m: float

sample(rng: Generator) → SoilState

soil_type: str = 'sand'

su_or_phi_cov: float = 0.1

su_or_phi_mean: float = 35.0

op3.uq.propagation.propagate_pisa_mc(*, diameter_m: float, embed_length_m: float, soil_priors: list[SoilPrior], n_samples: int = 500, seed: int = 42, correlated: bool = True) → ndarray

Run an MC sweep through pisa_pile_stiffness_6x6 and return an (n_samples, 6, 6) array of K matrices.

Parameters:

correlated – If True (default), all layers share the same realisation of the underlying standard normal so that “soft” draws have all layers softer simultaneously (perfectly correlated). If False, each layer is sampled independently. The published practice for site-specific design is to assume strong correlation.

op3.uq.propagation.summarise_samples(samples: ndarray) → dict

Reduce an (n, 6, 6) sample array to per-DOF summary statistics on the diagonal terms.

Returns:

Keys: ‘Kxx’, ‘Kyy’, ‘Kzz’, ‘Krxrx’, ‘Kryry’, ‘Krzrz’ Each value is a sub-dict with mean, std, p05, p50, p95.

Return type:

dict

Polynomial Chaos Expansion surrogate (Phase 5 / Task 5.2).

A minimal Hermite-polynomial PCE for surrogate modelling of expensive Op^3 responses (eigenvalue, DLC outputs) as a function of one or two standard-normal input parameters. The surrogate is built by pseudo-spectral projection over a Gauss-Hermite quadrature grid:

f(xi) ~ sum_{k=0}^{p} c_k H_k(xi)

where H_k are the probabilist Hermite polynomials and the coefficients are obtained as

c_k = E[ f * H_k ] / E[ H_k^2 ]

with the expectation taken over the standard-normal density. Hermite polynomial orthogonality gives E[H_k^2] = k!.

For 2D inputs the basis is the tensor product of 1D Hermite polynomials and the quadrature is a tensor product of 1D Gauss-Hermite nodes.

References

Wiener, N. (1938). “The Homogeneous Chaos”. Am. J. Math. 60(4),

897-936.

Xiu, D., & Karniadakis, G. E. (2002). “The Wiener-Askey polynomial

chaos for stochastic differential equations”. SIAM J. Sci. Comput. 24(2), 619-644.

class op3.uq.pce.HermitePCE(coeffs: 'np.ndarray', order: 'int', n_dim: 'int')

Bases: object

coeffs: ndarray

evaluate(xi: float | ndarray, xi2: float | ndarray | None = None) → ndarray

n_dim: int

order: int

op3.uq.pce.build_pce_1d(f: Callable[[float], float], order: int = 4, n_quad: int | None = None) → HermitePCE

Pseudo-spectral 1D Hermite PCE built on Gauss-Hermite quadrature. The input is the standard normal coordinate xi ~ N(0, 1); callers are responsible for mapping back to the physical parameter space.

op3.uq.pce.build_pce_2d(f: Callable[[float, float], float], order: int = 3, n_quad: int | None = None) → HermitePCE

op3.uq.pce.pce_mean_var(pce: HermitePCE) → tuple[float, float]

Closed-form mean and variance from Hermite PCE coefficients. Mean = c_0; Variance = sum_{k>=1} k! * c_k^2 .

op3.uq.pce.pce_sobol_2d(pce: HermitePCE) → dict

First-order and total Sobol indices from a 2D Hermite PCE.

For the tensor-product basis \(\Psi_{ij}(\xi_1, \xi_2) = H_i(\xi_1) H_j(\xi_2)\) the total variance decomposes as

\[V = \sum_{(i,j) \neq (0,0)} i!\,j!\,c_{ij}^2\]

with partial variances grouped by which input is active:

• \(V_1 = \sum_{i \geq 1, j = 0} i!\,c_{i0}^2\)

• \(V_2 = \sum_{i = 0, j \geq 1} j!\,c_{0j}^2\)

• \(V_{12} = \sum_{i \geq 1, j \geq 1} i!\,j!\,c_{ij}^2\)

First-order Sobol: \(S_i = V_i / V\). Total Sobol: \(S_i^T = (V_i + V_{12}) / V\).

Bayesian calibration of tower / foundation parameters (Phase 5 / Task 5.3).

Grid-based Bayesian inversion of a single scalar parameter (typically the tower fore-aft EI scaling factor or the foundation rotational stiffness multiplier) given a measured first natural frequency.

The grid approach is preferred over MCMC for the 1D problems Op^3 needs because:

1. The posterior is uni-modal and concentrated, so 200-500 grid points are enough for sub-percent posterior summaries.

2. No tuning parameters (chain length, burn-in, proposal width).

3. The forward model is fast (~10 ms per Op^3 eigen call) so the 200-point grid runs in 2 seconds end-to-end.

4. The result is fully reproducible – no random seed dependence.

The output is a posterior PDF on the calibration parameter, plus the posterior mean, std, and 5/50/95 percentiles.

class op3.uq.bayesian.BayesianPosterior(grid: 'np.ndarray', prior: 'np.ndarray', likelihood: 'np.ndarray', posterior: 'np.ndarray', mean: 'float', std: 'float', p05: 'float', p50: 'float', p95: 'float')

Bases: object

grid: ndarray

likelihood: ndarray

mean: float

p05: float

p50: float

p95: float

posterior: ndarray

prior: ndarray

std: float

op3.uq.bayesian.grid_bayesian_calibration(*, forward_model: Callable[[float], float], likelihood_fn: Callable[[float], float], grid: ndarray, prior: ndarray | None = None) → BayesianPosterior

1D grid Bayesian calibration.

Parameters:

• forward_model – Map from the calibration parameter (e.g. EI scale factor) to the predicted observable (e.g. f1 in Hz).

• likelihood_fn – Maps a predicted value to its likelihood under the measurement noise model. Use normal_likelihood(measured, sigma).

• grid – Discretisation of the parameter axis.

• prior – Prior PDF over the same grid (un-normalised is fine). If None, a uniform prior is used.

op3.uq.bayesian.normal_likelihood(measured: float, sigma: float) → Callable[[float], float]

Return L(predicted) = N(measured | predicted, sigma^2) callable.

Advanced UQ

Sequential Bayesian updating for scour depth estimation.

Propagates the posterior from one monitoring epoch to the next, so accumulated evidence progressively tightens the diagnosis. The posterior from epoch t becomes the prior for epoch t+1.

Usage

from op3.uq.sequential_bayesian import SequentialBayesianTracker

tracker = SequentialBayesianTracker(grid_resolution=200)

# Epoch 1: first field measurement tracker.update(freq_ratio=0.994, capacity_ratio=0.99, anomaly=False) print(tracker.summary())

# Epoch 2: six months later tracker.update(freq_ratio=0.991, capacity_ratio=0.97, anomaly=False) print(tracker.summary())

# Epoch 3: anomaly detected tracker.update(freq_ratio=0.985, capacity_ratio=0.92, anomaly=True) print(tracker.summary())

# Full trajectory trajectory = tracker.trajectory()

class op3.uq.sequential_bayesian.EpochResult(epoch: int, freq_ratio: float, capacity_ratio: float, anomaly: bool, posterior_mean: float, posterior_std: float, p05: float, p50: float, p95: float, recommended_action: str)

Bases: object

Diagnostic result at a single monitoring epoch.

anomaly: bool

capacity_ratio: float

epoch: int

freq_ratio: float

p05: float

p50: float

p95: float

posterior_mean: float

posterior_std: float

recommended_action: str

class op3.uq.sequential_bayesian.SequentialBayesianTracker(grid_resolution: int = 200, freq_sensitivity: float = 0.059, freq_exponent: float = 1.5, freq_sigma: float = 0.015, capacity_sigma: float = 0.05, anomaly_threshold: float = 0.39, anomaly_sharpness: float = 30.0)

Bases: object

Tracks the scour posterior across monitoring epochs.

At each epoch the previous posterior becomes the new prior, and the incoming sensor observations sharpen the estimate. The trajectory of posterior means over time reveals the degradation rate and enables remaining-useful-life estimation.

Parameters:

• grid_resolution (int) – Number of discrete scour-depth bins from 0 to 1 (S/D).

• freq_sensitivity (float) – Power-law coefficient for the frequency observation model.

• freq_exponent (float) – Power-law exponent for the frequency model.

• freq_sigma (float) – Gaussian likelihood width for the frequency channel.

• capacity_sigma (float) – Gaussian likelihood width for the capacity channel.

• anomaly_threshold (float) – Scour depth (S/D) above which the anomaly detector fires.

• anomaly_sharpness (float) – Steepness of the logistic step function for the anomaly channel.

property current_posterior: ndarray

reset() → None

Reset to uniform prior and clear history.

save(path: str | Path) → None

Save the full tracker state (trajectory + current posterior).

summary() → dict

Return the latest epoch’s diagnostic summary.

trajectory() → List[dict]

Return the full trajectory of posterior means over epochs.

update(freq_ratio: float, capacity_ratio: float, anomaly: bool) → EpochResult

Process one monitoring epoch and update the posterior.

The current posterior (from the previous epoch) is used as the prior, multiplied by the three channel likelihoods, and renormalised. The result becomes the prior for the next epoch.

Chapter 8 encoder <-> Op^3 UQ bridge (Task 20).

Exposes the 1,794-row OptumGX Monte Carlo database used by the dissertation Chapter 8 digital-twin encoder as a first-class Op^3 UQ data source. Previously the MC database was consumed via a side-channel CSV loaded directly by the encoder training script; this module turns it into a propagator that any downstream Op^3 stage (Bayesian calibration, PCE surrogate, DLC sensitivity) can consume.

Real CSV schema (SiteA integrated_database_1794.csv):

run, S_D, scour_m, su0, k_su, Hmax_kN, H_ratio, V_ratio, f1_Hz, f1_f0, fixity_proxy

1794 Monte Carlo runs sampled over:

scour_m in [0.0, 4.0] m su0 in [7.5, 27.8] kPa (surface undrained shear strength) k_su in [12, 35] kPa/m (depth gradient)

outputs:

f1_Hz first fore-aft bending frequency f1_f0 f1 normalised by the pristine (s=0) case Hmax_kN ultimate horizontal capacity at collapse H_ratio Hmax relative to pristine fixity_proxy (base rotational compliance indicator)

Use

from op3.uq.encoder_bridge import load_site_a_mc df = load_site_a_mc() # reads PHD/data/integrated_database_1794.csv # Sample the real joint distribution for Bayesian updates # or encoder training.

prior is a list of SoilPrior objects whose mean and COV are statistically consistent with the database, so propagate_pisa_mc(soil_priors=prior, ...) produces the same joint distribution the encoder was trained on. This decouples the encoder from the raw OptumGX runs and makes every Ch8 derivation traceable through the committed Op^3 pipeline.

op3.uq.encoder_bridge.bayesian_from_encoder(df: DataFrame, *, forward_model, observation_col: str = 'f1_Hz', parameter_col: str = 'G0_top_Pa', sigma: float = 0.005, n_grid: int = 101)

Treat one MC row as the “truth” observation and run an Op^3 Bayesian calibration over another parameter column. Useful for synthetic-truth verification of the encoder: if the encoder is consistent, the posterior mean should recover the row’s true parameter value.

Returns an op3.uq.bayesian.BayesianPosterior.

op3.uq.encoder_bridge.encoder_as_prior(df: DataFrame, columns: list[str], *, default_soil_type: str = 'sand', su_or_phi_default: float = 35.0) → list[SoilPrior]

Turn per-column statistics from the MC database into a list of SoilPrior objects, one per column. Each prior carries the empirical mean and COV computed from the column; the soil_type and strength default are user-selectable.

The resulting prior list is suitable for feeding into op3.uq.propagation.propagate_pisa_mc to reproduce the encoder training distribution through the deterministic PISA pipeline.

op3.uq.encoder_bridge.load_encoder_mc(csv_path: str | Path) → DataFrame

Load a generic Chapter 8 MC database. For the real SiteA case use load_site_a_mc() which auto-resolves the PHD path.

Returns a pandas DataFrame. The CSV is expected to have one row per MC realisation; any columns beyond the minimum set are preserved for downstream feature engineering.

op3.uq.encoder_bridge.load_site_a_mc() → DataFrame

Load the real 1794-sample SiteA integrated MC database from the PHD SSOT (F:/TREE_OF_THOUGHT/PHD/data/integrated_database_1794.csv).

This is the authoritative training set for the Chapter 8 digital twin encoder. The database spans scour [0, 4] m, su0 [7.5, 27.8] kPa, and k_su [12, 35] kPa/m, and records f1_Hz, Hmax_kN, and fixity_proxy outputs from 1794 real OptumGX + OpenSeesPy runs.

Multi-feature encoder training on real FEM data.

Trains the two-modality residual MLP encoder on the 1,794-sample OptumGX-OpenSeesPy Monte Carlo database. This script replaces any earlier proof-of-concept training that used synthetic or partially fabricated data.

Usage

python -m op3.uq.encoder_training –data PHD/data/integrated_database_1794.csv python -m op3.uq.encoder_training –data PHD/data/integrated_database_1794.csv –epochs 500

The trained model is saved as a PyTorch state dict and can be loaded by the Op³ web application for real-time inference.

class op3.uq.encoder_training.ResidualBlock(*args: Any, **kwargs: Any)

Bases: Module

forward(x)

class op3.uq.encoder_training.TrainingConfig(data_path: 'str' = '', epochs: 'int' = 300, batch_size: 'int' = 64, lr: 'float' = 0.001, weight_decay: 'float' = 0.0001, dropout_prob: 'float' = 0.3, latent_dim: 'int' = 64, hidden_dim: 'int' = 128, train_split: 'float' = 0.8, val_split: 'float' = 0.1, seed: 'int' = 42, output_dir: 'str' = 'models/')

Bases: object

batch_size: int = 64

data_path: str = ''

dropout_prob: float = 0.3

epochs: int = 300

hidden_dim: int = 128

latent_dim: int = 64

lr: float = 0.001

output_dir: str = 'models/'

seed: int = 42

train_split: float = 0.8

val_split: float = 0.1

weight_decay: float = 0.0001

class op3.uq.encoder_training.TwoModalityEncoder(*args: Any, **kwargs: Any)

Bases: Module

Two-modality residual MLP for structural state prediction.

Capacity features (6 dim) and dynamic features (3 dim) are each projected into a shared latent space, concatenated, fused through residual blocks, and mapped to three sensor-observable targets.

forward(x_cap: torch.Tensor, x_dyn: torch.Tensor, training: bool = False) → torch.Tensor

op3.uq.encoder_training.load_and_split(cfg: TrainingConfig) → Tuple

Load the MC database and split into train/val/test.

op3.uq.encoder_training.main()

op3.uq.encoder_training.train(cfg: TrainingConfig) → dict

Train the encoder and return metrics.