Technical reference

Mathematical formulation, unit conventions, coordinate systems, and sign conventions for the Op³ framework. This is the authoritative specification for anyone reviewing the code or reproducing the results in a publication.

1. Unit system 

Op³ uses strict SI units everywhere in public APIs:

Quantity	Unit	Symbol
Length	metre	`m`
Mass	kilogram	`kg`
Time	second	`s`
Force	newton	`N`
Pressure / stress / shear modulus	pascal	`Pa`
Stiffness (translation)	newton / metre	`N/m`
Stiffness (rotation)	newton-metre / radian	`Nm/rad`
Mass density	kilogram / cubic metre	`kg/m^3`
Angle	radian	`rad`

The only exception is frequency, which is reported in Hz (cycles/s) in line with published NREL / DNV / IEC tables. Internal OpenSees eigenvalues are angular frequency squared and are converted via \(f = \sqrt{\lambda}/(2\pi)\).

CSV files carrying geotechnical data use mixed conventions inherited from the OptumGX / SACS export formats:

depth_m – metres (SI)
k_ini_kN_per_m – kilonewtons / metre (multiplied by 1e3 on load)
p_ult_kN_per_m – kilonewtons / metre
D_total_kJ – kilojoules
w_z – dimensionless weight
su – pascal (SI)
phi – degrees (converted to radians internally when needed)

The unit invariance test 2.13 in Verification & Validation verifies that the entire pipeline gives bit-identical frequencies when rebuilt in mm / N / tonne units.

2. Coordinate system 

Op³ uses a right-handed Cartesian coordinate system that matches the OpenFAST v5 / SubDyn convention:

  +z  (upward, from mudline toward hub)
  ^
  |
  +--------> +x  (downwind, nominal wind direction)
  /
 /
+y  (to the left of the rotor when viewed from behind)

The tower base sits on or above the mudline, typically at z = TowerBsHt read from the OpenFAST ElastoDyn file.
The mudline is at z = 0 for monopiles and at z = -WaterDepth for offshore structures measured from MSL.
The hub is at z = TowerHt + Twr2Shft + NacCMzn above the tower base (v0.4.0+ with the rigid CM offset).
Tower foreaft bending is in the x-z plane; side-side is in y-z.

3. Degree-of-freedom numbering 

OpenSees in 3D with -ndf 6 uses this DOF ordering at every node:

DOF 1 : translation in x (Ux)     foreaft tower displacement
DOF 2 : translation in y (Uy)     side-side tower displacement
DOF 3 : translation in z (Uz)     axial
DOF 4 : rotation about x (Rx)     tower side-side bending moment
DOF 5 : rotation about y (Ry)     tower foreaft bending moment
DOF 6 : rotation about z (Rz)     torsion

The 6x6 head-stiffness matrix K is indexed in this order:

          Ux    Uy    Uz    Rx    Ry    Rz
     +---------------------------------------
Ux   | Kxx    0     0     0     K_x_ry  0
Uy   | 0    Kyy    0    K_y_rx  0       0
Uz   | 0     0    Kzz    0      0       0
Rx   | 0   K_y_rx  0   Krxrx    0       0
Ry   | K_x_ry 0    0     0    Kryry     0
Rz   | 0     0     0     0      0      Krzrz

The off-diagonal coupling terms K[0,4] and K[1,3] are the lateral-rocking coupling. For a monopile in the Op³ convention K[0,4] < 0 and K[1,3] > 0.

4. Foundation modes – mathematical definitions 

Mode A (FIXED)

Pure rigid constraint at the tower base node via ops.fix. No spring elements, no compliance. Used as the limiting case for Mode B (K_diag -> infinity).

Mode B (STIFFNESS_6X6)

A lumped 6x6 elastic relation between forces and displacements at the tower base:

