简介：面向海洋工程方向学习者与结构设计人员的随机波浪计算代码包，围绕JONSWAP谱随机波浪模拟、小振幅理论下的波浪速度场求解，以及Morison方程计算小尺寸结构物所受波浪力三条主线展开，可服务于海洋平台、立管等结构的耐波性分析入门练习与课程作业验证。包内共2个文件，以1个MATLAB脚本和1张结果说明图为主，压缩包约39KB；脚本对应随机波列生成、速度场计算与波浪力求解流程，图片用于辅助核对计算过程与输出结果。目前已有880人学习下载。读者可按JONSWAP谱参数调整海况并复现算例，理解Morison方程中速度项与加速度项的作用，也可在脚本基础上改写并用于课题验证，对照图示检查结果、降低调试成本。

1. 随机波浪速度及波浪力计算：先把“谱—速度—力”这条链跑通

做随机波浪模拟时，最容易翻车的不是谱公式，而是从随机波浪速度到随机波浪力这条链上的单位、相位和分段。这套“随机波浪速度及波浪力计算”资源，适合手里已经有设计波高、谱峰周期和水深，需要快速生成水面时间序列、节点速度/加速度，并进一步用 Morison 方程估算小尺度构件波浪力的场景。它不解决选型，但能把“谱参数怎么落到时间步、速度场怎么对深度、力怎么积分”这条线跑通。新手可以照着复现基准算例，熟手重点看参数边界和避坑章节。

2. 随机波浪模拟：JONSWAP 谱离散成 2048 个频点的参数表

2.1 为什么随机波浪不能只给一个波高和周期

单色波只有一个频率，水面是干净的正弦，速度、加速度和 Morison 力都围绕同一个相位转。真实海况的能量分布在一段频率上，波面高低不平，拖曳项里的u|u|又是非线性项，单色波会把峰值和相位关系算得很“漂亮”，但和随机波工况差得很远。随机波浪模拟的常见做法是给一条谱，把谱离散成有限个谐波，每个谐波随机相位，再叠加成时间序列。

选谱先看数据条件。只有有效波高Hs和谱峰周期Tp时，JONSWAP 是常用入口，谱峰增强因子gamma常取 3.3；如果资源里给的是 P-M 谱，等价于把谱峰增强压平，不要拿 JONSWAP 的gamma去套 P-M 的公式。Hs和谱零阶矩的关系是Hs = 4*sqrt(m0)，其中m0 = ∫S(ω)dω。很多算例对不上，不是谱公式错，而是把Hs当成了振幅，或者积分变量用了f却套了ω的谱密度。

2.2 频率划分、相位随机化与时间窗

频率下限不能贴零。低频分量决定长波能量，上限决定加速度和拖曳项的高频细节。我一般取ω_min = 0.5ω_p、ω_max = 4ω_p，频点数N=2048，这样谱峰附近有足够采样，尾部也不至于被硬截断得太突然。N不是越大越好，频点过密会让相邻谐波周期太接近，时间序列看起来更“随机”，但计算量也上去了；对常规Tp=8~12 s的工况，2048 点够用。

相位用[0, 2π)均匀随机，但必须固定随机种子。否则每次跑出来的Hs、F_max都有小差异，排查时分不清是参数改了还是随机相位换了。时间窗也要够长，1200 s只能做快速检查，做极值统计要跑多种子，别拿单次F_max当设计值。采样频率fs不能只按最高谐波频率的两倍来卡，Morison 拖曳项是非线性的，fs=20 Hz对Tp=9 s的算例比较稳妥，后面还要做时间步减半检查。

2.3 Python 实现：从谱密度到水面 elevation

下面这段是谱生成和波面叠加的最小闭环。S先按未归一化形状算，再用目标Hs缩放，这样换gamma或频率范围时，波高不会跟着跑偏。

