简介面向海洋工程随机波浪分析场景的配套计算资料包适合需要将Jonswap谱、小振幅理论与Morison方程落到代码中的学习者重点解决随机波浪序列生成、波浪速度场求解及小尺寸结构物波浪力估算问题。压缩包共2个文件包含1个MATLAB脚本和1个PNG图示整体约39KB脚本对应核心计算流程图示用于辅助理解与结果核验。已有880人学习说明该主题在课程设计、科研入门和工程复核中有一定参考需求。读者可据此梳理从随机波浪谱到水质点水平、垂向速度分布再到Morison方程计算波浪力的完整链路并借助图示检查参数设置、结果趋势与实现思路适合作为海洋结构物耐波性分析的入门示例和代码参考。1. 随机波浪速度及波浪力计算从谱到时程别把随机当成随便随机波浪速度及波浪力计算最容易翻车的地方不是公式而是把谱参数、相位和积分步长当成三件互不相干的事。我见过不少计算书把有效波高、谱峰周期填得很漂亮速度时程却拿规则波乘一个随机相位最后桩柱力峰值和疲劳幅值都对不上。这个方向解决的是给定水深、波高、周期和桩径用谱离散重构波面再由线性波理论得到水质点速度、加速度沿水深积分拖曳-惯性力输出随机波浪力时程和统计量。新手可以按后续代码跑通最小闭环熟手更该盯频率截断、分段积分和系数取值边界。适合海上结构初设、疲劳谱分析、课程设计复现。2. 谱离散与时程重构把随机波面变成可积分的速度场随机波浪速度及波浪力计算的第一步不是写力公式而是把一条目标谱变成时间域里可重复的波面、速度和加速度。谱只告诉你能量随频率怎么分布真正送到力公式里的是每个频率分量的幅值、相位、波数和沿水深的垂向衰减。这里最怕两种做法一是只给一个有效波高和谱峰周期直接生成一条“看起来像波浪”的正弦叠加二是频率范围拍脑袋截断低频少一段、高频多一段反算有效波高就对不上。先立住谱离散这条线后面的力积分才有意义。2.1 水深、有效波高和谱峰周期怎么落到频率网格水深 h、有效波高 Hs、谱峰周期 Tp 是输入但频率网格不是随便取。总时长 T_total 决定最低角频率 ω_min2π/T_total频率分辨率 Δω 也由它决定最高角频率一般取 3 倍谱峰角频率或者看累计能量覆盖到 99% 以上再截。频率分量数 Nf 不需要和时程采样点数一样但 Nf 太少会让谱峰附近变得很粗糙速度峰值会被削平。我的习惯是先用 2048 个频率分量做最小闭环再把 Nf 和 T_total 各放大一倍看统计量是否稳定。有效波高和谱之间还有一个缩放动作。常用单峰谱的原始积分面积不一定刚好等于 (Hs/4)^2所以生成谱后要用 m0ΣS(ω)Δω 反算再把整条谱乘上 (Hs/4)^2/m0。这个动作不做后面反算 4 倍波面标准差时就会对不上读者会以为速度公式错了其实是谱能量没归一。谱峰增强因子 γ 常见取 3.3 左右γ 越大谱峰越尖时程里大浪群和安静段的对比越明显做疲劳分析时这个参数别随手改改完要重新看统计样本。2.2 用色散关系生成随机波浪速度及加速度时程色散关系 ω²gk tanh(kh) 是连接频率和波数的黑匣子深水、浅水都不能直接拿 ω²/g 一把梭。我的做法是用牛顿迭代初值按水深做分区浅水取 ω/√(gh)深水取 ω²/g迭代 30 次基本够用。得到每个频率分量的 k 后波面写成 Σ a_i cos(ω_i t - k_i x φ_i)速度在 z 高程处乘上 cosh[k_i(zh)]/sinh(k_i h)加速度再对时间求导相位从 cos 变 -sin。下面这段代码是最小闭环先把波面、速度、加速度三条时程跑出来。import numpy as np g 9.81 h 30.0 Hs 3.0 Tp 8.0 gamma 3.3 T_total 3600.0 dt 0.1 omega_p 2 * np.pi / Tp omega_min 2 * np.pi / T_total omega_max 3.0 * omega_p Nf 2048 omega np.linspace(omega_min, omega_max, Nf) domega omega[1] - omega[0] def single_peak_spectrum(w, wp, gamma): # 常用单峰谱形式sigma 按峰前峰后分段 sigma np.where(w wp, 0.07, 0.09) r np.exp(-((w - wp) ** 2) / (2 * sigma ** 2 * wp ** 2)) S 0.0081 * g ** 2 / w ** 5 * np.exp(-1.25 * (wp / w) ** 4) * gamma ** r return S S single_peak_spectrum(omega, omega_p, gamma) m0 np.sum(S) * domega S S * (Hs / 4.0) ** 2 / m0 # 缩放到目标有效波高 def solve_k(w, h, g9.81): # 浅水/深水初值分区再牛顿迭代 k np.where(w ** 2 * h / g 1.0, w / np.sqrt(g * h), w ** 2 / g) for _ in range(30): f g * k * np.tanh(k * h) - w ** 2 df g * np.tanh(k * h) g * k * h / np.cosh(k * h) ** 2 k k - f / df return k k solve_k(omega, h, g) rng np.random.default_rng(202406) phase rng.uniform(0, 2 * np.pi, Nf) amp np.sqrt(2 * S * domega) t np.arange(0, T_total, dt) z 0.0 # 水面附近水质点算桩柱力时要换成多个 z 节点 eta np.zeros_like(t) u np.zeros_like(t) du np.zeros_like(t) for i in range(Nf): theta omega[i] * t - k[i] * 0.0 phase[i] eta amp[i] * np.cos(theta) coef omega[i] * np.cosh(k[i] * (z h)) / np.sinh(k[i] * h) u amp[i] * coef * np.cos(theta) du -amp[i] * omega[i] * coef * np.sin(theta) print(Hs, 4 * np.std(eta), np.max(np.abs(u)), np.max(np.abs(du)))逻辑说明先按总时长和谱峰周期切频率网格再用单峰谱算能量密度缩放使积分面积对应目标有效波高。色散关系用牛顿迭代求波数相位用固定种子随机抽取保证同一组输入每次结果可复现。循环里同时累加波面、水平速度和水平加速度速度系数只跟频率、波数和水深有关跟时间无关所以可以提前算成数组。参数说明h 是水深单位 mHs 是有效波高单位 mTp 是谱峰周期单位 sgamma 控制谱峰尖锐程度T_total 决定最低频和统计样本长度dt 是时间步长必须满足最高频率分量的采样要求omega_max 取 3 倍谱峰角频率是常用截断Nf 是频率分量数z 是水质点高程z0 表示水面附近算桩柱力时要从水底到水面布节点。随机相位不是玄学种子固定后结果可复现换种子只影响单条时程的峰值位置不影响谱统计量。2.3 频率截断与时间窗能量覆盖率和周期分辨率频率截断最容易出现“谱看着对、时程不对”。ω_min 太低没有意义因为总时长只有 3600 s最低频分量一个周期就占满全窗ω_max 太低会把高频小波动全丢掉速度峰值偏小拖曳项因为带绝对值平方偏低会更明显。我一般看累计能量把 S(ω)Δω 从 ω_min 累加到 ω_max覆盖到 99% 以上再停如果只覆盖 95%先别急着算力把 ω_max 抬高再看反算 Hs。时间窗也要够长至少覆盖 100 个以上谱峰周期否则统计有效波高和力峰值都受单条时程的随机性影响。时间步长 dt 不是越小越好。dt 太小会生成巨大的数组后面沿水深积分时内存吃紧dt 太大又会让速度峰值被削平加速度对时间求导后误差更大。一个实用检查是看 4 倍波面标准差和输入 Hs 的误差控制在 5% 以内再往下走。还要看速度时程的最大值和均方根换 dt 减半后峰值变化小于 2%这个时间步长才算过关。注意先检查 4*np.std(eta) 与 Hs 的误差再检查换 dt 和换种子后的统计量别一上来就盯着单条力峰值。3. 拖曳-惯性力方程分段积分随机波浪力从公式到代码随机波浪速度及波浪力计算走到力这一步核心是小直径构件沿水深分段积分。单位长度上的力分成拖曳项和惯性项拖曳项与 |u|u 成正比惯性项与 du/dt 成正比。拖曳项是非线性的速度正负换向时力也换向绝对值乘速度不能漏惯性项是线性的但相位和速度差 90 度写代码时 sin 和 cos 一错力峰值就整体跑偏。工程上常见的做法是沿水深布 60 到 120 个节点用梯形积分把单位长度力累加成总力再用时间步长收敛和分段数收敛两道检查兜底。3.1 拖曳项与惯性项的量级判断什么时候不能线性化拖曳项和惯性项谁主导不用等算完才知道。把速度幅值写成 aω拖曳项量级约 0.5ρCD D (aω)²惯性项量级约 ρCM πD²/4·aω²两者比值大致正比于 (4CD a)/(πCM D)。波幅 a 相对桩径 D 越大拖曳越占优桩径越大、波幅越小惯性越占优。惯性主导时把拖曳项线性化还能接受误差主要落在小力段拖曳主导时|u|u 的尖峰和换向非常明显线性化会把力峰值和疲劳幅值同时带偏这时候老老实实时域积分。我一般先做一张量级表给定 Hs、Tp、D算几个代表频率的速度幅值估拖曳项和惯性项的最大值。若拖曳项不到惯性项的 20%可以先用线性化做快速估算若超过 50%不要省非线性积分。这个判断还能帮你检查 CD、CM 有没有填反CD 影响拖曳项CM 影响惯性项两者改一个力量级变化的方向不一样先看谁主导再调参数比盲目试数快。3.2 沿水深分段积分随机波浪力的最小代码下面的代码假定第 2 章已经算好了每个 z 节点上的 u_tz 和 du_tz形状是“时间 × 水深节点”。拖曳项先算单位长度力惯性项也算单位长度力再沿水深做梯形积分。节点从 -h 到 0水面节点和水底节点都要有水底速度按线性波理论趋近于零但不要直接删掉否则积分区间短一截。import numpy as np rho 1025.0 D 1.2 CD 1.1 CM 1.9 h 30.0 z_nodes np.linspace(-h, 0.0, 61) # 从水底到水面含两端 def force_time_series(u_tz, du_tz, z_nodes, D, CD, CM, rho): dz np.diff(z_nodes) # 单位长度拖曳力和惯性力形状均为 (时间, 水深节点) drag_density 0.5 * rho * CD * D * np.abs(u_tz) * u_tz inertia_density rho * CM * np.pi * D ** 2 / 4.0 * du_tz # 对每个时间步沿水深做梯形积分 F_drag np.sum(0.5 * (drag_density[:, 1:] drag_density[:, :-1]) * dz[None, :], axis1) F_inertia np.sum(0.5 * (inertia_density[:, 1:] inertia_density[:, :-1]) * dz[None, :], axis1) return F_drag F_inertia, F_drag, F_inertia # 实际使用时把下面两行替换成第 2 章按每个 z 节点算出的矩阵 # u_tz ... # du_tz ... # F, Fd, Fi force_time_series(u_tz, du_tz, z_nodes, D, CD, CM, rho) # print(F.max(), F.min(), F.std())逻辑说明拖曳项先做 np.abs(u)*u保证速度换向时力也换向惯性项直接乘加速度。沿水深方向用相邻节点平均值乘 dz等价于梯形积分。代码把拖曳力和惯性力分开返回方便看谁主导也方便做系数敏感性。实际跑的时候u_tz 和 du_tz 的列顺序必须和 z_nodes 一一对应水面节点不要漏。参数说明rho 是水密度工程上常用 1025 kg/m³ 左右D 是桩径单位 mCD 是拖曳系数常见取值在 1.0 到 1.2 附近CM 是惯性系数常见取值在 1.8 到 2.0 附近且 CM1附加质量系数z_nodes 从 -h 到 0dz 是相邻节点间距。分段数先取 60再翻倍看总力峰值变化变化小于 2% 再固定。注意CM 不是附加质量系数本身填表时别把 0.9 直接填进 CM否则惯性项会小一半。3.3 CD、CM 与桩径 D 的取值边界和敏感性CD、CM 不是越“保守”越好改大了拖曳项和惯性项都跟着涨但两者对峰值和疲劳的贡献不一样。拖曳主导的工况CD 改 10% 可能让力峰值改 10% 以上惯性主导的工况CM 改 10% 影响更直接。我的血泪经验是先把 CD、CM 各做 ±10% 的敏感性表再看桩径 D 变化时拖曳项和惯性项谁涨得快。D 在拖曳项里是一次方在惯性项里是二次方桩径变大时惯性项占比会明显上升原来拖曳主导的工况可能变成惯性主导这时候再沿用旧系数就会偏。取值边界还要看表面状态和流速范围。粗糙度、附着物、群桩遮蔽都会让 CD 往上走惯性系数在非定常流动里也会偏离常用区间。做初设可以用常见区间扫一遍做详细校核时至少给出一组下限、一组上限别只交一个中间值。随机波浪力本身有统计离散性系数再只给单点后面复核很难判断差异是系数造成的还是随机种子造成的。4. 方向谱与群桩叠加随机波浪力不能只算单根桩单桩算通以后很多人直接把单桩力乘桩数当群桩合力这是随机波浪力计算里很典型的偷懒。群桩之间既有波向投影带来的相位差也有上游桩对下游桩的遮蔽还有方向谱把能量分散到不同波向角。单方向谱只给一条主波向所有桩的相位差按同一波向算方向谱每个频率分量带一个波向角相位差要在 x cosβ y sinβ 上投影。两者在合力峰值和疲劳幅值上差别不小尤其是桩距接近一个波长的时候。4.1 方向离散与波向角单方向谱和方向谱的差别单方向谱适合主波向很明确、桩距不大的初算。做法是把每个频率分量的波向角固定成主波向相位项里只保留 kx cosβ。方向谱则要把能量按方向分布拆开常见做法是每个频率分量配若干方向角方向分布函数在主波向附近最大两侧衰减。方向离散数不用一上来取 72 个先取 12 到 24 个方向看群桩合力峰值和单方向谱差多少再决定要不要加密。方向谱最容易错在归一化。方向分布函数对方向积分应该等于 1离散后要检查 ΣD(β)Δβ 是否接近 1不接近就先把方向分布归一再乘频率谱。另一个坑是波向角的定义有人用“来波方向”有人用“传播方向”两者差 180 度群桩相位差会整体反号。我的习惯是代码里统一写成传播方向注释写清楚输出图里标一个箭头别让后面的人猜。4.2 群桩相位差下的波浪力叠加群桩相位差的核心是每个频率分量在不同桩位的相位不同。对第 j 根桩相位项写成 ωt - k(x_j cosβ y_j sinβ) φ。单桩力时程不能直接平移因为速度、加速度的幅值也随频率变化正确做法是把每个频率分量的贡献先算到每根桩再叠加。下面这段代码只演示相位差实际力要把每个频率分量的幅值、速度系数、加速度系数带进去。import numpy as np piles np.array([[0.0, 0.0], [6.0, 0.0], [0.0, 6.0], [6.0, 6.0]]) # x,y 单位 m beta np.deg2rad(30.0) # 传播方向相对 x 轴 omega np.array([0.6, 0.8, 1.0]) k np.array([0.08, 0.12, 0.18]) # 实际应由色散关系逐频率求解 phase np.array([0.2, 1.1, 2.4]) t np.linspace(0, 60, 600) def phase_sum(t, x, y, beta, omega, k, phase): theta np.zeros_like(t) for i in range(len(omega)): proj x * np.cos(beta[i] if np.ndim(beta) else beta) y * np.sin(beta[i] if np.ndim(beta) else beta) theta np.cos(omega[i] * t - k[i] * proj phase[i]) return theta for j, (x, y) in enumerate(piles): th phase_sum(t, x, y, beta, omega, k, phase) print(j, th.max(), th.min())逻辑说明每根桩的相位差由位置在波向上的投影决定x cosβ y sinβ 越大相位滞后越多。单方向谱时 beta 是标量方向谱时 beta 是每个频率分量对应的方向角数组投影要逐分量算。代码里只累加余弦相位实际力还要乘各分量的速度幅值和力系数。参数说明piles 是群桩坐标beta 是波向角omega、k、phase 是频率分量参数t 是时间向量。桩距接近半波长时相位差最容易让合力出现抵消或增强做群桩合力时至少扫 3 个波向角不要只算主波向。4.3 遮蔽效应不要拍脑袋经验系数和验证工况遮蔽效应是上游桩让下游桩流速降低常见做法是给下游桩乘一个 0.8 到 1.0 的遮蔽系数但系数不是固定的。桩距小、波向和桩列方向一致时遮蔽明显波向斜穿桩列时下游桩可能只被遮一部分。我的做法是先算无遮蔽的群桩合力再按桩距和波向给一组遮蔽系数上下限看合力峰值变化范围。如果遮蔽系数一改合力峰值变化超过 15%说明这个工况对遮蔽很敏感不能只给单点值。验证工况可以设三组单桩、双桩顺波向、双桩垂直波向。单桩用来检查相位和力公式双桩顺波向看遮蔽和相位差叠加双桩垂直波向看波向投影。三组都跑通再上四桩、九桩。群桩合力不是单桩力乘桩数也不是简单取最大值相加而是先把每根桩的时程算出来再在同一时间轴上叠加。叠加前检查时间轴是否一致别一根桩用 0.1 s 步长、另一根用 0.2 s 步长插值错了后面全乱。工况桩距波向主要检查单桩无主波向力峰值与解析规则波退化双桩顺波向0.5 倍波长沿桩列相位差与遮蔽系数双桩垂直波向0.5 倍波长垂直桩列波向投影是否写反5. 随机波浪力计算避坑五条从异常现象倒推原因的记录随机波浪力计算的坑大多不是公式不会而是中间某个环节悄悄偏了最后只在力峰值上表现出来。下面五条按“现象→原因→解决”写都是我在最小闭环里反复遇到的。每条都先用一个异常现象开头再倒推原因最后给可执行的检查动作。读者可以拿自己的结果对号入座别等交图前才发现。5.1 反算有效波高对不上谱缩放和频率截断的锅现象输入 Hs3.0 m生成波面后 4 倍标准差只有 2.6 m速度峰值也偏小。原因通常有两个谱没有缩放到目标有效波高或者 ω_max 截得太低高频能量没覆盖够。前者是 m0 没归一到 (Hs/4)^2后者是累计能量覆盖不到 99%。解决先算 m0ΣSΔω乘缩放因子后再生成时程再把 ω_max 从 2 倍谱峰频率逐步抬到 4 倍看 4 倍标准差是否稳定。若抬到 4 倍还差 10% 以上检查频率分辨率 Δω 是否太粗谱峰附近有没有被漏掉。还有一种隐蔽情况是 T_total 太短最低频分量一个周期就占满全窗谱的低频尾巴被强行截断。解决是把 T_total 加到至少 100 个谱峰周期再看统计量。别用单条时程的最大波高去反算 Hs随机波浪的有效波高是统计量单条时程的最大值本来就会波动。5.2 速度与加速度相位错位别把sin和cos写反现象力时程的峰值出现在波峰附近但理论上惯性项应该在波峰和波谷之间换向拖曳项才在波峰附近最大。原因多半是速度用 cos、加速度也用 cos或者加速度的负号漏了。线性波理论里速度与波面同相加速度超前 90 度写成 -sin 那一项不能少。解决拿一个单频规则波做退化画波面、速度、加速度三条曲线检查速度过零时加速度是否最大加速度过零时速度是否最大。单频都对不上随机叠加肯定错。另一个常见错位是沿水深系数用了 sinh 和 cosh 写反。水面附近速度最大水底趋近于零画一条沿水深的幅值剖面就能看出来。若水底速度反而比水面大别怀疑随机相位先查双曲函数。5.3 力峰值忽大忽小分段积分与时间步长没做收敛现象同一组谱参数只换随机种子力峰值差 20% 以上或者分段数从 30 加到 60峰值差 10%。原因可能是统计样本不够也可能是分段积分和时间步长没收敛。解决先固定种子把水深节点从 30 加到 60、120看总力峰值变化再把 dt 减半看峰值变化。两项都收敛到 2% 以内再换 3 个种子看统计分布。若换种子仍然差很多说明 T_total 太短单条时程的随机性太大不是公式问题。分段积分还要注意水底和水面节点。水底节点速度接近零但拖曳项在底部小惯性项也小漏掉底部一小段影响不大水面节点漏掉速度最大的一段直接丢掉力峰值会明显偏小。节点从 -h 到 0 布包含两端别用 -hdz 到 0-dz。5.4 CD/CM一改结果差一倍先看拖曳和惯性谁主导现象CD 从 1.0 改到 1.2力峰值涨了 30%CM 从 1.8 改到 2.0峰值只涨 5%。原因不是代码错而是这个工况拖曳主导。反过来大直径桩惯性主导时CM 改一点峰值就明显变。解决先输出拖曳力和惯性力的分量最大值算两者比值再决定调哪个系数。拖曳主导的工况CD 的敏感性表要做得细惯性主导的工况CM 和桩径 D 的敏感性要优先看。别把两个系数一起改否则分不清是谁的贡献。还要检查 CM 有没有把附加质量系数直接填进去。CM1附加质量系数附加质量系数 0.9 对应 CM1.9。填 0.9 会让惯性项小一半力峰值直接翻车。5.5 群桩合力比单桩相加还小方向谱、相位差和遮蔽系数混用现象四桩合力峰值比四根单桩峰值之和还小很多或者群桩合力曲线出现不合理的深谷。原因可能是相位差算错、方向谱没归一或者遮蔽系数乘在了合力上而不是单桩上。解决先关掉遮蔽只算相位差看群桩合力是否在单桩力之和的合理范围内再加方向谱检查方向分布积分是否为 1最后加遮蔽遮蔽系数只乘下游桩的单桩力不要乘总合力。波向角定义也要统一来波方向和传播方向差 180 度相位差会整体反号。群桩叠加还要检查时间轴。每根桩的力时程必须在同一时间向量上叠加不能各自用自己的 t。若桩位坐标单位是米、波数单位是 1/米投影项单位一致若坐标用了厘米相位差会差两个数量级合力曲线会变成一团乱麻。6. 验证与交付随机波浪力结果怎么让别人敢用6.1 四条验证线和一个交付习惯随机波浪力结果交出去之前我一般跑四条验证线。第一条是谱矩反算用生成的波面算 4 倍标准差和输入 Hs 比误差控制在 5% 以内第二条是单频退化把谱换成一条单频分量波面、速度、加速度、力时程都能和手算规则波对上第三条是收敛检查水深节点翻倍、时间步长减半力峰值变化小于 2%第四条是随机种子换 3 个种子看力峰值和均方根值的离散范围别只交一条时程。验证线做法合格线谱矩反算4×std(eta) 对输入 Hs误差小于 5%单频退化单分量规则波手算对比相位和幅值一致分段与步长收敛节点 30/60/120dt 减半峰值变化小于 2%随机种子3 个种子统计峰值离散范围可解释交付时我习惯附一张参数表h、Hs、Tp、γ、T_total、dt、Nf、z 节点数、CD、CM、D、波向角、遮蔽系数上下限。表里不写“常用值”写实际用的数和取值理由。别人复核时先看这张表再看反算 Hs 和收敛检查基本能判断结果敢不敢用。这个习惯让我少返工很多次力峰值对不上时先翻参数表往往能发现是某个系数或节点数在捣鬼。我自己的教训是先把单频规则波退化跑通再上随机谱每次换 CD/CM 先做量级表不直接看峰值。希望帮到你。本文还有配套的精品资源点击获取