\[\begin{split}\begin{bmatrix} F_x \\ F_y \\ F_z \\ M_x \\ M_y \\ M_z \end{bmatrix} = \mathbf{K}_{6\times 6} \begin{bmatrix} u_x \\ u_y \\ u_z \\ \psi_x \\ \psi_y \\ \psi_z \end{bmatrix}\end{split}\]

K must be symmetric and positive-definite (enforced by the V&V test C4). The current Op³ implementation uses an OpenSees zeroLength element with six uniaxialMaterial “Elastic” instances on the diagonal; full off-diagonal coupling via zeroLengthND is a v0.4 extension.

Mode C (DISTRIBUTED_BNWF)

Distributed Winkler springs along the embedded length of the pile. At depth \(z\), the local reaction is:

\[p(z, v) = k(z) \cdot v\]

where k(z) is read from the spring profile CSV. The head stiffness is obtained by Winkler integration:

\[\begin{split}K_{xx} &= \int_0^L k(z)\,dz \\ K_{x,\psi_y} &= \int_0^L k(z) \cdot z\,dz \\ K_{\psi_y\psi_y} &= \int_0^L k(z) \cdot z^2\,dz + \int_0^L k_m(z)\,dz\end{split}\]

where \(k_m(z)\) is the distributed moment stiffness (zero if not provided). Scour relief is applied by multiplying \(k(z)\) by \(\sqrt{(z-s)/z}\) for \(z > s\) (and zero below).

Mode D (DISSIPATION_WEIGHTED)

The novel Op³ contribution. The elastic Mode C springs are multiplied by a dimensionless weighting function:

\[w(D, D_{\max}, \alpha, \beta) = \beta + (1 - \beta) \left(1 - \frac{D}{D_{\max}}\right)^\alpha\]

where \(D(z)\) is the cumulative plastic dissipation from an OptumGX analysis and \(D_{\max}\) is the layer maximum. The weighted stiffness is:

\[k^D_i = k^{\rm el}_i \cdot w(D_i, D_{\max}, \alpha, \beta)\]

Parameters \(\alpha\) (default 1.0) and \(\beta\) (default 0.05) are the only free knobs. See docs/MODE_D_DISSIPATION_WEIGHTED.md for the full formulation, limits, and validation gates.

5. PISA conic shape function 

The PISA framework (Burd 2020 / Byrne 2020) represents soil reactions via a 4-parameter conic curve:

\[\frac{y}{y_u} = \frac{1 + n(1 - x/x_u) - \sqrt{\left(1 + n(1 - x/x_u)\right)^2 - 4\,n\,\frac{x}{x_u}\,(1-n)}}{2(1-n)}\]

with four component-specific parameters:

\(k\) – dimensionless initial slope
\(n\) – curvature exponent (0 = bilinear, 1 = elastic-perfectly-plastic)
\(x_u\) – ultimate normalised displacement
\(y_u\) – ultimate normalised reaction

In Op³ the four parameters are depth-dependent (v1.0.0-rc1+):

\[\begin{split}k(z/D) &= k_1 + k_2 \cdot (z/D) \\ n(z/D) &= n_1 + n_2 \cdot (z/D) \\ x_u(z/D) &= x_{u,1} + x_{u,2} \cdot (z/D) \\ y_u(z/D) &= y_{u,1} + y_{u,2} \cdot (z/D)\end{split}\]

For base components the variable is \(L/D\) (the full pile slenderness) rather than \(z/D\). Coefficients are pinned to the published calibrations:

Sand – Burd 2020 Table 5 (PISA_SAND in op3/standards/pisa.py)
Clay – Byrne 2020 Table 4 second-stage (PISA_CLAY)

The normalisations (Byrne 2020 Table 1) are:

\[\begin{split}\bar{p} &= \frac{p}{\sigma'_v \cdot D} &\bar{v} &= \frac{v \cdot G}{D \cdot \sigma'_v} \\ \bar{m} &= \frac{m}{\sigma'_v \cdot D^2} &\bar{\psi} &= \frac{\psi \cdot G}{\sigma'_v} \\ \bar{H_B} &= \frac{H_B}{\sigma'_v \cdot D^2} &\bar{M_B} &= \frac{M_B}{\sigma'_v \cdot D^3}\end{split}\]

