简介面向海洋工程与船舶工程方向的学生、研究者及工程技术人员这份资源聚焦随机波浪的数值模拟与结构物受力计算借助Jonswap谱生成符合实际海况统计特性的随机波序列采用小振幅理论推求水粒子速度场并应用Morison方程求解小尺寸构件的波浪力形成从波谱输入到力输出的完整求解链条。资源包为zip压缩格式共2个文件大小仅39KB。其中1个m文件为MATLAB脚本集成了随机波浪序列构造、速度场计算以及波浪力求解的核心程序另1个png为说明图用于展示计算流程与结果对照便于快速理清算法脉络并验证输出。目前已有879人学习下载主体为初涉随机波浪课题的研究生与工程师。借助这些文件可快速获得可运行的随机波浪力计算原型理解谱参数选取、速度分解与力叠加等关键步骤为后续开展更复杂的波流耦合或结构响应分析提供直接参考。1. 随机波浪速度及波浪力计算从谱到力的完整链路随机波浪模拟不是把几个正弦波凑在一起那么轻巧。我拆过不少波浪力计算脚本最后发现一个共性的现象波高和周期都对准了算出来的力却偏出一倍多。问题往往不在Morison方程本身而在一开始的水质点速度建模——高频段能量的截断、随机相位的种子、频谱形状参数每一项都直接影响速度峰值进而作用于拖曳力。这篇笔记把“随机波浪速度及波浪力计算”这套资源的实现链路拆开从JONSWAP谱生成波面时历到用线性波理论算不同水深的速度和加速度再到用Morison方程算柱体波浪力最后列出几个我实际踩过的参数坑。适合做海上风电基础、桥梁墩柱、浮式结构初步设计的工程师新手能顺着代码跑通熟手则可以直接拿参数边界做复核。2. 随机波浪模拟从谱密度到时间序列2.1 为什么用谱叠加稳态随机过程的本质海面波浪是一个随机过程不是单靠设计波高和周期就能描述的。工程上常说的有效波高Hs来自对波面时历的统计而随机波浪模拟的目标就是生成一条和目标频谱一致的波面时历这样后续计算速度和力才有实际意义。频谱描述了能量随频率的分布比如JONSWAP谱、Pierson-Moskowitz谱都是常见形态。反直觉的一点是即使两条时历的Hs完全一样只要频谱形状和相位不同算出的波浪力极值也可能相差很大。原因是拖曳力正比于速度平方而速度的峰顶分布对高频能量和相位分布极其敏感。所以谱叠加不是“随便叠几个波”而是一个需要认真做参数标定的过程。谱叠加法的理论依据是一个零均值的平稳随机过程它的时域实现可以用无穷多个简谐波叠加近似。振幅按频谱密度确定相位取均匀随机数。随着分量数趋近无穷叠加后时历的各阶统计特征趋于目标谱。实际计算中我们只能取有限个频率分量所以分量数的选取直接决定模拟质量。先写一个JONSWAP谱密度函数这是最常用的有限风区谱import numpy as np def jonswap(f, Hs, Tp, gamma3.3): 生成JONSWAP谱的谱密度值 :param f: 频率数组 (Hz) :param Hs: 有效波高 (m) :param Tp: 谱峰周期 (s) :param gamma: 峰增强因子, 默认3.3 :return: 谱密度 S(f) (m^2*s) g 9.81 fp 1.0 / Tp sigma np.where(f fp, 0.07, 0.09) alpha 0.0081 * (g / (2 * np.pi * fp))**4 * Hs**2 exp_a np.exp(-1.25 * (f / fp)**-4) exp_b np.exp(-((f - fp)**2) / (2 * sigma**2 * fp**2)) S alpha * g**2 / (2*np.pi)**4 / f**5 * exp_a * gamma**exp_b return S这段代码有几个关键参数。alpha是Phillips常数这里用经验式由Hs标定但严格来说这个写法算出的谱面积并不一定等于Hs^2/16因为gamma和峰形会改变面积。常见做法是在alpha上再做一次缩放让谱面积精确对应Hs。我会在后面用谱面积修正这一步很关键。sigma是峰形宽度在峰频左侧取0.07右侧取0.09用来描述谱峰的宽度差异。gamma控制谱峰陡峭程度典型在1到7之间取值。默认3.3是JONSWAP标称值适用于有限水深风浪。如果模拟的是成熟涌浪gamma应取1此时谱就退化为Pierson-Moskowitz形式。2.2 随机相位叠加法实现频率离散与随机相位谱是连续函数模拟前要先离散化。常用的做法是把频率范围分成M个区间在每个区间取一个代表频率用该处的谱密度代表整个区间能量。每个频率分量的振幅取sqrt(2*S(f)*df)这个系数来源于谱密度与振幅平方的关系只适用于单边谱。由于代表频率中点固定、相位随机也有人把它叫做随机相位叠加法是实现随机波浪时历最常见的一种。代表频率的位置一般取每个区间中点也可以按能量等分来取。高频端谱密度小但计数密度高低频端谱密度大如果简单用等间距频率容易丢失谱峰附近的分辨率。更稳妥的做法是让频率间距不等但工程上多数代码为了简单还是用等间距所以M要取足够大。下面是一个完整的波面时历生成函数def simulate_eta(freq_range, Hs, Tp, M, duration, fs, gamma3.3, seed42): 用随机相位叠加法生成随机波浪波面时历 :param freq_range: [fmin, fmax] 模拟频率范围 (Hz) :param Hs: 有效波高 (m) :param Tp: 谱峰周期 (s) :param M: 频率分量数 :param duration: 模拟时长 (s) :param fs: 采样频率 (Hz) :param gamma: 峰增强因子 :param seed: 随机种子 :return: t, eta fmin, fmax freq_range df (fmax - fmin) / M f fmin (np.arange(M) 0.5) * df S jonswap(f, Hs, Tp, gamma) m0 np.sum(S) * df Hs_actual 4 * np.sqrt(m0) S S * (Hs / Hs_actual)**2 amp np.sqrt(2 * S * df) rng np.random.default_rng(seed) phase rng.uniform(0, 2*np.pi, M) t np.arange(0, duration, 1/fs) eta np.zeros_like(t) for k in range(M): eta amp[k] * np.cos(2*np.pi*f[k]*t phase[k]) return t, eta代码逻辑分几步先按等分频率生成代表频率数组再算出谱密度做一次能量修正然后乘以sqrt(2df)得到振幅最后逐个叠加余弦波。注意这一版频率点是固定的只有相位随机这是随机相位叠加法如果想要更好的随机性可以在频率区间内再做一次均匀随机抖动也就是频率随机叠加法。两者都能用但频率随机法生成的时历更不容易出现周期包络代价是实测频谱会略有偏差。这里有几个参数直接决定模拟质量。M建议取400到800太低会让波面出现明显的周期性太高则计算量增大。duration和fs决定了总点数模拟时长至少要覆盖一百个波浪周期比如Tp8s时至少要800秒否则极值统计不可靠。seed用于复现结果设计报告里换一个种子复核是常见的习惯。2.3 参数标定JONSWAP谱的\gamma取值与方向扩散gamma不是越陡越好它描述的是风浪成长中的谱峰集中度。对于有限风区gamma大约在3.3附近对于风区很短的年轻风浪gamma可能更高到5左右对于涌浪或衰减后的波浪gamma接近1。我们一般在代码里先固定gamma3.3但实际项目里如果有实测频谱最好用最小二乘法拟合得到gamma标定值。另一个容易被忽略的参数是方向扩散。上面的代码只模拟了单向长峰波。真实海面是短峰波能量分布在多个方向上。如果结构是细长柱体方向的影响还不太大如果是大尺寸结构或浮体方向谱的作用就不可忽视。最简单的方向分布函数是cos^2形式D(theta) (2/pi)*cos^2(theta)把谱密度乘以方向分布后再把每个方向分量叠加。这样得到的波面时历会明显更“乱”统计也更接近真实海面。对于单柱Morison力计算我一般先做单向波确认整个流程正确后再升级到方向谱。因为单向波算出的力通常比多向波略偏大相当于保守估计这在对标物模时需要注意。下面给一组我用过的标定参数供参考参数典型值说明Hs3.0 m有效波高设计波高统计值Tp8.0 s谱峰周期常取谱峰对应的谱密度最高点周期gamma3.3谱峰增强因子涌浪取1fmin0.005 Hz低于此频率的能量通常可忽略但太低会增大计算量fmax0.6 Hz高频截断由实际风浪谱衰减决定M500频率分量数至少400duration1200 s足够覆盖150个波浪周期fs4 Hz采样频率必须高于关注最高频率的2倍这组参数在计算中能保持频谱形状稳定力极值不会因为模拟时长太短而飘得太厉害。如果只是为了演示可以把M和duration降下来但这样得到的结果只能看趋势不能直接作为设计值。3. 水质点速度与加速度线性理论的计算细节3.1 线性波理论的速度公式与水深修正波面模拟完成后下一步是求水质点速度。线性波理论Airy波在中等水深和深水条件下表现稳定也是大多数Morison力计算的基础。对固定水深h、静水面z0、z轴向上水平速度u和垂直速度w的单个分量为u a * ω * cosh(k(zh)) / sinh(kh) * cos(kx - ωt φ) w a * ω * sinh(k(zh)) / sinh(kh) * sin(kx - ωt φ)其中k是波数满足色散关系ω^2 g k tanh(kh)a是每个分量的振幅ω2πf。这个公式里最关键的是cosh和sinh的比值。在深水区kh很大cosh(k(zh))/sinh(kh)约等于exp(kz)速度随水深指数衰减在浅水区这个比值变成一个线性函数速度沿水深更均匀这也是浅水区和深水区速度剖面差异的根源。写代码时不要直接用深水波数也不要忽略色散关系直接取kω^2/g。很多翻车案例都是因为省了这一步结果深水区看似差不多浅水区速度偏大20%以上。牛顿迭代求解色散关系是最快的方式初值取深水波数ω^2/g即可。下面是一个独立的波数求解函数def wave_number(omega, depth, g9.81): 用定点迭代求解色散关系 omega^2 g*k*tanh(k*h) :param omega: 角频率 (rad/s) :param depth: 水深 (m) :return: 波数 k (1/m) k omega**2 / g for _ in range(50): k_new omega**2 / (g * np.tanh(k * depth)) if np.abs(k_new - k) 1e-10: break k k_new return k这个迭代法的收敛速度很快对于水深特别浅的极端情况初值选深水解可能导致迭代震荡这时可以限制每次迭代的步长或者改用二分法。不过常规海洋工程水深下定点迭代已经够用。注意水深h和z的符号在代码里必须统一我习惯把海底定义为-h这样zh在海底为0、水面为h避免符号混乱。3.2 从频域到时间域速度时历的生成上一章生成的波面时历是许多频率分量的叠加速度同样可以按频率叠加只不过每个分量要乘一个传递函数。水平速度对波面振幅的传递函数为H_u(f, z) ω * cosh(k(zh)) / sinh(kh)注意相位关系在波浪行进方向水平速度与波面同相垂向速度与波面差90度。这意味着垂向速度在波峰处为0在波面过零处最大。很多刚开始写代码的人把垂直速度相位搞错会造成垂向力相位完全不对。下面是生成速度与加速度时历的函数def linear_kinematics(t, f, phase, amp, h, z, g9.81): 基于线性理论计算给定水深点的水平、垂向速度与加速度 :param t: 时间数组 :param f: 频率分量数组 (Hz) :param phase: 各分量相位 (rad) :param amp: 各分量振幅 (m) :param h: 水深 (m) :param z: 目标水深坐标 (m, 静水面为0, 向上为正) :return: u, w, du_dt, dw_dt omega 2 * np.pi * f u np.zeros_like(t) w np.zeros_like(t) du_dt np.zeros_like(t) dw_dt np.zeros_like(t) for k in range(len(f)): K wave_number(omega[k], h, g) theta -omega[k] * t phase[k] cosh_ratio np.cosh(K * (z h)) / np.sinh(K * h) sinh_ratio np.sinh(K * (z h)) / np.sinh(K * h) u amp[k] * omega[k] * cosh_ratio * np.cos(theta) w amp[k] * omega[k] * sinh_ratio * np.sin(theta) du_dt amp[k] * omega[k]**2 * cosh_ratio * np.sin(theta) dw_dt amp[k] * omega[k]**2 * sinh_ratio * np.cos(theta) return u, w, du_dt, dw_dt代码里theta取了x0表示作用点位置在原点实际结构坐标可以通过加kx来平移。水平加速度比水平速度相位领先90度所以cos变成sin垂向加速度比垂向速度领先90度所以sin变成cos。加速度在频域相当于乘iω但在时域实现时直接用解析导数更稳。这个函数的时间复杂度是M乘N如果M500、N4000就是两百万次循环Python跑完大约几秒可以接受。但如果要计算整个桩身沿深度几十个点就要用数组化重写避免三层循环。后面第4.3节我会给出沿桩身积分的思路。3.3 流速剖面插值不同水深处的值波浪力沿桩身是变化的不能只取水面一点的速度。通常做法是在z方向取若干代表点比如从海底到海面均匀取20个点每点算一条速度时历然后按Morison方程计算每点单位长度力再沿深度做数值积分。这样做会得到一个总的力时历而不是简单用表面速度乘以面积。近海底处垂直速度因为sinh(0)0而趋于0但水平速度并不为0。势流理论下海底水平速度仍存在实际流体中边界层内速度会衰减但对Morison力影响不大。所以不要直接把海底点的速度设成0除非你的模型里额外加了边界层修正。插值时要注意如果柱体穿过多个水深范围海洋结构高程坐标要统一。一般把海底定义为z-h静水面为z0向上为正。z0以上是自由水面不能随便取点因为线性理论在自由水面上仍然适用但z0处已超出静水面波浪表面实际可能淹没或暴露这时候要结合瞬时波面做湿表面判断。更粗略的做法是只计算到静水面然后用波峰和波谷两条路径包络工程上也能接受。4. 波浪力计算Morison方程与适用范围4.1 Morison方程拆解惯性力拖曳力直径相对波长较小的细长圆柱体波浪场基本不被柱体扰动此时波浪力可以用Morison方程计算。总水平力由两部分组成F F_drag F_inertia其中拖曳力正比于速度绝对值乘速度方向固定惯性力正比于加速度。对于单位长度柱体f(z,t) 0.5 * ρ * C_D * D * |u(z,t)| * u(z,t) ρ * C_M * π * D^2 / 4 * du/dt(z,t)第一项是速度平方项第二项包含了惯性质量和附加质量效应。公式里的D是柱体直径ρ是海水密度C_M等于1加上附加质量系数C_a。注意这个公式在速度方向变化时拖曳力要用|u|u来处理否则力的方向会乱。适用范围通常用D/L和KC数判断。D/L不应太大一般要求D/L0.2否则需要考虑绕射效应改用大尺度结构物上的Froude-Krylov公式。KC数是衡量拖曳力和惯性力相对重要性的无量纲数KC较小时惯性力主导KC很大时拖曳力主导。对于常见的导管架桩腿D通常在1到2米波长为50到100米时D/L在0.02到0.04Morison方程完全够用。下面是一个单位长度力计算函数可以直接接入前面的速度时历def morison_per_length(u, du_dt, D, rho1025.0, Cd0.8, Cm1.8): 计算单位长度柱体的Morison波浪力时历 :param u: 水平速度时历 (m/s) :param du_dt: 水平加速度时历 (m/s^2) :param D: 柱体直径 (m) :param rho: 海水密度 (kg/m^3) :param Cd: 拖曳力系数 :param Cm: 惯性力系数 :return: 单位长度力时历 (N/m) force_drag 0.5 * rho * Cd * D * np.abs(u) * u force_inertia rho * Cm * (np.pi * D**2 / 4.0) * du_dt return force_drag force_inertia4.2 力系数C_D、C_M的选取与敏感性C_D和C_M是Morison方程里最敏感、也最“有黑匣子味”的两个参数。规范上一般给出范围C_D在0.6到1.2之间C_M在1.5到2.0之间。实际取值要看结构形式、表面粗糙度和振动响应。对光滑圆管C_D通常取0.7左右对数生锈或有海生物附着的表面C_D会明显增大甚至达到1.5。C_M则更依赖KC数。在KC数较小时流动分离不明显C_M接近2KC数增大后尾流效应让惯性力系数下降大约在1.3到1.6之间。工程上如果没有明确数据我先取C_D0.8、C_M1.8后面再用不同取值做敏感性包络。这两个系数对力极值的敏感性不同在速度峰值处拖曳力占主导在波面过零处加速度达到峰值惯性力占主导。如果计算出的合力时历在过零点附近出现大尖峰多半是惯性力系数C_M过大或者加速度相位错误。反过来的情况也遇到过拖曳力系数设得太大导致峰值出现在速度最大时整个时历波形看起来特别胖。我一般会做一个系数敏感性表算三组值工况C_DC_M极值变化基准0.81.8—光滑管0.61.6下降约10%左右粗糙管1.22.0上升约15%左右这个表不用做太多能看出力的极值跨度和包络就够了。报告里写清楚系数依据比随便给一个数更禁得起推敲。4.3 波浪力时历与极值统计把单位长度力沿水深积分就得到整根柱体上的水平力。常用梯形积分或辛普森积分。先在不同水深点计算速度然后调用morison_per_length得到每个节点的力最后积分from scipy.integrate import simpson def total_force_on_pile(z_nodes, u_per_node, du_per_node, D, Cd0.8, Cm1.8): 沿水深积分计算整根桩的波浪力时历 :param z_nodes: 水深节点坐标数组 (m) :param u_per_node: 每个节点的水平速度时历, shape (节点数, 时间点数) :param du_per_node: 每个节点的水平加速度时历 :return: 总力时历 (N) f_per_node np.zeros_like(u_per_node) for i in range(len(z_nodes)): f_per_node[i] morison_per_length(u_per_node[i], du_per_node[i], D, CdCd, CmCm) force_total simpson(f_per_node, xz_nodes, axis0) return force_total, f_per_node如果柱体是倾斜的Morison方程还要考虑沿轴向的切向力如果是锥形管或变直径管D要按水深分段取。这段代码先按等直径直立柱处理变直径时只要把D改成一个按深度取值的数组即可。积分之后得到总力时历然后统计最大极值和最小极值。注意极值分布会受模拟时长影响如果只算一条时历最大力可能是这条时历里碰巧出现的一个大波造成的。所以在初步设计里我会用这个极值去做结构截面校核但最终设计一定要做多种子极值分析具体见第5章。5. 避坑指南随机波浪模拟与波浪力计算的五个常见问题5.1 模拟序列不平稳频率截断与能量泄漏现象生成的波面时历前500秒和后500秒的标准差差了20%以上甚至出现明显的周期包络。原因频率分量数M太少或频率下限截断过高让低频率分量的周期长度接近甚至超过模拟时长造成波面在时历两端没有充分采样。另一个可能性是相位种子和固定频率区间的中心重合度太高导致各分量相位差有规律波面产生了伪周期。解决把M提高到500以上频率下限降到0.005Hz甚至更低并检查波面时历的自相关函数。自相关函数应该随滞后时间迅速衰减如果衰减后又回升说明序列里有周期性成分。这时换一个随机种子再试若序列形态变化太大多半还是分量数不够。5.2 水质点速度与加速度的高频抖动现象波面看起来还算光滑但速度时历尤其是靠近水面点的速度像有高频噪声一样抖动剧烈。原因高频分量的波数k很大cosh(k(zh))在z接近0时增长很快导致高频速度传递函数被放大。虽然波面谱本身在高频段已经很小但经过传递函数放大后速度谱的高频尾部仍然有能量表现为速度时历上的小尺度波动。解决在计算速度前对谱做一次高频截断或滤波。通常把模拟使用的fmax作为截断上限实际计算速度时再用一次略低的上限比如f_velocity_max min(fmax, 2*fp)。还可以对速度时历做低通滤波但更好的办法是从频域相位上直接避免高频分量过大。5.3 Morison方程力的尖峰与系数敏感现象总力时历在波面过零附近出现非常窄的双向尖峰看起来和规则波浪力完全不一样。原因过零时加速度最大惯性力也最大。如果加速度项在时域用差分计算高频分量会被差分运算放大产生数值尖刺。规范做法是在频域直接乘iω得到加速度或者用前面的解析导数方法不要用中心差分。解决检查加速度生成方式确保使用传递函数乘以角频率的平方。同时检查C_M取值在KC较大时C_M应适当降低。把C_M从2.0降到1.6尖峰峰值往往能下降15%左右这个变化正常。5.4 浅水区线性理论失效非线性波的取舍现象在浅水或近岸场地按线性理论算出的波浪力比实测波力大很多速度剖面和波面形状也不对波峰处没有尖瘦形状反而很平。原因线性波理论假设波陡很小浅水里波高和水深比值较大时非线性效应变得十分显著自由表面条件需要保留二阶甚至更高阶项。此时波峰处水质点速度并非线性理论给出的余弦分布。解决用相对水深和波陡判断当kh大于π/2且波陡H/L小于1/50线性理论够用当kh小于π/2或波陡大于1/25改用Stokes五阶波或流函数波。随机波浪模拟中可以先把每个频率分量按线性叠加生成波面再对每个大波做瞬时自由表面修正工程上复杂一些。若只是初步筛选可以在线性结果基础上乘一个1.1到1.2的经验放大系数但最终设计不能这么做。5.5 极值统计偏差随机种子与模拟时长现象同一个Hs和Tp只是换了一下随机种子最大波浪力从580kN变成720kN差了20%以上。原因随机波浪时历只是海况的一种实现极值本身就是随机变量。短时历模拟往往只能覆盖有限个波浪周期最大波高没有机会充分显现。模拟时长不够长时极值完全取决于该时段内偶尔出现的特大波。解决把模拟时长提高到至少3小时10800s并生成20到50个独立种子对每个种子提取极值然后用Gumbel分布拟合用外插得到100年重现期极值。这个流程虽然计算量大但比仅仅拉长单次模拟要稳健得多。我一般会把所有种子的极值画成概率图如果线性拟合度好说明极值分布已收敛。6. 验证与进阶用能量谱闭合校验计算结果整个流程跑完后别急着出结果先做一次能量谱闭合校验。方法很简单对生成的波面时历做傅里叶变换用Welch方法估计实测谱然后和输入的目标JONSWAP谱画在同一张图上。两者的面积即0阶矩应该一致Hs误差在1%以内谱峰频率和宽度也应该大致重合。下面是一段校验代码from scipy.signal import welch def check_spectrum(eta, fs, f_input, S_input): 比较生成波面的实测谱与目标谱 freqs, S_est welch(eta, fsfs, npersegmin(len(eta)//4, 4096), noverlaplen(eta)//16) df_est freqs[1] - freqs[0] Hs_est 4 * np.sqrt(np.sum(S_est) * df_est) df_target f_input[1] - f_input[0] Hs_input 4 * np.sqrt(np.sum(S_input) * df_target) print(f目标Hs: {Hs_input:.2f} m, 实测Hs: {Hs_est:.2f} m) return freqs, S_est速度谱的校验也类似。理论上速度谱等于波面谱乘以传递函数的平方S_u(f) S_eta(f) * (ω * cosh(k(zh))/sinh(kh))^2对比实测速度时历的谱和这个理论谱可以快速发现传递函数代码里的乘除错误。我见过一个脚本把cosh和sinh写反了波面谱校验通过了但速度谱完全不对结果力小了一半。从那以后我每次做随机波浪模拟最后都会强制走一遍波面和速度的双重谱闭合校验否则力和弯矩的极值我都不敢对外报。希望这份笔记里的实现链路和参数坑能帮你把随机波浪速度及波浪力计算这条路走顺。本文还有配套的精品资源点击获取 
长河君 · 2026/10/21 8:05:01 · 阅读 1.2k 随机波浪速度与波浪力计算：从JONSWAP谱到Morison方程的完整实现
简介面向海洋工程与船舶工程方向的学生、研究者及工程技术人员这份资源聚焦随机波浪的数值模拟与结构物受力计算借助Jonswap谱生成符合实际海况统计特性的随机波序列采用小振幅理论推求水粒子速度场并应用Morison方程求解小尺寸构件的波浪力形成从波谱输入到力输出的完整求解链条。资源包为zip压缩格式共2个文件大小仅39KB。其中1个m文件为MATLAB脚本集成了随机波浪序列构造、速度场计算以及波浪力求解的核心程序另1个png为说明图用于展示计算流程与结果对照便于快速理清算法脉络并验证输出。目前已有879人学习下载主体为初涉随机波浪课题的研究生与工程师。借助这些文件可快速获得可运行的随机波浪力计算原型理解谱参数选取、速度分解与力叠加等关键步骤为后续开展更复杂的波流耦合或结构响应分析提供直接参考。1. 随机波浪速度及波浪力计算从谱到力的完整链路随机波浪模拟不是把几个正弦波凑在一起那么轻巧。我拆过不少波浪力计算脚本最后发现一个共性的现象波高和周期都对准了算出来的力却偏出一倍多。问题往往不在Morison方程本身而在一开始的水质点速度建模——高频段能量的截断、随机相位的种子、频谱形状参数每一项都直接影响速度峰值进而作用于拖曳力。这篇笔记把“随机波浪速度及波浪力计算”这套资源的实现链路拆开从JONSWAP谱生成波面时历到用线性波理论算不同水深的速度和加速度再到用Morison方程算柱体波浪力最后列出几个我实际踩过的参数坑。适合做海上风电基础、桥梁墩柱、浮式结构初步设计的工程师新手能顺着代码跑通熟手则可以直接拿参数边界做复核。2. 随机波浪模拟从谱密度到时间序列2.1 为什么用谱叠加稳态随机过程的本质海面波浪是一个随机过程不是单靠设计波高和周期就能描述的。工程上常说的有效波高Hs来自对波面时历的统计而随机波浪模拟的目标就是生成一条和目标频谱一致的波面时历这样后续计算速度和力才有实际意义。频谱描述了能量随频率的分布比如JONSWAP谱、Pierson-Moskowitz谱都是常见形态。反直觉的一点是即使两条时历的Hs完全一样只要频谱形状和相位不同算出的波浪力极值也可能相差很大。原因是拖曳力正比于速度平方而速度的峰顶分布对高频能量和相位分布极其敏感。所以谱叠加不是“随便叠几个波”而是一个需要认真做参数标定的过程。谱叠加法的理论依据是一个零均值的平稳随机过程它的时域实现可以用无穷多个简谐波叠加近似。振幅按频谱密度确定相位取均匀随机数。随着分量数趋近无穷叠加后时历的各阶统计特征趋于目标谱。实际计算中我们只能取有限个频率分量所以分量数的选取直接决定模拟质量。先写一个JONSWAP谱密度函数这是最常用的有限风区谱import numpy as np def jonswap(f, Hs, Tp, gamma3.3): 生成JONSWAP谱的谱密度值 :param f: 频率数组 (Hz) :param Hs: 有效波高 (m) :param Tp: 谱峰周期 (s) :param gamma: 峰增强因子, 默认3.3 :return: 谱密度 S(f) (m^2*s) g 9.81 fp 1.0 / Tp sigma np.where(f fp, 0.07, 0.09) alpha 0.0081 * (g / (2 * np.pi * fp))**4 * Hs**2 exp_a np.exp(-1.25 * (f / fp)**-4) exp_b np.exp(-((f - fp)**2) / (2 * sigma**2 * fp**2)) S alpha * g**2 / (2*np.pi)**4 / f**5 * exp_a * gamma**exp_b return S这段代码有几个关键参数。alpha是Phillips常数这里用经验式由Hs标定但严格来说这个写法算出的谱面积并不一定等于Hs^2/16因为gamma和峰形会改变面积。常见做法是在alpha上再做一次缩放让谱面积精确对应Hs。我会在后面用谱面积修正这一步很关键。sigma是峰形宽度在峰频左侧取0.07右侧取0.09用来描述谱峰的宽度差异。gamma控制谱峰陡峭程度典型在1到7之间取值。默认3.3是JONSWAP标称值适用于有限水深风浪。如果模拟的是成熟涌浪gamma应取1此时谱就退化为Pierson-Moskowitz形式。2.2 随机相位叠加法实现频率离散与随机相位谱是连续函数模拟前要先离散化。常用的做法是把频率范围分成M个区间在每个区间取一个代表频率用该处的谱密度代表整个区间能量。每个频率分量的振幅取sqrt(2*S(f)*df)这个系数来源于谱密度与振幅平方的关系只适用于单边谱。由于代表频率中点固定、相位随机也有人把它叫做随机相位叠加法是实现随机波浪时历最常见的一种。代表频率的位置一般取每个区间中点也可以按能量等分来取。高频端谱密度小但计数密度高低频端谱密度大如果简单用等间距频率容易丢失谱峰附近的分辨率。更稳妥的做法是让频率间距不等但工程上多数代码为了简单还是用等间距所以M要取足够大。下面是一个完整的波面时历生成函数def simulate_eta(freq_range, Hs, Tp, M, duration, fs, gamma3.3, seed42): 用随机相位叠加法生成随机波浪波面时历 :param freq_range: [fmin, fmax] 模拟频率范围 (Hz) :param Hs: 有效波高 (m) :param Tp: 谱峰周期 (s) :param M: 频率分量数 :param duration: 模拟时长 (s) :param fs: 采样频率 (Hz) :param gamma: 峰增强因子 :param seed: 随机种子 :return: t, eta fmin, fmax freq_range df (fmax - fmin) / M f fmin (np.arange(M) 0.5) * df S jonswap(f, Hs, Tp, gamma) m0 np.sum(S) * df Hs_actual 4 * np.sqrt(m0) S S * (Hs / Hs_actual)**2 amp np.sqrt(2 * S * df) rng np.random.default_rng(seed) phase rng.uniform(0, 2*np.pi, M) t np.arange(0, duration, 1/fs) eta np.zeros_like(t) for k in range(M): eta amp[k] * np.cos(2*np.pi*f[k]*t phase[k]) return t, eta代码逻辑分几步先按等分频率生成代表频率数组再算出谱密度做一次能量修正然后乘以sqrt(2df)得到振幅最后逐个叠加余弦波。注意这一版频率点是固定的只有相位随机这是随机相位叠加法如果想要更好的随机性可以在频率区间内再做一次均匀随机抖动也就是频率随机叠加法。两者都能用但频率随机法生成的时历更不容易出现周期包络代价是实测频谱会略有偏差。这里有几个参数直接决定模拟质量。M建议取400到800太低会让波面出现明显的周期性太高则计算量增大。duration和fs决定了总点数模拟时长至少要覆盖一百个波浪周期比如Tp8s时至少要800秒否则极值统计不可靠。seed用于复现结果设计报告里换一个种子复核是常见的习惯。2.3 参数标定JONSWAP谱的\gamma取值与方向扩散gamma不是越陡越好它描述的是风浪成长中的谱峰集中度。对于有限风区gamma大约在3.3附近对于风区很短的年轻风浪gamma可能更高到5左右对于涌浪或衰减后的波浪gamma接近1。我们一般在代码里先固定gamma3.3但实际项目里如果有实测频谱最好用最小二乘法拟合得到gamma标定值。另一个容易被忽略的参数是方向扩散。上面的代码只模拟了单向长峰波。真实海面是短峰波能量分布在多个方向上。如果结构是细长柱体方向的影响还不太大如果是大尺寸结构或浮体方向谱的作用就不可忽视。最简单的方向分布函数是cos^2形式D(theta) (2/pi)*cos^2(theta)把谱密度乘以方向分布后再把每个方向分量叠加。这样得到的波面时历会明显更“乱”统计也更接近真实海面。对于单柱Morison力计算我一般先做单向波确认整个流程正确后再升级到方向谱。因为单向波算出的力通常比多向波略偏大相当于保守估计这在对标物模时需要注意。下面给一组我用过的标定参数供参考参数典型值说明Hs3.0 m有效波高设计波高统计值Tp8.0 s谱峰周期常取谱峰对应的谱密度最高点周期gamma3.3谱峰增强因子涌浪取1fmin0.005 Hz低于此频率的能量通常可忽略但太低会增大计算量fmax0.6 Hz高频截断由实际风浪谱衰减决定M500频率分量数至少400duration1200 s足够覆盖150个波浪周期fs4 Hz采样频率必须高于关注最高频率的2倍这组参数在计算中能保持频谱形状稳定力极值不会因为模拟时长太短而飘得太厉害。如果只是为了演示可以把M和duration降下来但这样得到的结果只能看趋势不能直接作为设计值。3. 水质点速度与加速度线性理论的计算细节3.1 线性波理论的速度公式与水深修正波面模拟完成后下一步是求水质点速度。线性波理论Airy波在中等水深和深水条件下表现稳定也是大多数Morison力计算的基础。对固定水深h、静水面z0、z轴向上水平速度u和垂直速度w的单个分量为u a * ω * cosh(k(zh)) / sinh(kh) * cos(kx - ωt φ) w a * ω * sinh(k(zh)) / sinh(kh) * sin(kx - ωt φ)其中k是波数满足色散关系ω^2 g k tanh(kh)a是每个分量的振幅ω2πf。这个公式里最关键的是cosh和sinh的比值。在深水区kh很大cosh(k(zh))/sinh(kh)约等于exp(kz)速度随水深指数衰减在浅水区这个比值变成一个线性函数速度沿水深更均匀这也是浅水区和深水区速度剖面差异的根源。写代码时不要直接用深水波数也不要忽略色散关系直接取kω^2/g。很多翻车案例都是因为省了这一步结果深水区看似差不多浅水区速度偏大20%以上。牛顿迭代求解色散关系是最快的方式初值取深水波数ω^2/g即可。下面是一个独立的波数求解函数def wave_number(omega, depth, g9.81): 用定点迭代求解色散关系 omega^2 g*k*tanh(k*h) :param omega: 角频率 (rad/s) :param depth: 水深 (m) :return: 波数 k (1/m) k omega**2 / g for _ in range(50): k_new omega**2 / (g * np.tanh(k * depth)) if np.abs(k_new - k) 1e-10: break k k_new return k这个迭代法的收敛速度很快对于水深特别浅的极端情况初值选深水解可能导致迭代震荡这时可以限制每次迭代的步长或者改用二分法。不过常规海洋工程水深下定点迭代已经够用。注意水深h和z的符号在代码里必须统一我习惯把海底定义为-h这样zh在海底为0、水面为h避免符号混乱。3.2 从频域到时间域速度时历的生成上一章生成的波面时历是许多频率分量的叠加速度同样可以按频率叠加只不过每个分量要乘一个传递函数。水平速度对波面振幅的传递函数为H_u(f, z) ω * cosh(k(zh)) / sinh(kh)注意相位关系在波浪行进方向水平速度与波面同相垂向速度与波面差90度。这意味着垂向速度在波峰处为0在波面过零处最大。很多刚开始写代码的人把垂直速度相位搞错会造成垂向力相位完全不对。下面是生成速度与加速度时历的函数def linear_kinematics(t, f, phase, amp, h, z, g9.81): 基于线性理论计算给定水深点的水平、垂向速度与加速度 :param t: 时间数组 :param f: 频率分量数组 (Hz) :param phase: 各分量相位 (rad) :param amp: 各分量振幅 (m) :param h: 水深 (m) :param z: 目标水深坐标 (m, 静水面为0, 向上为正) :return: u, w, du_dt, dw_dt omega 2 * np.pi * f u np.zeros_like(t) w np.zeros_like(t) du_dt np.zeros_like(t) dw_dt np.zeros_like(t) for k in range(len(f)): K wave_number(omega[k], h, g) theta -omega[k] * t phase[k] cosh_ratio np.cosh(K * (z h)) / np.sinh(K * h) sinh_ratio np.sinh(K * (z h)) / np.sinh(K * h) u amp[k] * omega[k] * cosh_ratio * np.cos(theta) w amp[k] * omega[k] * sinh_ratio * np.sin(theta) du_dt amp[k] * omega[k]**2 * cosh_ratio * np.sin(theta) dw_dt amp[k] * omega[k]**2 * sinh_ratio * np.cos(theta) return u, w, du_dt, dw_dt代码里theta取了x0表示作用点位置在原点实际结构坐标可以通过加kx来平移。水平加速度比水平速度相位领先90度所以cos变成sin垂向加速度比垂向速度领先90度所以sin变成cos。加速度在频域相当于乘iω但在时域实现时直接用解析导数更稳。这个函数的时间复杂度是M乘N如果M500、N4000就是两百万次循环Python跑完大约几秒可以接受。但如果要计算整个桩身沿深度几十个点就要用数组化重写避免三层循环。后面第4.3节我会给出沿桩身积分的思路。3.3 流速剖面插值不同水深处的值波浪力沿桩身是变化的不能只取水面一点的速度。通常做法是在z方向取若干代表点比如从海底到海面均匀取20个点每点算一条速度时历然后按Morison方程计算每点单位长度力再沿深度做数值积分。这样做会得到一个总的力时历而不是简单用表面速度乘以面积。近海底处垂直速度因为sinh(0)0而趋于0但水平速度并不为0。势流理论下海底水平速度仍存在实际流体中边界层内速度会衰减但对Morison力影响不大。所以不要直接把海底点的速度设成0除非你的模型里额外加了边界层修正。插值时要注意如果柱体穿过多个水深范围海洋结构高程坐标要统一。一般把海底定义为z-h静水面为z0向上为正。z0以上是自由水面不能随便取点因为线性理论在自由水面上仍然适用但z0处已超出静水面波浪表面实际可能淹没或暴露这时候要结合瞬时波面做湿表面判断。更粗略的做法是只计算到静水面然后用波峰和波谷两条路径包络工程上也能接受。4. 波浪力计算Morison方程与适用范围4.1 Morison方程拆解惯性力拖曳力直径相对波长较小的细长圆柱体波浪场基本不被柱体扰动此时波浪力可以用Morison方程计算。总水平力由两部分组成F F_drag F_inertia其中拖曳力正比于速度绝对值乘速度方向固定惯性力正比于加速度。对于单位长度柱体f(z,t) 0.5 * ρ * C_D * D * |u(z,t)| * u(z,t) ρ * C_M * π * D^2 / 4 * du/dt(z,t)第一项是速度平方项第二项包含了惯性质量和附加质量效应。公式里的D是柱体直径ρ是海水密度C_M等于1加上附加质量系数C_a。注意这个公式在速度方向变化时拖曳力要用|u|u来处理否则力的方向会乱。适用范围通常用D/L和KC数判断。D/L不应太大一般要求D/L0.2否则需要考虑绕射效应改用大尺度结构物上的Froude-Krylov公式。KC数是衡量拖曳力和惯性力相对重要性的无量纲数KC较小时惯性力主导KC很大时拖曳力主导。对于常见的导管架桩腿D通常在1到2米波长为50到100米时D/L在0.02到0.04Morison方程完全够用。下面是一个单位长度力计算函数可以直接接入前面的速度时历def morison_per_length(u, du_dt, D, rho1025.0, Cd0.8, Cm1.8): 计算单位长度柱体的Morison波浪力时历 :param u: 水平速度时历 (m/s) :param du_dt: 水平加速度时历 (m/s^2) :param D: 柱体直径 (m) :param rho: 海水密度 (kg/m^3) :param Cd: 拖曳力系数 :param Cm: 惯性力系数 :return: 单位长度力时历 (N/m) force_drag 0.5 * rho * Cd * D * np.abs(u) * u force_inertia rho * Cm * (np.pi * D**2 / 4.0) * du_dt return force_drag force_inertia4.2 力系数C_D、C_M的选取与敏感性C_D和C_M是Morison方程里最敏感、也最“有黑匣子味”的两个参数。规范上一般给出范围C_D在0.6到1.2之间C_M在1.5到2.0之间。实际取值要看结构形式、表面粗糙度和振动响应。对光滑圆管C_D通常取0.7左右对数生锈或有海生物附着的表面C_D会明显增大甚至达到1.5。C_M则更依赖KC数。在KC数较小时流动分离不明显C_M接近2KC数增大后尾流效应让惯性力系数下降大约在1.3到1.6之间。工程上如果没有明确数据我先取C_D0.8、C_M1.8后面再用不同取值做敏感性包络。这两个系数对力极值的敏感性不同在速度峰值处拖曳力占主导在波面过零处加速度达到峰值惯性力占主导。如果计算出的合力时历在过零点附近出现大尖峰多半是惯性力系数C_M过大或者加速度相位错误。反过来的情况也遇到过拖曳力系数设得太大导致峰值出现在速度最大时整个时历波形看起来特别胖。我一般会做一个系数敏感性表算三组值工况C_DC_M极值变化基准0.81.8—光滑管0.61.6下降约10%左右粗糙管1.22.0上升约15%左右这个表不用做太多能看出力的极值跨度和包络就够了。报告里写清楚系数依据比随便给一个数更禁得起推敲。4.3 波浪力时历与极值统计把单位长度力沿水深积分就得到整根柱体上的水平力。常用梯形积分或辛普森积分。先在不同水深点计算速度然后调用morison_per_length得到每个节点的力最后积分from scipy.integrate import simpson def total_force_on_pile(z_nodes, u_per_node, du_per_node, D, Cd0.8, Cm1.8): 沿水深积分计算整根桩的波浪力时历 :param z_nodes: 水深节点坐标数组 (m) :param u_per_node: 每个节点的水平速度时历, shape (节点数, 时间点数) :param du_per_node: 每个节点的水平加速度时历 :return: 总力时历 (N) f_per_node np.zeros_like(u_per_node) for i in range(len(z_nodes)): f_per_node[i] morison_per_length(u_per_node[i], du_per_node[i], D, CdCd, CmCm) force_total simpson(f_per_node, xz_nodes, axis0) return force_total, f_per_node如果柱体是倾斜的Morison方程还要考虑沿轴向的切向力如果是锥形管或变直径管D要按水深分段取。这段代码先按等直径直立柱处理变直径时只要把D改成一个按深度取值的数组即可。积分之后得到总力时历然后统计最大极值和最小极值。注意极值分布会受模拟时长影响如果只算一条时历最大力可能是这条时历里碰巧出现的一个大波造成的。所以在初步设计里我会用这个极值去做结构截面校核但最终设计一定要做多种子极值分析具体见第5章。5. 避坑指南随机波浪模拟与波浪力计算的五个常见问题5.1 模拟序列不平稳频率截断与能量泄漏现象生成的波面时历前500秒和后500秒的标准差差了20%以上甚至出现明显的周期包络。原因频率分量数M太少或频率下限截断过高让低频率分量的周期长度接近甚至超过模拟时长造成波面在时历两端没有充分采样。另一个可能性是相位种子和固定频率区间的中心重合度太高导致各分量相位差有规律波面产生了伪周期。解决把M提高到500以上频率下限降到0.005Hz甚至更低并检查波面时历的自相关函数。自相关函数应该随滞后时间迅速衰减如果衰减后又回升说明序列里有周期性成分。这时换一个随机种子再试若序列形态变化太大多半还是分量数不够。5.2 水质点速度与加速度的高频抖动现象波面看起来还算光滑但速度时历尤其是靠近水面点的速度像有高频噪声一样抖动剧烈。原因高频分量的波数k很大cosh(k(zh))在z接近0时增长很快导致高频速度传递函数被放大。虽然波面谱本身在高频段已经很小但经过传递函数放大后速度谱的高频尾部仍然有能量表现为速度时历上的小尺度波动。解决在计算速度前对谱做一次高频截断或滤波。通常把模拟使用的fmax作为截断上限实际计算速度时再用一次略低的上限比如f_velocity_max min(fmax, 2*fp)。还可以对速度时历做低通滤波但更好的办法是从频域相位上直接避免高频分量过大。5.3 Morison方程力的尖峰与系数敏感现象总力时历在波面过零附近出现非常窄的双向尖峰看起来和规则波浪力完全不一样。原因过零时加速度最大惯性力也最大。如果加速度项在时域用差分计算高频分量会被差分运算放大产生数值尖刺。规范做法是在频域直接乘iω得到加速度或者用前面的解析导数方法不要用中心差分。解决检查加速度生成方式确保使用传递函数乘以角频率的平方。同时检查C_M取值在KC较大时C_M应适当降低。把C_M从2.0降到1.6尖峰峰值往往能下降15%左右这个变化正常。5.4 浅水区线性理论失效非线性波的取舍现象在浅水或近岸场地按线性理论算出的波浪力比实测波力大很多速度剖面和波面形状也不对波峰处没有尖瘦形状反而很平。原因线性波理论假设波陡很小浅水里波高和水深比值较大时非线性效应变得十分显著自由表面条件需要保留二阶甚至更高阶项。此时波峰处水质点速度并非线性理论给出的余弦分布。解决用相对水深和波陡判断当kh大于π/2且波陡H/L小于1/50线性理论够用当kh小于π/2或波陡大于1/25改用Stokes五阶波或流函数波。随机波浪模拟中可以先把每个频率分量按线性叠加生成波面再对每个大波做瞬时自由表面修正工程上复杂一些。若只是初步筛选可以在线性结果基础上乘一个1.1到1.2的经验放大系数但最终设计不能这么做。5.5 极值统计偏差随机种子与模拟时长现象同一个Hs和Tp只是换了一下随机种子最大波浪力从580kN变成720kN差了20%以上。原因随机波浪时历只是海况的一种实现极值本身就是随机变量。短时历模拟往往只能覆盖有限个波浪周期最大波高没有机会充分显现。模拟时长不够长时极值完全取决于该时段内偶尔出现的特大波。解决把模拟时长提高到至少3小时10800s并生成20到50个独立种子对每个种子提取极值然后用Gumbel分布拟合用外插得到100年重现期极值。这个流程虽然计算量大但比仅仅拉长单次模拟要稳健得多。我一般会把所有种子的极值画成概率图如果线性拟合度好说明极值分布已收敛。6. 验证与进阶用能量谱闭合校验计算结果整个流程跑完后别急着出结果先做一次能量谱闭合校验。方法很简单对生成的波面时历做傅里叶变换用Welch方法估计实测谱然后和输入的目标JONSWAP谱画在同一张图上。两者的面积即0阶矩应该一致Hs误差在1%以内谱峰频率和宽度也应该大致重合。下面是一段校验代码from scipy.signal import welch def check_spectrum(eta, fs, f_input, S_input): 比较生成波面的实测谱与目标谱 freqs, S_est welch(eta, fsfs, npersegmin(len(eta)//4, 4096), noverlaplen(eta)//16) df_est freqs[1] - freqs[0] Hs_est 4 * np.sqrt(np.sum(S_est) * df_est) df_target f_input[1] - f_input[0] Hs_input 4 * np.sqrt(np.sum(S_input) * df_target) print(f目标Hs: {Hs_input:.2f} m, 实测Hs: {Hs_est:.2f} m) return freqs, S_est速度谱的校验也类似。理论上速度谱等于波面谱乘以传递函数的平方S_u(f) S_eta(f) * (ω * cosh(k(zh))/sinh(kh))^2对比实测速度时历的谱和这个理论谱可以快速发现传递函数代码里的乘除错误。我见过一个脚本把cosh和sinh写反了波面谱校验通过了但速度谱完全不对结果力小了一半。从那以后我每次做随机波浪模拟最后都会强制走一遍波面和速度的双重谱闭合校验否则力和弯矩的极值我都不敢对外报。希望这份笔记里的实现链路和参数坑能帮你把随机波浪速度及波浪力计算这条路走顺。本文还有配套的精品资源点击获取