import numpy as np g = 9.81 Hs_target = 3.0 # 目标有效波高，m Tp = 9.0 # 谱峰周期，s gamma = 3.3 # JONSWAP 谱峰增强因子 h = 20.0 # 水深，m N = 2048 # 频点数 wp = 2.0 * np.pi / Tp omega = np.linspace(0.5 * wp, 4.0 * wp, N, endpoint=False) dw = omega[1] - omega[0] # JONSWAP 谱形，先不归一化 sigma = np.where(omega <= wp, 0.07, 0.09) S = (g**2 / omega**5) * np.exp(-1.25 * (wp / omega)**4) S *= gamma ** np.exp(-((omega - wp)**2) / (2.0 * sigma**2 * wp**2)) # 用目标 Hs 反算缩放系数：Hs = 4*sqrt(m0) m0 = np.sum(S) * dw S *= (Hs_target / 4.0)**2 / m0 # 每个谐波振幅：a_i = sqrt(2*S(ω_i)*dω) a = np.sqrt(2.0 * S * dw) rng = np.random.default_rng(202405) phi = rng.uniform(0.0, 2.0 * np.pi, N) T_total = 1200.0 fs = 20.0 t = np.arange(0.0, T_total, 1.0 / fs) eta = np.zeros_like(t) for ai, wi, phii in zip(a, omega, phi): eta += ai * np.cos(wi * t + phii) print("Hs_target =", Hs_target) print("Hs_est =", 4.0 * eta.std())

代码里S的单位是m²/(rad/s)，所以积分用dw，不能用df。a的单位是m，eta是自由水面相对静水位的抬升。phi固定种子是为了复现；T_total和fs决定时间序列长度和分辨率，fs=20 Hz对应dt=0.05 s，对Tp=9 s的随机波浪和后续拖曳力积分都够用。最后打印的Hs_est用4*std(eta)估算，和Hs_target差几个百分点正常，差 20% 以上先查谱积分和缩放。

2.4 谱形校验：有效波高和谱峰周期怎么反算

不要只看波面图“像不像海浪”。谱矩反算是最便宜的校验：m0 = ΣS·dw，m2 = ΣS·ω²·dw，Hs = 4√m0，Tm02 = 2π√(m0/m2)。Tm02不是Tp，但它对频率分布很敏感，能抓出“谱峰位置对了、尾部却错了”的情况。时间序列这边再看4*std(eta)和Tm02的估计值，和谱矩结果对得上，才继续算速度。

如果Hs对、Tm02明显偏小，常见原因是高频尾部截得太早，或者N太小导致谱峰附近欠采样。把ω_max从4ω_p放到5ω_p再看，若Tm02明显变化，说明原来的频率范围不够。若Hs本身不对，先看S缩放前m0是否为正、omega是否从零开始附近有异常值。

3. 随机波浪速度与加速度：有限水深下每个节点的水平、垂直分量

3.1 线性波理论里的速度势与双曲函数

随机波浪速度按线性波叠加：每个谐波有波数k_i，由色散关系ω_i² = g k_i tanh(k_i h)解出。坐标取z=0为静水位，z=-h为海底，x为传播方向，相位θ_i = k_i x - ω_i t + φ_i。波面η = Σa_i cosθ_i，水平速度u = Σa_iω_i cosh[k_i(z+h)]/sinh(k_i h) cosθ_i，垂直速度w = Σa_iω_i sinh[k_i(z+h)]/sinh(k_i h) sinθ_i。加速度是速度对时间求导，水平加速度ax = Σa_iω_i² cosh[k_i(z+h)]/sinh(k_i h) sinθ_i，垂直加速度az = -Σa_iω_i² sinh[k_i(z+h)]/sinh(k_i h) cosθ_i。

这里最容易错的是sinh和cosh的位置。水平速度在海底不为零，垂直速度在海底必须为零；如果你在z=-h处看到w明显不为零，基本就是双曲函数写反了。另一个边界是水面：z=0处水平速度达到该谐波的最大值，垂直速度也达到最大值但相位差 90 度。先把这两个边界试通，再谈随机叠加。

