Mathematical formulation, unit conventions, coordinate systems, and sign conventions for the Op3 framework. This is the authoritative specification for anyone reviewing the code or reproducing the results in a publication.
Contents
- 1. Unit system
- 2. Coordinate system
- 3. Degree-of-freedom numbering
- 4. Foundation modes -- mathematical definitions
- 5. PISA conic shape function
- 6. Effective head stiffness under eccentric load
- 7. Hardin-Drnevich modulus reduction
- 8. Tower bending frequency formulas (analytical references)
- 9. RNA mass and inertia convention
- 10. Hermite polynomial chaos
- 11. Grid Bayesian calibration
- 12. Key references
Op3 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 / metreD_total_kJ-- kilojoulesw_z-- dimensionless weightsu-- pascal (SI)phi-- degrees (converted to radians internally when needed)
The unit invariance test 2.13 in :doc:`verification` verifies that the entire pipeline gives bit-identical frequencies when rebuilt in mm / N / tonne units.
Op3 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 = TowerBsHtread from the OpenFAST ElastoDyn file. - The mudline is at
z = 0for monopiles and atz = -WaterDepthfor offshore structures measured from MSL. - The hub is at
z = TowerHt + Twr2Shft + NacCMznabove the tower base (v0.4.0+ with the rigid CM offset). - Tower foreaft bending is in the
x-zplane; side-side is iny-z.
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 Op3
convention K[0,4] < 0 and K[1,3] > 0.
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).
A lumped 6x6 elastic relation between forces and displacements at the tower base:
\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}
K must be symmetric and positive-definite (enforced by the V&V
test C4). The current Op3 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.
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:
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
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).
The novel Op3 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.
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 Op3 the four parameters are depth-dependent (v1.0.0-rc1+):
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)
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_SANDinop3/standards/pisa.py) - Clay -- Byrne 2020 Table 4 second-stage (
PISA_CLAY)
The normalisations (Byrne 2020 Table 1) are:
\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}
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.
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:
u_x &= \frac{K_{ryry} - h K_{x,ry}}{\det} H \\
\det &= K_{xx} K_{ryry} - K_{x,ry}^2
so the effective head stiffness is:
k_{H,\rm init} = \frac{\det}{K_{ryry} - h \cdot K_{x,ry}}
This is implemented in :func:`op3.standards.pisa.effective_head_stiffness` and is the correct comparator to the published field-test k_Hinit values.
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):
\gamma_{\rm ref}(PI = 0) &\approx 1 \times 10^{-4} \\
\gamma_{\rm ref}(PI = 200) &\approx 5 \times 10^{-3}
with linear-in-log interpolation between the knots digitised from Figure 5 of the paper.
Op3 is verified against Euler-Bernoulli closed-form cantilever results.
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.
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.
\delta = \frac{P L^3}{3 EI}
used in the V&V test C for a bit-exact comparison.
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}
The Op3 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. Op3
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.
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)
Op3 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 Op3 needs because the posterior is uni-modal and 200-500 points are sufficient for sub-percent accuracy with no tuning.
| 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.