For sand, \(\sigma'_v\) is the effective vertical stress (approximated as \(\gamma'_{\rm eff} \cdot z\) with a default \(\gamma' = 10 \text{ kN/m}^3\)). For clay the normalising pressure is the undrained shear strength \(s_u\).

6. Effective head stiffness under eccentric load 

The McAdam 2020 / Byrne 2020 field-test k_Hinit metric is the secant slope of H vs ground-level displacement with the load applied at height \(h\) above the seabed. For the 2x2 (x, ry) block of K, the ground-level displacement is:

\[\begin{split}u_x &= \frac{K_{ryry} - h K_{x,ry}}{\det} H \\ \det &= K_{xx} K_{ryry} - K_{x,ry}^2\end{split}\]

so the effective head stiffness is:

\[k_{H,\rm init} = \frac{\det}{K_{ryry} - h \cdot K_{x,ry}}\]

This is implemented in op3.standards.pisa.effective_head_stiffness() and is the correct comparator to the published field-test k_Hinit values.

7. Hardin-Drnevich modulus reduction 

Cyclic soil degradation follows the modified hyperbolic:

\[\frac{G}{G_{\max}} = \frac{1}{1 + (\gamma / \gamma_{\rm ref})^a}\]

with curvature exponent \(a = 1\) (default) or \(a = 0.92\) (Darendeli 2001). At \(\gamma = \gamma_{\rm ref}\) the ratio is exactly 0.5 by definition.

Reference strain is plasticity-index-dependent via Vucetic & Dobry (1991):

\[\begin{split}\gamma_{\rm ref}(PI = 0) &\approx 1 \times 10^{-4} \\ \gamma_{\rm ref}(PI = 200) &\approx 5 \times 10^{-3}\end{split}\]

with linear-in-log interpolation between the knots digitised from Figure 5 of the paper.

8. Tower bending frequency formulas (analytical references)

Op³ is verified against Euler-Bernoulli closed-form cantilever results.

Cantilever without tip mass 

\[f_n = \frac{\beta_n^2}{2\pi L^2} \sqrt{\frac{EI}{m_L}}\]

with \(\beta_1 L = 1.87510\) for the first mode, \(\beta_2 L = 4.69409\) for the second, etc. m_L is the mass per unit length.

Cantilever with tip mass (Rayleigh)

\[f_1 = \frac{1}{2\pi} \sqrt{\frac{3 EI}{L^3 (M_{\rm tip} + 0.2235\, m_L L)}}\]

from Blevins (1979) Table 8-1. The 0.2235 coefficient is the Rayleigh modal-mass fraction for the first mode of a uniform cantilever.

Static tip deflection 

\[\delta = \frac{P L^3}{3 EI}\]

used in the V&V test C for a bit-exact comparison.

9. RNA mass and inertia convention 

The rotor-nacelle assembly (RNA) is placed at a rigid offset from the tower top via ops.rigidLink("beam", ...):

\[\text{offset} = (N_{CM,xn}, 0, \text{Twr2Shft} + N_{CM,zn})\]

read from the OpenFAST ElastoDyn file. The rigid CM node carries:

Translational mass \(m_{RNA} = m_{\rm hub} + m_{\rm nac} + n_{\rm blades} \cdot m_{\rm blade}\)
Roll inertia (x-axis) \(\approx (I_{\rm hub} + I_{\rm nac,yaw})/2\)
Pitch inertia (y-axis) \(= I_{\rm hub}\)
Yaw inertia (z-axis) \(= I_{\rm nac,yaw}\)

10. Hermite polynomial chaos 

The Op³ PCE uses probabilist Hermite polynomials \(He_k(\xi)\) orthogonal under the standard normal density \(\phi(\xi) = (2\pi)^{-1/2} e^{-\xi^2/2}\):

\[\mathbb{E}[He_i(\xi) He_j(\xi)] = \delta_{ij} \cdot i!\]

A 1D PCE of order p is the expansion:

\[f(\xi) \approx \sum_{k=0}^{p} c_k He_k(\xi)\]

with pseudo-spectral projection coefficients:

\[c_k = \frac{1}{k!} \int_{-\infty}^{\infty} f(\xi) He_k(\xi) \phi(\xi)\,d\xi \approx \frac{1}{k!} \sum_i w_i f(\xi_i) He_k(\xi_i)\]

where \((\xi_i, w_i)\) are Gauss-Hermite quadrature nodes and weights. Critical implementation detail: numpy.polynomial.hermite_e.hermegauss returns weights for \(\int f(x) e^{-x^2/2} dx\), missing the \(1/\sqrt{2\pi}\) standard-normal density factor. Op³ divides the weights by \(\sqrt{2\pi}\) internally – see the comment in op3/uq/pce.py.

Mean and variance from the coefficients (closed form, no resampling):

\[\mathbb{E}[f] = c_0 \qquad \text{Var}[f] = \sum_{k=1}^{p} k! \cdot c_k^2\]

Sobol indices (2D):

\[S_1 = V_1 / V, \quad S_2 = V_2 / V, \quad S_{12} = V_{12}/V\]

where \(V_i = \sum_{k_i \geq 1, k_j = 0} i! j! c_{ij}^2\).

11. Grid Bayesian calibration 

For a scalar calibration parameter \(\theta\) with forward model \(g(\theta)\) and measured observable \(y_{\rm meas}\), the posterior is:

\[p(\theta | y) \propto p(y | \theta) \cdot p(\theta)\]

with Gaussian likelihood:

\[p(y | \theta) = \frac{1}{\sigma \sqrt{2\pi}} \exp\left(-\frac{(y_{\rm meas} - g(\theta))^2}{2\sigma^2}\right)\]

Op³ evaluates \(g(\theta)\) on a user-supplied grid, multiplies by the prior, normalises via trapezoidal integration, and reports posterior mean, std, and quantiles. The grid approach is preferred over MCMC for the 1D problems Op³ needs because the posterior is uni-modal and 200-500 points are sufficient for sub-percent accuracy with no tuning.

12. Key references 

Topic	Citation
PISA sand	Burd et al. (2020), Geotechnique 70(11), 1048-1066
PISA clay	Byrne et al. (2020), Geotechnique 70(11), 1030-1047
Dunkirk field tests	McAdam et al. (2020), Geotechnique 70(11), 986-998
Ground characterisation	Zdravkovic et al. (2020), Geotechnique 70(11), 945-962
OC6 Phase II	Bergua et al. (2021), NREL/TP-5000-79989
NREL 5 MW	Jonkman et al. (2009), NREL/TP-500-38060
OC3 monopile	Jonkman & Musial (2010), NREL/TP-500-47535
IEA 15 MW	Gaertner et al. (2020), NREL/TP-5000-75698
Hardin-Drnevich	Hardin & Drnevich (1972), J. Soil Mech. Found. Div. 98(7)
Vucetic-Dobry	Vucetic & Dobry (1991), J. Geotech. Eng. 117(1)
Houlsby-Byrne caisson	Houlsby & Byrne (2005), Proc. ICE Geotech. Eng. 158(2/3)
Hermite PCE	Xiu & Karniadakis (2002), SIAM J. Sci. Comput. 24(2)
Rayleigh cantilever	Blevins (1979), Formulas for Natural Frequency, Table 8-1
Phoon-Kulhawy COVs	Phoon & Kulhawy (1999), Can. Geotech. J. 36(4)

See paper/paper.bib for the full BibTeX entries.