wave-theory
SKILL.md
Wave Theory SME Skill
Comprehensive ocean wave theory expertise including wave mechanics, spectral analysis, wave statistics, and irregular sea modeling for offshore engineering applications.
When to Use This Skill
Use wave theory knowledge when:
- Wave spectra - JONSWAP, Pierson-Moskowitz, scatter diagrams
- Wave statistics - Significant wave height, spectral parameters
- Irregular seas - Generate time series from spectra
- Wave kinematics - Particle velocities and accelerations
- Wave transformation - Shoaling, refraction, diffraction
- Extreme values - Design wave estimation
Core Knowledge Areas
1. Regular Wave Theory
Linear (Airy) Wave Theory:
import numpy as np
def airy_wave_properties(
H: float,
T: float,
d: float,
g: float = 9.81
) -> dict:
"""
Calculate Airy wave properties.
Valid for: H/L < 0.14, d/L > 0.5 (deep water) or d/L < 0.05 (shallow)
Args:
H: Wave height (m)
T: Wave period (s)
d: Water depth (m)
g: Gravity (m/s²)
Returns:
Wave properties dictionary
"""
# Wave frequency
omega = 2 * np.pi / T
# Dispersion relation: ω² = gk·tanh(kd)
# Solve iteratively for wave number k
from scipy.optimize import fsolve
def dispersion(k):
return omega**2 - g * k * np.tanh(k * d)
k0 = omega**2 / g # Deep water approximation
k = fsolve(dispersion, k0)[0]
# Wave length
L = 2 * np.pi / k
# Wave celerity (phase speed)
C = omega / k
# Group velocity
n = 0.5 * (1 + 2*k*d / np.sinh(2*k*d)) # Shoaling coefficient
Cg = n * C
# Deep water classification
if d / L > 0.5:
regime = "Deep water"
elif d / L < 0.05:
regime = "Shallow water"
else:
regime = "Intermediate"
return {
'height_m': H,
'period_s': T,
'depth_m': d,
'wavelength_m': L,
'wave_number': k,
'frequency_rad_s': omega,
'frequency_hz': omega / (2*np.pi),
'celerity_m_s': C,
'group_velocity_m_s': Cg,
'regime': regime,
'd_over_L': d / L,
'H_over_L': H / L,
'steepness': H / L
}
# Example
wave = airy_wave_properties(H=8, T=12, d=1500)
print(f"Wave Properties (H={wave['height_m']}m, T={wave['period_s']}s):")
print(f" Wavelength: {wave['wavelength_m']:.1f} m")
print(f" Regime: {wave['regime']} (d/L = {wave['d_over_L']:.3f})")
print(f" Celerity: {wave['celerity_m_s']:.2f} m/s")
print(f" Steepness: {wave['steepness']:.4f}")
Wave Kinematics:
def wave_particle_kinematics(
z: float,
H: float,
T: float,
d: float,
t: float = 0,
x: float = 0,
g: float = 9.81
) -> dict:
"""
Calculate wave particle velocities and accelerations.
Args:
z: Vertical position (0 at SWL, negative below)
H: Wave height (m)
T: Wave period (s)
d: Water depth (m)
t: Time (s)
x: Horizontal position (m)
g: Gravity (m/s²)
Returns:
Particle kinematics
"""
# Wave properties
wave = airy_wave_properties(H, T, d, g)
k = wave['wave_number']
omega = wave['frequency_rad_s']
# Amplitude
a = H / 2
# Hyperbolic functions
cosh_kz_d = np.cosh(k * (z + d))
sinh_kz_d = np.sinh(k * (z + d))
cosh_kd = np.cosh(k * d)
sinh_kd = np.sinh(k * d)
# Wave phase
phase = k * x - omega * t
# Horizontal velocity
u = (omega * a * cosh_kz_d / sinh_kd) * np.cos(phase)
# Vertical velocity
w = (omega * a * sinh_kz_d / sinh_kd) * np.sin(phase)
# Horizontal acceleration
ax = -(omega**2 * a * cosh_kz_d / sinh_kd) * np.sin(phase)
# Vertical acceleration
az = (omega**2 * a * sinh_kz_d / sinh_kd) * np.cos(phase)
# Dynamic pressure
p_dynamic = g * a * (cosh_kz_d / cosh_kd) * np.cos(phase)
return {
'horizontal_velocity': u,
'vertical_velocity': w,
'horizontal_acceleration': ax,
'vertical_acceleration': az,
'dynamic_pressure': p_dynamic,
'total_velocity': np.sqrt(u**2 + w**2),
'total_acceleration': np.sqrt(ax**2 + az**2)
}
# Example: Surface velocity (z=0)
kinematics = wave_particle_kinematics(z=0, H=8, T=12, d=1500, t=0, x=0)
print(f"Surface Particle Kinematics:")
print(f" Horizontal velocity: {kinematics['horizontal_velocity']:.2f} m/s")
print(f" Vertical velocity: {kinematics['vertical_velocity']:.2f} m/s")
print(f" Total velocity: {kinematics['total_velocity']:.2f} m/s")
2. Wave Spectra
JONSWAP Spectrum:
def jonswap_spectrum(
frequencies: np.ndarray,
Hs: float,
Tp: float,
gamma: float = 3.3,
alpha: float = None
) -> np.ndarray:
"""
Calculate JONSWAP wave spectrum.
S(f) = α g² (2π)^-4 f^-5 exp[-5/4(f/fp)^-4] γ^exp[-(f-fp)²/(2σ²fp²)]
Args:
frequencies: Frequency array (Hz)
Hs: Significant wave height (m)
Tp: Peak period (s)
gamma: Peak enhancement factor (3.3 for North Sea)
alpha: Phillips constant (calculated if None)
Returns:
Spectral density S(f) (m²/Hz)
"""
g = 9.81
fp = 1 / Tp # Peak frequency (Hz)
# Calculate alpha if not provided
if alpha is None:
# Relationship: Hs = 4*sqrt(m0)
# For JONSWAP: alpha ≈ 5.061 * Hs² / Tp⁴ * (1 - 0.287*ln(γ))
alpha = 5.061 * Hs**2 / Tp**4 * (1 - 0.287 * np.log(gamma))
# Sigma parameter
sigma = np.where(frequencies <= fp, 0.07, 0.09)
# Pierson-Moskowitz spectrum
S_PM = alpha * g**2 * (2*np.pi)**(-4) * frequencies**(-5) * \
np.exp(-1.25 * (frequencies / fp)**(-4))
# Peak enhancement
r = np.exp(-(frequencies - fp)**2 / (2 * sigma**2 * fp**2))
gamma_factor = gamma ** r
# JONSWAP spectrum
S = S_PM * gamma_factor
return S
# Example: Generate JONSWAP spectrum
freq = np.linspace(0.01, 0.5, 500)
S = jonswap_spectrum(freq, Hs=8.5, Tp=12.0, gamma=3.3)
# Verify Hs
m0 = np.trapz(S, freq)
Hs_calc = 4 * np.sqrt(m0)
print(f"JONSWAP Spectrum:")
print(f" Input Hs: 8.5 m")
print(f" Calculated Hs: {Hs_calc:.2f} m")
print(f" Peak frequency: {1/12:.4f} Hz")
Pierson-Moskowitz Spectrum:
def pierson_moskowitz_spectrum(
frequencies: np.ndarray,
Hs: float,
Tp: float = None,
U19_5: float = None
) -> np.ndarray:
"""
Calculate Pierson-Moskowitz spectrum (fully developed sea).
Args:
frequencies: Frequency array (Hz)
Hs: Significant wave height (m)
Tp: Peak period (s) - optional
U19_5: Wind speed at 19.5m height (m/s) - optional
Returns:
Spectral density S(f) (m²/Hz)
"""
g = 9.81
if Tp is not None:
# Use peak period
fp = 1 / Tp
elif U19_5 is not None:
# Calculate from wind speed
fp = 0.877 * g / (2 * np.pi * U19_5)
else:
raise ValueError("Must provide either Tp or U19_5")
# Phillips constant
alpha = 0.0081 # For fully developed seas
# P-M spectrum
S = alpha * g**2 * (2*np.pi)**(-4) * frequencies**(-5) * \
np.exp(-1.25 * (frequencies / fp)**(-4))
return S
# Example
S_PM = pierson_moskowitz_spectrum(freq, Hs=8.5, Tp=12.0)
m0_PM = np.trapz(S_PM, freq)
Hs_PM = 4 * np.sqrt(m0_PM)
print(f"P-M Spectrum Hs: {Hs_PM:.2f} m")
3. Wave Statistics
Spectral Parameters:
def calculate_spectral_parameters(
S: np.ndarray,
frequencies: np.ndarray
) -> dict:
"""
Calculate spectral wave parameters.
Args:
S: Wave spectrum (m²/Hz)
frequencies: Frequency array (Hz)
Returns:
Spectral parameters
"""
# Spectral moments
m0 = np.trapz(S, frequencies)
m1 = np.trapz(S * frequencies, frequencies)
m2 = np.trapz(S * frequencies**2, frequencies)
m4 = np.trapz(S * frequencies**4, frequencies)
# Significant wave height
Hs = 4 * np.sqrt(m0)
# Mean period
Tm01 = m0 / m1
# Zero-crossing period
Tz = np.sqrt(m0 / m2)
# Peak period (from spectrum maximum)
peak_idx = np.argmax(S)
Tp = 1 / frequencies[peak_idx]
# Spectral width
epsilon = np.sqrt(1 - m2**2 / (m0 * m4))
# Wave steepness
k_mean = 2 * np.pi / (9.81 * Tz**2 / (2*np.pi)) # Deep water approx
steepness = k_mean * Hs / 2
return {
'm0': m0,
'm1': m1,
'm2': m2,
'm4': m4,
'Hs': Hs,
'Tp': Tp,
'Tz': Tz,
'Tm01': Tm01,
'spectral_width': epsilon,
'steepness': steepness
}
# Example
params = calculate_spectral_parameters(S, freq)
print(f"Spectral Parameters:")
print(f" Hs: {params['Hs']:.2f} m")
print(f" Tp: {params['Tp']:.2f} s")
print(f" Tz: {params['Tz']:.2f} s")
print(f" Spectral width: {params['spectral_width']:.3f}")
Wave Height Distribution:
def rayleigh_distribution(
H: np.ndarray,
Hs: float
) -> np.ndarray:
"""
Rayleigh distribution for wave heights in irregular seas.
P(H) = probability that wave height exceeds H
Args:
H: Wave height array (m)
Hs: Significant wave height (m)
Returns:
Exceedance probability
"""
# Rayleigh parameter
H_rms = Hs / np.sqrt(2)
# Exceedance probability
P = np.exp(-(H / H_rms)**2)
return P
def significant_wave_statistics(Hs: float) -> dict:
"""
Calculate wave statistics from Hs using Rayleigh distribution.
Args:
Hs: Significant wave height (m)
Returns:
Wave statistics
"""
H_rms = Hs / np.sqrt(2)
# Various statistical wave heights
H_mean = H_rms * np.sqrt(np.pi / 2)
H_1_10 = H_rms * np.sqrt(2 * np.log(10)) # Average of highest 1/10
H_1_100 = H_rms * np.sqrt(2 * np.log(100)) # Average of highest 1/100
H_max_1000 = H_rms * np.sqrt(2 * np.log(1000)) # Most probable max in 1000 waves
return {
'Hs': Hs,
'H_mean': H_mean,
'H_rms': H_rms,
'H_1_10': H_1_10,
'H_1_100': H_1_100,
'H_max_1000': H_max_1000
}
# Example
stats = significant_wave_statistics(Hs=8.5)
print(f"Wave Statistics (Hs = {stats['Hs']} m):")
print(f" Mean height: {stats['H_mean']:.2f} m")
print(f" H_1/10: {stats['H_1_10']:.2f} m")
print(f" H_1/100: {stats['H_1_100']:.2f} m")
print(f" H_max (in 1000 waves): {stats['H_max_1000']:.2f} m")
4. Time Series Generation
Generate Irregular Wave Time Series:
def generate_irregular_wave_time_series(
S: np.ndarray,
frequencies: np.ndarray,
duration: float,
dt: float,
random_seed: int = None
) -> tuple[np.ndarray, np.ndarray]:
"""
Generate irregular wave elevation time series from spectrum.
Args:
S: Wave spectrum (m²/Hz)
frequencies: Frequency array (Hz)
duration: Duration (s)
dt: Time step (s)
random_seed: Random seed for reproducibility
Returns:
(time, elevation) arrays
"""
if random_seed is not None:
np.random.seed(random_seed)
# Time array
time = np.arange(0, duration, dt)
# Initialize elevation
eta = np.zeros_like(time)
# Frequency resolution
df = frequencies[1] - frequencies[0]
# Generate wave components
for i, f in enumerate(frequencies):
if S[i] > 0:
# Amplitude from spectrum
amplitude = np.sqrt(2 * S[i] * df)
# Random phase
phase = np.random.uniform(0, 2*np.pi)
# Wave component
omega = 2 * np.pi * f
eta += amplitude * np.cos(omega * time + phase)
return time, eta
# Example: Generate 1 hour of wave data
t, elevation = generate_irregular_wave_time_series(
S, freq,
duration=3600,
dt=0.1,
random_seed=42
)
# Verify statistics
Hs_timeseries = 4 * np.std(elevation)
print(f"Time Series Statistics:")
print(f" Target Hs: {params['Hs']:.2f} m")
print(f" Generated Hs: {Hs_timeseries:.2f} m")
print(f" Duration: {len(t) * 0.1 / 3600:.2f} hours")
5. Wave Scatter Diagrams
Create Wave Scatter Diagram:
def create_wave_scatter_diagram(
Hs_bins: np.ndarray,
Tp_bins: np.ndarray,
location_data: dict
) -> np.ndarray:
"""
Create wave scatter diagram (probability table).
Args:
Hs_bins: Hs bin edges (m)
Tp_bins: Tp bin edges (s)
location_data: Historical wave data or hindcast
Returns:
Probability matrix (sum = 1.0)
"""
# This is typically based on hindcast data
# Simplified example using lognormal distribution
n_Hs = len(Hs_bins) - 1
n_Tp = len(Tp_bins) - 1
scatter = np.zeros((n_Hs, n_Tp))
# Simplified: Tp roughly proportional to sqrt(Hs)
for i in range(n_Hs):
Hs_mid = (Hs_bins[i] + Hs_bins[i+1]) / 2
# Expected Tp for this Hs (empirical: Tp ≈ 3.6*sqrt(Hs))
Tp_expected = 3.6 * np.sqrt(Hs_mid)
# Hs occurrence (Weibull distribution)
from scipy.stats import weibull_min
p_Hs = weibull_min.pdf(Hs_mid, c=2, scale=2.5)
# Tp distribution given Hs (normal around expected)
from scipy.stats import norm
for j in range(n_Tp):
Tp_mid = (Tp_bins[j] + Tp_bins[j+1]) / 2
p_Tp_given_Hs = norm.pdf(Tp_mid, loc=Tp_expected, scale=1.5)
scatter[i, j] = p_Hs * p_Tp_given_Hs
# Normalize to probabilities
scatter /= scatter.sum()
return scatter
# Example
Hs_bins = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
Tp_bins = np.array([0, 4, 6, 8, 10, 12, 14, 16])
scatter = create_wave_scatter_diagram(Hs_bins, Tp_bins, {})
# Annual hours
annual_hours = scatter * 8760
print(f"Wave Scatter Diagram:")
print(f" Total probability: {scatter.sum():.4f}")
print(f" Most probable sea state: Hs={Hs_bins[np.unravel_index(scatter.argmax(), scatter.shape)[0]]}-{Hs_bins[np.unravel_index(scatter.argmax(), scatter.shape)[0]+1]}m")
6. Extreme Value Analysis
Design Wave from Return Period:
def calculate_extreme_wave_height(
return_period_years: float,
Hs_annual_max: np.ndarray = None,
distribution: str = 'weibull'
) -> dict:
"""
Calculate design wave height for given return period.
Args:
return_period_years: Return period (years)
Hs_annual_max: Array of annual maximum Hs values
distribution: 'weibull' or 'gumbel'
Returns:
Extreme wave height statistics
"""
from scipy.stats import weibull_min, gumbel_r
if Hs_annual_max is None:
# Example data: 25 years of annual maxima
np.random.seed(42)
Hs_annual_max = weibull_min.rvs(c=2, scale=10, size=25)
# Fit distribution
if distribution == 'weibull':
params = weibull_min.fit(Hs_annual_max)
c, loc, scale = params
dist = weibull_min(c, loc, scale)
elif distribution == 'gumbel':
loc, scale = gumbel_r.fit(Hs_annual_max)
dist = gumbel_r(loc, scale)
else:
raise ValueError("Unknown distribution")
# Exceedance probability for return period
exceedance_prob = 1 / return_period_years
# Extreme value
Hs_extreme = dist.ppf(1 - exceedance_prob)
# Confidence intervals (simplified)
Hs_lower = dist.ppf(1 - exceedance_prob - 0.1)
Hs_upper = dist.ppf(1 - exceedance_prob + 0.1)
return {
'return_period_years': return_period_years,
'Hs_extreme': Hs_extreme,
'Hs_lower_bound': Hs_lower,
'Hs_upper_bound': Hs_upper,
'distribution': distribution,
'exceedance_probability': exceedance_prob
}
# Example: 100-year return period
extreme_100yr = calculate_extreme_wave_height(
return_period_years=100,
distribution='weibull'
)
print(f"100-Year Wave:")
print(f" Hs: {extreme_100yr['Hs_extreme']:.2f} m")
print(f" Range: {extreme_100yr['Hs_lower_bound']:.2f} - {extreme_100yr['Hs_upper_bound']:.2f} m")
Complete Examples
Example 1: Complete Wave Analysis
def complete_wave_analysis(
Hs: float,
Tp: float,
depth: float,
duration: float = 3600,
output_dir: str = 'reports/wave_analysis'
) -> dict:
"""
Complete wave analysis: spectrum, time series, statistics.
Args:
Hs: Significant wave height (m)
Tp: Peak period (s)
depth: Water depth (m)
duration: Time series duration (s)
output_dir: Output directory
Returns:
Complete analysis results
"""
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from pathlib import Path
output_path = Path(output_dir)
output_path.mkdir(parents=True, exist_ok=True)
# 1. Generate spectrum
freq = np.linspace(0.01, 0.5, 500)
S = jonswap_spectrum(freq, Hs, Tp)
# 2. Calculate spectral parameters
params = calculate_spectral_parameters(S, freq)
# 3. Generate time series
t, eta = generate_irregular_wave_time_series(S, freq, duration, dt=0.1)
# 4. Wave statistics
wave_stats = significant_wave_statistics(Hs)
# 5. Regular wave properties (using Tp)
regular_wave = airy_wave_properties(Hs, Tp, depth)
# 6. Create visualizations
fig = make_subplots(
rows=2, cols=2,
subplot_titles=(
'JONSWAP Spectrum',
'Wave Elevation Time Series',
'Wave Height Distribution',
'Wave Steepness'
)
)
# Plot 1: Spectrum
fig.add_trace(
go.Scatter(x=freq, y=S, name='S(f)', line=dict(color='blue')),
row=1, col=1
)
# Plot 2: Time series (first 10 minutes)
t_plot = t[:6000]
eta_plot = eta[:6000]
fig.add_trace(
go.Scatter(x=t_plot, y=eta_plot, name='η(t)', line=dict(width=1)),
row=1, col=2
)
# Plot 3: Wave height distribution
H_array = np.linspace(0, Hs*2, 100)
P_exceedance = rayleigh_distribution(H_array, Hs)
fig.add_trace(
go.Scatter(
x=H_array, y=P_exceedance,
name='Rayleigh',
line=dict(color='red')
),
row=2, col=1
)
# Plot 4: Steepness vs frequency
steepness_freq = (2*np.pi*freq)**2 / 9.81 * np.sqrt(S)
fig.add_trace(
go.Scatter(x=freq, y=steepness_freq, name='Steepness'),
row=2, col=2
)
fig.update_layout(height=800, showlegend=True, title_text=f'Wave Analysis (Hs={Hs}m, Tp={Tp}s)')
fig.write_html(output_path / 'wave_analysis.html')
# Export summary
summary = {
'input': {
'Hs': Hs,
'Tp': Tp,
'depth': depth
},
'spectral_params': params,
'statistics': wave_stats,
'regular_wave': regular_wave,
'time_series': {
'duration_s': duration,
'timestep_s': 0.1,
'points': len(t)
}
}
import json
with open(output_path / 'wave_summary.json', 'w') as f:
json.dump(summary, f, indent=2, default=str)
print(f"✓ Wave analysis complete")
print(f" Output: {output_dir}")
return summary
# Example
analysis = complete_wave_analysis(
Hs=8.5,
Tp=12.0,
depth=1500,
duration=3600
)
Resources
- Shore Protection Manual: US Army Corps of Engineers
- Ocean Waves and Oscillating Systems: J. Falnes
- Water Wave Mechanics for Engineers and Scientists: R.G. Dean & R.A. Dalrymple
- DNV-RP-C205: Environmental Conditions and Environmental Loads
- ISO 19901-1: Metocean design and operating considerations
Use this skill for all wave analysis in DigitalModel!
Weekly Installs
10
Repository
vamseeachanta/workspace-hubFirst Seen
Jan 24, 2026
Security Audits
Installed on
gemini-cli8
claude-code8
codex8
opencode8
antigravity7
windsurf7