3.2 采样频率、水深和坐标方向的硬约束

色散关系要用全水深形式，不要提前判深水。h/L不大时tanh(kh)不能近似为 1，否则长波速度和加速度都会偏大。波数求解用牛顿迭代，初值取k0=ω²/g，一般 20~30 次就收敛。频率点里不能有ω=0，否则ω²/g和1/sinh(kh)都会出问题；上面的谱从0.5ω_p开始，已经避开了零频。

坐标方向要写死在注释里：z向上为正，海底是-h，静水位是0。如果资源里某个函数把z定义成从海底向上的正坐标，调用前必须换算成z-h，否则cosh(k(z+h))会变成cosh(kz)，越往深处速度越大，属于典型的“公式没错、坐标错了”。采样上，fs=20 Hz能覆盖Tp=9 s的常规工况，但如果你的Tp很小、或者构件直径大导致拖曳项更尖，先做fs=20/40/80 Hz三档对比，峰值差超过 5% 就说明时间步还不够。

3.3 批量计算代码：从 eta 到 u、w、ax、az

下面把波数求解和节点运动学写成函数。输出形状是(nz, nt)，第一维是深度节点，第二维是时间，方便后面按深度积分。

import numpy as np def solve_k(omega, h, g=9.81): """用牛顿迭代解色散关系：omega^2 = g*k*tanh(k*h)""" k = omega**2 / g for _ in range(30): f = g * k * np.tanh(k * h) - omega**2 df = g * (np.tanh(k * h) + k * h / np.cosh(k * h)**2) k -= f / df return k def wave_kinematics(t, a, omega, phi, h, z, x=0.0): """返回 u, w, ax, az，形状均为 (nz, nt)""" k = solve_k(omega, h) theta = k[:, None] * x - omega[:, None] * t[None, :] + phi[:, None] kh = k * h u = np.zeros((len(z), len(t))) w = np.zeros_like(u) ax = np.zeros_like(u) az = np.zeros_like(u) for iz, zi in enumerate(z): cz = np.cosh(k * (zi + h)) / np.sinh(kh) sz = np.sinh(k * (zi + h)) / np.sinh(kh) u[iz] = np.sum(a[:, None] * omega[:, None] * cz[:, None] * np.cos(theta), axis=0) w[iz] = np.sum(a[:, None] * omega[:, None] * sz[:, None] * np.sin(theta), axis=0) ax[iz] = np.sum(a[:, None] * omega[:, None]**2 * cz[:, None] * np.sin(theta), axis=0) az[iz] = -np.sum(a[:, None] * omega[:, None]**2 * sz[:, None] * np.cos(theta), axis=0) return u, w, ax, az z = np.linspace(-h, 0.0, 21) u, w, ax, az = wave_kinematics(t, a, omega, phi, h, z) print(u.shape, w.shape) # (21, 24000)

solve_k的初值用omega**2/g，迭代次数固定 30 次，常规水深下足够；如果h很小、ω很大，tanh饱和得快，也能收敛。wave_kinematics里theta的形状是(Nf, Nt)，每次深度循环只算该层的cz、sz，再对频率维求和。z从-h到0，共 21 个节点，u.shape是(21, 24000)。如果只要某一层，把z改成单元素数组即可。x=0表示在原点看时间序列；沿传播方向布置多个点，就改x并保留同一组phi。

3.4 输出什么：时间序列、包络和极值表

算完不要只留一个u数组。至少输出三类结果：指定深度的u/w/ax/az时间序列，用sqrt(u²+w²)算的速度包络，以及按max/min/mean/std汇总的极值表。Morison 力对u|u|和ax的相位很敏感，包络能快速看出哪个时间段速度最大，极值表则用来和后续力峰值对齐。若资源里有多个构件节点，建议把深度、u_rms、ax_rms、u_max、ax_max放在一张表里，后面查“力为什么大”时不用回头重跑。

4. 波浪力计算：Morison

本文还有配套的精品资源，点击获取