简介：面向海洋工程学习者与结构设计人员的随机波浪计算资料，聚焦如何由Jonswap谱生成符合统计特性的随机波浪序列，并在小振幅理论下求解水质点水平与垂向速度场，再借助Morison方程估算波浪对小尺寸结构物的作用力。压缩包共2个文件、约39KB，包含1个MATLAB脚本和1个PNG图示：脚本对应随机波浪模拟、速度场计算与波浪力求解的完整流程，图示则用于辅助理解计算步骤或核对结果。已有880人学习，适合具备一定水动力基础、希望把随机波浪谱、速度场与Morison方程串联起来的中高级读者。通过这份资料，可快速获得可运行的算法脚本与结果示意，理清从谱模型到波浪荷载的计算链路，并为海洋结构物耐波性分析与结果验证提供参考。
1. 随机波浪速度及波浪力计算：为什么先算谱再算力
随机波浪速度及波浪力计算，最容易出现的误判是把设计波高当成一个规则波，拿单个波去算力。现场海面是一串相位随机、频率混合的波，构件上的速度、加速度和 Morison 力都是随机过程。先把波浪谱离散成一组频率分量，再叠加出波面与水质点速度，最后代进 Morison 方程做时程统计，这条路径比“设计波法”更接近真实荷载。适合做海上结构物、近岸构件、细长桩柱初步动力分析的人；如果构件直径与波长比已经很大，这套细长杆件思路就该停，换绕射模型。下面按“谱—速度—力—统计—避坑—验证”的顺序拆开。
2. 从波浪谱到水质点速度：相位叠加的最小实现与三个参数
2.1 随机波浪不是“一个波”：谱、相位和频率离散各管什么
随机波浪的速度计算，第一步不是写 Morison 方程，而是把海面描述成一个谱。常见做法是采用 JONSWAP 谱或 Pierson-Moskowitz 谱，用有效波高 Hs 和谱峰周期 Tp 控制能量水平和峰值位置，JONSWAP 再用谱峰增强因子 γ 控制谱峰尖锐度。谱面积 m0 与 Hs 的关系是 Hs=4√m0，这个关系决定了你归一化谱的时候不能拍脑袋。相位在 0 到 2π 之间均匀随机，频率按等间隔离散，每个频率分量对应一个振幅 √(2S(f)Δf) 和一个随机相位，波面就是这些余弦分量的叠加。
水质点速度不是直接把波面乘一个系数。线性波理论下，水平速度的传递函数与 cosh(k(z+h))/sinh(kh) 有关，k 由色散关系 ω²=gk tanh(kh) 求出。深水近似 k=ω²/g 在浅水会明显偏大，速度会整体算高，这是后面避坑章要展开的一条。频率离散范围要覆盖主要能量，一般从 0.5/Tp 到 4/Tp，Nf 取 256 到 1024。相位种子只影响单次实现，不影响大量样本的统计量，但单次峰值会跳，所以做极值统计时不能只取一个种子。
| 输入参数 | 作用 | 常见取值/注意点 |
|---|---|---|
| Hs | 控制谱面积，决定波高量级 | 由设计海况给出，单位 m |
| Tp | 控制谱峰频率 fp=1/Tp | 单位 s，决定主要波长 |
| γ | JONSWAP 谱峰增强因子 | 常用 1 到 7，取 3.3 附近较多 |
| h | 水深 | 有限水深必须用完整色散关系 |
| Nf | 频率离散数 | 256～1024，太低会漏掉高频能量 |
| seed | 随机相位种子 | 单次展示可固定，统计要换种子 |
2.2 用 Python 搭一个可复现的随机波面与速度场
下面这段代码把谱归一化、色散关系、随机相位叠加放在一个函数里，输入时间数组、Hs、Tp、水深 h 和高度 z，输出波面 η、水平速度 u 和水平加速度 ax。谱采用形状函数再缩放到 m0=Hs²/16，避免不同 JONSWAP 写法里 α 系数不一致带来的翻车。
import numpy as np def jonswap_shape(f, fp, gamma=3.3): # 未归一化谱形，后面按目标 m0 缩放 sigma = np.where(f <= fp, 0.07, 0.09) r = np.exp(-((f - fp) ** 2) / (2.0 * sigma ** 2 * fp ** 2)) return f ** -5 * np.exp(-1.25 * (fp / f) ** 4) * gamma ** r def dispersion_k(omega, h, g=9.81): # 牛顿迭代求有限水深色散关系 omega^2 = g*k*tanh(k*h) k = omega ** 2 / g # 深水初值 for _ in range(30): th = np.tanh(k * h) fun = g * k * th - omega ** 2 dfun = g * (th + k * h * (1.0 - th ** 2)) k_new = k - fun / dfun if np.max(np.abs(k_new - k)) < 1e-12: k = k_new break k = k_new return k def random_wave_kinematics(t, Hs, Tp, h, z, gamma=3.3, nf=512, seed=42): # t: 时间数组 s；z: 相对静水面高度 m，海底为 -h rng = np.random.default_rng(seed) fp = 1.0 / Tp f = np.linspace(0.5 / Tp, 4.0 / Tp, nf) df = f[1] - f[0] S = jonswap_shape(f, fp, gamma) S = S / np.trapz(S, f) * (Hs ** 2 / 16.0) # 归一化到 m0 = Hs^2/16 amp = np.sqrt(2.0 * S * df) phase = rng.uniform(0.0, 2.0 * np.pi, nf) omega = 2.0 * np.pi * f k = dispersion_k(omega, h) eta = np.zeros_like(t) u = np.zeros_like(t) ax = np.zeros_like(t) for i in range(nf): theta = omega[i] * t + phase[i] eta += amp[i] * np.cos(theta) coef = amp[i] * omega[i] * np.cosh(k[i] * (z + h)) / np.sinh(k[i] * h) u += coef * np.cos(theta) ax += -coef * omega[i] * np.sin(theta) return eta, u, ax # 调用示例：Tp=8 s 时 fmax=0.5 Hz，采样频率至少 1 Hz，这里取 1.25 Hz fs = 10.0 / 8.0 t = np.arange(0.0, 600.0, 1.0 / fs) eta, u, ax = random_wave_kinematics(t, Hs=3.0, Tp=8.0, h=30.0, z=-5.0, seed=7)这段代码的逻辑是：先构造谱形，再按 Hs 缩放谱面积；频率从 0.5/Tp 到 4/Tp 等间隔取，每个频率分量的振幅由单侧谱密度的离散能量决定；相位由种子生成，保证单次结果可复现；色散关系用牛顿迭代求 k，避免深水近似在浅水把速度算大；速度传递函数里 z 必须大于等于 -h，z=-h 时 cosh(0)=1，分母 sinh(kh) 不为零，海底速度仍可算。调用示例里采样频率 fs 取 1.25 Hz，对应 Nyquist 0.625 Hz，高于 fmax=0.5 Hz，避免高频分量混叠。时间长度 600 s 只是示例，做极值统计时要根据目标超越概率增加。
2.3 频率截断、采样频率和相位种子：三个必须写进配置的输入
频率上限 fmax 对速度影响中等，对加速度影响很大，因为加速度传递函数多乘了一个 ω。fmax 取 4/Tp 是常见起点，但如果 Tp 很小或者构件对高频敏感，截断太早会让加速度峰值偏低。反过来，fmax 取太大而采样频率不够，会出现混叠，速度时程看起来“毛刺”很重，Morison 惯性项会被污染。我一般把 fs 设到 8 到 10 倍 fp，并检查 fmax 是否小于 fs/2。
频率下限 fmin 也不能忽略。取 0.5/Tp 是为了避开 0 Hz 附近的奇异性，但长周期能量如果被切掉，波面会缺少慢变分量，速度的低频趋势会变弱。对一般风浪海况，0.5/Tp 够用；如果谱很宽或者有涌浪成分，可以把 fmin 降到 0.3/Tp，同时增加 Nf。Nf 不是越大越好，Nf 大、df 小，单个分量振幅 √(2SΔf) 变小，但分量数多，计算量线性增加。做蒙特卡洛时 Nf 取 512 通常够，做单次高精度展示可以取 1024。
相位种子要写进配置，不能每次调用都随机。单次图看相位种子影响很大，峰值可能差百分之十几；但 200 个种子以上的峰值分布会稳定。我的习惯是：展示用固定 seed，统计用 seed 序列，并在报告里写清楚样本数和种子范围。频率截断、采样频率、相位种子这三个参数不写进配置，后面复现时就是一笔糊涂账。
3. 用 Morison 方程把速度变成波浪力：拖曳项和惯性项怎么分开算
3.1 Morison 方程适用边界：KC 数、D/L 和 Cd/Cm 怎么选
Morison 方程把单位长度波浪力拆成拖曳力和惯性力：拖曳力 0.5ρCdD|u|u，惯性力 ρCm(πD²/4)ax。ρ 是水密度，D 是构件直径，Cd 是拖曳力系数，Cm 是惯性力系数，u 是水平速度，ax 是水平加速度。这个公式适合细长构件，常见判断是 D/L 小于 0.2 左右，L 用主要波对应的波长；D/L 再大，绕射效应不能忽略，Morison 的细长杆件假设就站不住。KC 数定义为 umaxT/D，KC 大表示拖曳力相对重要，KC 小表示惯性力主导。
Cd 和 Cm 不是常数。圆柱在稳态流里 Cd 随雷诺数变，在波浪往复流里还受 KC 数和表面粗糙度影响。初步计算常用 Cd=1.0、Cm=2.0 作为起点，但这不是万能值。表面粗糙、海生物附着会把 Cd 推高，C m 在 KC 较小时也会偏离 2.0。我的做法是：先用一组基准系数算量级，再做 Cd、Cm 的敏感性扫描，至少看 ±20% 变化对峰值力的影响。如果峰值力对系数变化非常敏感，说明选型阶段不能只报一个数，要把系数区间和对应的统计分位数一起给出。
| 参数 | 含义 | 初步取值/注意 |
|---|---|---|
| ρ | 水密度 | 常用 1025 kg/m³ 左右 |
| D | 构件直径 | 单位 m，沿高度可能变化 |
| Cd | 拖曳力系数 | 基准 1.0 附近，粗糙度会推高 |
| Cm | 惯性力系数 | 基准 2.0 附近，KC 小时需复核 |
| KC | umaxT/D | 判断拖曳/惯性主导 |
| D/L | 直径/波长 | 大于 0.2 要警惕绕射 |
3.2 把速度时程代进 Morison 方程：单位长度力和沿水深积分
单位长度力可以直接用上一节的 u 和 ax 算。下面代码先算某高度 z 处的单位长度力，再沿水深分段积分得到整根构件的总力。注意拖曳项里是 u*np.abs(u)，不是 u**2，否则方向信息丢了，力会永远为正，这是血泪经验里最常见的一条。
rho = 1025.0 D = 0.8 Cd = 1.0 Cm = 2.0 def morison_force_per_meter(u, ax, rho=1025.0, D=0.8, Cd=1.0, Cm=2.0): # 单位长度波浪力 N/m，拖曳项保留方向 f_drag = 0.5 * rho * Cd * D * u * np.abs(u) f_inertia = rho * Cm * np.pi * D ** 2 / 4.0 * ax return f_drag + f_inertia # 沿水深分段积分，假设构件从海底 z=-h 到 z=-1.0 m z_nodes = np.linspace(-30.0, -1.0, 31) dz = z_nodes[1] - z_nodes[0] F_total = np.zeros_like(t) for z in z_nodes: _, u_z, ax_z = random_wave_kinematics( t, Hs=3.0, Tp=8.0, h=30.0, z=z, seed=7 ) f_z = morison_force_per_meter(u_z, ax_z, rho=rho, D=D, Cd=Cd, Cm=Cm) F_total += f_z * dz这段代码的逻辑是：每个高度 z 调用同一个 seed，保证相位一致，速度场在空间上是同一组波分量；单位长度力由拖曳项和惯性项相加；沿水深用矩形积分累加，dz 是分段长度。参数说明：D 沿高度不变时可以直接用；如果构件变直径，要把 D 写成 z 的函数，在循环里逐段替换；z 不能低于 -h，否则 cosh 的参数为负，虽然数学上偶函数，但物理上已经出了水底；分段数 31 只是示例，D/L 小、水深不大时够用，做收敛检查可以加密到 61、121 看总力峰值是否稳定。积分得到的是总力 N，如果要弯矩，还要乘力臂再做一次累加。
3.3 拖曳项和惯性项的相位差：最大速度不等于最大力
拖曳项和 u|u| 同相，惯性项和加速度同相，而加速度与速度在线性波里相差约 90 度。结果是总力峰值不会正好出现在最大速度时刻，也不会正好出现在最大加速度时刻，而是两者按幅值加权后的某个中间相位。浅水、大直径、长周期情况下惯性项可能主导；小直径、强流或 KC 数大时拖曳项更突出。做时程图时，我习惯把 u、ax、拖曳力、惯性力、总力画在同一时间轴上，看峰值错位是否合理。如果总力峰值和最大速度完全重合，而惯性项又很大，大概率是 ax 的符号或相位写错了。
| 主导情形 | 判断线索 | 计算上要注意 |
|---|---|---|
| 拖曳主导 | KC 大、D 小、流速强 | u |
| 惯性主导 | KC 小、D 大、周期长 | ax 的相位和符号要核对 |
| 两者相当 | 峰值错位明显 | 分别输出两项再做统计 |
4. 把随机性算成统计量：有效波高、谱峰周期与波浪力极值
4.1 从时程到谱矩：Hs、Tp 和谱宽参数怎么复核
算完波面和力，不能只看一条时程曲线。先复核波面统计：用 Welch 法估计谱，积分得到 m0，再算 Hs=4√m0；谱峰频率 fp 对应 Tp=1/fp。谱矩 m0、m1、m2 还能给出谱宽参数，检查离散谱是否过窄或过宽。下面代码用 scipy.signal.welch 做谱估计，nperseg 取 4096 或时间序列长度的一半，避免分段太短导致频率分辨率不够。
from scipy.signal import welch def wave_stats_from_eta(eta, dt): fs = 1.0 / dt f, S = welch(eta, fs=fs, nperseg=min(4096, len(eta)), detrend='constant') m0 = np.trapz(S, f) m1 = np.trapz(S * f, f) m2 = np.trapz(S * f ** 2, f) Hs_calc = 4.0 * np.sqrt(m0) fp = f[np.argmax(S)] Tp_calc = 1.0 / fp eps = np.sqrt(max(0.0, 1.0 - m1 ** 2 / (m0 * m2))) return Hs_calc, Tp_calc, eps, f, S Hs_calc, Tp_calc, eps, f_spec, S_spec = wave_stats_from_eta(eta, t[1] - t[0])这段代码的逻辑是：Welch 法把时程分段加窗平均，得到单侧谱密度；m0 积分得到总能量，Hs_calc 应该接近输入的 Hs；fp 取谱峰对应频率，Tp_calc 应该接近输入的 Tp；eps 是谱宽参数，越接近 0 谱越窄。参数说明：nperseg 太小，频率分辨率差，Tp 会跳；nperseg 太大，平均次数少，谱线毛刺多。如果 Hs_calc 和输入差超过 5%，先查谱归一化和采样频率，不要急着调 Morison 系数。
4.2 用蒙特卡洛看相位随机性对峰值力的影响
单次时程的峰值力会随相位种子跳，做设计值不能只取一个种子。下面代码跑 200 个种子，每个种子算一条总力时程，取绝对值最大值，最后看均值、标准差和 95% 分位。样本数 200 是起步，做超越概率更低的极值时要增加。
n_real = 200 peak_forces = [] for seed in range(n_real): eta_s, u_s, ax_s = random_wave_kinematics( t, Hs=3.0, Tp=8.0, h=30.0, z=-5.0, seed=seed ) F_s = morison_force_per_meter(u_s, ax_s, rho=rho, D=D, Cd=Cd, Cm=Cm) peak_forces.append(np.max(np.abs(F_s))) peak_forces = np.array(peak_forces) mean_peak = peak_forces.mean() std_peak = peak_forces.std(ddof=1) p95 = np.percentile(peak_forces, 95)这段代码的逻辑是：每个种子生成独立的随机相位，保持谱和传递函数不变，只改变相位实现；对每一条单位长度力时程取绝对值最大；最后统计峰值分布。参数说明：n_real 越大，分位数越稳定；如果只看均值，50 个种子可能够，但看 95% 分位建议 200 以上；峰值力对 z 敏感，这里只取了 z=-5 m，做整根构件要把沿水深积分放进循环。注意不要用同一个种子算不同 z 再拼成总力，那样空间相位不一致，会人为削弱峰值。
4.3 极值统计：样本时长、超越概率和设计值
峰值力的分布常用 Rayleigh 或 Weibull 拟合，但拟合前要看样本是否独立。相邻峰值之间如果间隔太短，样本不独立，极值会被低估。我的做法是：先按零上跨率或峰值间隔做粗筛，再用经验分位数和拟合分布互相校验。样本时长要覆盖足够多的峰，600 s 在 Tp=8 s 时大约 75 个峰，做 95% 分位勉强，做 99% 分位就太短。增加时长比增加种子更有效，因为单条时程内的峰值相关结构更真实。
| 样本数/时长 | 用途 | 注意 |
|---|---|---|
| 50 个种子 | 看均值、粗标准差 | 不看低超越概率分位 |
| 200 个种子 | 95% 分位起步 | 检查峰值直方图 |
| 3600 s 单条 | 极值分布拟合 | 检查峰独立性和非平稳性 |
| 多海况组合 | 设计值 | 每个海况单独统计再组合 |
5. 避坑与排查：随机波浪力计算里最容易翻车的五个地方
5.1 五个高频踩坑记录
现象一：波面看起来正常，速度却整体偏大。原因通常是把深水色散关系 k=ω²/g 用到了浅水，水深 h 较小、kh 不满足深水条件时，深水 k 偏小，cosh(k(z+h))/sinh(kh) 的传递函数会被放大，速度随之偏高。解决方法是老老实实解 ω²=gk tanh(kh)，用牛顿迭代或二分法都行，迭代初值可以用深水 k，但收敛后要检查 kh 是否合理。每次换水深，先把几个代表频率的 k 和波长打印出来，和手算深水值对比，偏差大就说明浅水效应不能忽略。
现象二：波浪力峰值随相位种子大幅跳变。原因多半是样本太短、频率截断不足，或者只取了一个相位实现就下结论。随机相位叠加的单次峰值本来就有离散性，Tp 附近能量集中时更明显。解决方法是固定谱和传递函数，跑 100 到 200 个种子，看峰值均值和分位数的收敛曲线；如果 95% 分位在 100 个种子后还在明显漂，就加到 300 个。展示单次时程可以固定 seed，但报告里要写明这是单次实现，不是设计值。
现象三：拖曳力符号或量级不对。常见错误是把 u*np.abs(u) 写成 u**2，速度反向时力仍然为正，总力时程的负半周被削掉，峰值统计偏小。另一个错误是单位换算，D 用 mm、u 用 m/s，力会差三个量级。解决方法是先取一个简单正弦速度代进去，手算拖曳力峰值，和代码输出对比；再画 u 和 f_drag 的时程，看两者是否同相、反向时是否变号。单位统一用 m、s、kg、N，进公式前把 D 从 mm 换成 m。
现象四：惯性力在浅水段异常大。原因可能是加速度用差分算，采样频率不够，或者对速度做了低通滤波但没同步处理加速度。线性波里 ax 的传递函数多一个 ω，高频噪声会被放大。解决方法是直接用解析传递函数算 ax，不要对速度做数值差分；如果必须差分，采样频率至少取 8 到 10 倍 fp，并在差分前确认没有混叠。画 ax 的频谱，看高频段是否异常抬高，抬高就查 fs 和 fmax。
现象五：谱峰周期对不上。原因可能是频率分辨率太粗，谱峰落在两个频点之间，argmax 取到了旁边；也可能是谱归一化时把 f 轴单位搞错，Hz 和 rad/s 混用。解决方法是把频率分辨率 df 降到 0.005 Hz 以下，或者对谱峰做抛物线插值；检查 ω=2πf 是否在色散关系里用对。Tp 对不上会直接影响波长和 D/L 判断，进而影响 Morison 适用性，所以这个复核不能省。
5.2 排查顺序：从谱到力的四步核对
第一步核对谱：输入 Hs、Tp 和生成的波面反算 Hs、Tp，偏差超过 5% 就先修谱归一化和频率分辨率。第二步核对速度：取一个固定频率分量，手算传递函数和代码输出对比，确认有限水深 k 和 cosh/sinh 没写反。第三步核对力：用正弦 u 和 ax 代 Morison，检查拖曳项方向、惯性项符号和单位。第四步核对统计：换种子看峰值收敛，换样本时长看极值分位是否稳定。四步里任何一步没闭合，不要往下调 Cd、Cm，否则只是把错误拟合得更像。
6. 验证与提速：用谱矩和收敛曲线判断结果能不能信
验证随机波浪力计算，我一般先看三个数：反算 Hs 与输入 Hs 的偏差、反算 Tp 与输入 Tp 的偏差、峰值力随种子数的变异系数。前两个偏差控制在 5% 以内，第三个在 200 个种子后应该趋于平稳。下面这段代码画峰值力的滑动均值和滑动标准差，用来判断样本数够不够。
import numpy as np def convergence_curve(peaks, window=20): peaks = np.asarray(peaks) n = len(peaks) run_mean = np.array([peaks[:i].mean() for i in range(window, n + 1)]) run_std = np.array([peaks[:i].std(ddof=1) for i in range(window, n + 1)]) return run_mean, run_std run_mean, run_std = convergence_curve(peak_forces, window=20) cv_last = run_std[-1] / run_mean[-1]这段代码的逻辑是：从第 window 个样本开始，逐步增加样本算滑动均值和标准差，最后看变异系数 cv_last。参数说明：window 取 20 是经验值，太小曲线抖，太大看不出早期收敛；如果 cv_last 大于 0.1，说明峰值离散还大，要么加种子，要么加样本时长。提速方面，不要一上来就把 nf 和种子数拉满，先用 nf=256、50 个种子跑通流程，再用收敛曲线决定加到多少。沿水深积分时，z 节点从 31 加到 61，如果总力峰值变化小于 2%，就可以停。
| 检查项 | 通过标准 | 不通过先查 |
|---|---|---|
| Hs 反算偏差 | 小于 5% | 谱归一化、采样频率 |
| Tp 反算偏差 | 小于 5% | 频率分辨率、argmax 插值 |
| 峰值力变异系数 | 200 种子后小于 0.1 | 样本数、频率截断 |
| 沿水深积分收敛 | 31→61 节点峰值差小于 2% | dz、z 范围 |
| D/L 适用性 | 小于 0.2 左右 | 是否该换绕射模型 |
我现在的习惯是：每换一个海况，先跑一条 600 s 时程看波形和力量级，再跑 200 个种子看峰值分布，最后用收敛曲线决定要不要加样本。这个顺序能早点发现谱归一化或色散关系的错，避免在 Morison 系数上反复试。希望帮到你。
本文还有配套的精品资源，点击获取
