随机波浪模拟与波浪力计算：JONSWAP谱到Morison方程

简介面向海洋工程方向学习者与结构设计人员的随机波浪计算代码包围绕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频点数N2048这样谱峰附近有足够采样尾部也不至于被硬截断得太突然。N不是越大越好频点过密会让相邻谐波周期太接近时间序列看起来更“随机”但计算量也上去了对常规Tp8~12 s的工况2048 点够用。相位用[0, 2π)均匀随机但必须固定随机种子。否则每次跑出来的Hs、F_max都有小差异排查时分不清是参数改了还是随机相位换了。时间窗也要够长1200 s只能做快速检查做极值统计要跑多种子别拿单次F_max当设计值。采样频率fs不能只按最高谐波频率的两倍来卡Morison 拖曳项是非线性的fs20 Hz对Tp9 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, endpointFalse) 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的单位是meta是自由水面相对静水位的抬升。phi固定种子是为了复现T_total和fs决定时间序列长度和分辨率fs20 Hz对应dt0.05 s对Tp9 s的随机波浪和后续拖曳力积分都够用。最后打印的Hs_est用4*std(eta)估算和Hs_target差几个百分点正常差 20% 以上先查谱积分和缩放。2.4 谱形校验有效波高和谱峰周期怎么反算不要只看波面图“像不像海浪”。谱矩反算是最便宜的校验m0 ΣS·dwm2 ΣS·ω²·dwHs 4√m0Tm02 2π√(m0/m2)。Tm02不是Tp但它对频率分布很敏感能抓出“谱峰位置对了、尾部却错了”的情况。时间序列这边再看4*std(eta)和Tm02的估计值和谱矩结果对得上才继续算速度。检查量公式期望谱零阶矩m0 ΣS·dw与目标Hs对应谱二阶矩m2 ΣS·ω²·dw用于Tm02对照有效波高Hs 4√m0与输入偏差小于 5%时间序列波高4*std(eta)与谱矩结果接近平均跨零周期Tm02 2π√(m0/m2)不要求等于Tp但量级合理如果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)解出。坐标取z0为静水位z-h为海底x为传播方向相位θ_i k_i x - ω_i t φ_i。波面η Σa_i cosθ_i水平速度u Σa_iω_i cosh[k_i(zh)]/sinh(k_i h) cosθ_i垂直速度w Σa_iω_i sinh[k_i(zh)]/sinh(k_i h) sinθ_i。加速度是速度对时间求导水平加速度ax Σa_iω_i² cosh[k_i(zh)]/sinh(k_i h) sinθ_i垂直加速度az -Σa_iω_i² sinh[k_i(zh)]/sinh(k_i h) cosθ_i。这里最容易错的是sinh和cosh的位置。水平速度在海底不为零垂直速度在海底必须为零如果你在z-h处看到w明显不为零基本就是双曲函数写反了。另一个边界是水面z0处水平速度达到该谐波的最大值垂直速度也达到最大值但相位差 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(zh))会变成cosh(kz)越往深处速度越大属于典型的“公式没错、坐标错了”。采样上fs20 Hz能覆盖Tp9 s的常规工况但如果你的Tp很小、或者构件直径大导致拖曳项更尖先做fs20/40/80 Hz三档对比峰值差超过 5% 就说明时间步还不够。3.3 批量计算代码从 eta 到 u、w、ax、az下面把波数求解和节点运动学写成函数。输出形状是(nz, nt)第一维是深度节点第二维是时间方便后面按深度积分。import numpy as np def solve_k(omega, h, g9.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, x0.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), axis0) w[iz] np.sum(a[:, None] * omega[:, None] * sz[:, None] * np.sin(theta), axis0) ax[iz] np.sum(a[:, None] * omega[:, None]**2 * cz[:, None] * np.sin(theta), axis0) az[iz] -np.sum(a[:, None] * omega[:, None]**2 * sz[:, None] * np.cos(theta), axis0) 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改成单元素数组即可。x0表示在原点看时间序列沿传播方向布置多个点就改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放在一张表里后面查“力为什么大”时不用回头重跑。输出维度用途u, w(nz, nt)拖曳项和相对速度ax, az(nz, nt)惯性项速度包络(nz, nt)找峰值时段极值表nz行对比不同深度和工况4. 波浪力计算Morison本文还有配套的精品资源点击获取