MATLAB船舶波浪仿真全流程：从PM谱随机波生成到波浪力与垂荡运动Simulink建模

简介面向海洋工程、船舶水动力学与计算力学方向的学者和工程师这份MATLAB源码包围绕“船舶在波浪中的运动”这一主题提供了从波浪建模到受力分析、再到运动仿真的完整代码方案。资源仅7个文件以6个m脚本为主、1个txt说明文件为辅压缩包大小2KB轻量便捷。m脚本分别实现波浪模型初始化、线性波浪模拟、JONSWAP/PM谱随机波浪生成、波浪力计算及主仿真流程txt文件则可用于参数记录或使用说明。目前已有1555人学习下载实用性获得一定关注。通过阅读和运行这些代码读者可以快速掌握利用傅里叶级数构造线性波、基于海洋工程常用谱生成随机波以及计算波浪对船体作用力的核心方法并以此为起点扩展至Simulink船舶运动仿真或六自由度运动分析兼具教学与二次开发价值。1. 为什么把船舶波浪仿真跑在 MATLAB 里第一次把 PM 谱转成波面时我盯着结果看了十分钟——设定 4 米有义波高合成出来的波面最大值只有 2.6 米。不是代码写错是频率离散和相位种子出了问题。这类问题在船舶波浪仿真里非常典型理论公式谁都会背但真正要把线性波叠加、随机波谱、波浪力积分和运动方程串成一条能跑的链路中间隔着一堆工程细节。这套资源里main.m、bomian.m、linearWaveSimulation.m、waveModelInit.m、specturmPM.m、waveForce.m六个文件覆盖的正是这条链路从波浪场生成到船舶受力再到运动响应仿真。适合三种人一是海洋工程专业做课程设计或毕业设计的学生二是刚转入射流体力学的算法工程师三是做船舶操纵性仿真但一直用商用软件、想搞清楚底层计算逻辑的从业者。MATLAB 的优势在于矩阵运算、信号处理和 Simulink 无缝衔接改参数观察响应不需要重新编译这也是这类仿真首选它的原因。2. 波浪场建模从线性波叠加到 PM 谱随机波生成2.1 线性波理论为什么是起点线性波理论Airy 波理论的核心假设是波高远小于波长自由表面边界条件可以做线性化处理。在这个前提下任意不规则波面可以看作无数规则余弦波的线性叠加[ \eta(x,t) \sum_{i1}^{N} a_i \cos(k_i x - \omega_i t \varepsilon_i) ]其中 (a_i) 是第 (i) 个分量的振幅由目标谱密度决定(k_i) 是对应频率 (\omega_i) 的波数需要满足色散关系 (\omega^2 gk\tanh(kh))(\varepsilon_i) 是随机相位通常服从 ([0, 2\pi)) 均匀分布。这套叠加法也叫长峰波模型它假设所有波浪沿同一方向传播适合初步耐波性分析但不适合研究斜浪和短峰海况。linearWaveSimulation.m做的事情就是把公式翻译成 MATLAB 代码。写的时候有几个容易坑的地方一是频率范围必须覆盖谱的有效能量区间通常取 (0.3\omega_p) 到 (3\omega_p)(\omega_p) 是谱峰频率二是频率分量个数不能太少少于 50 个合成波面会出现明显的周期性三是随机相位种子必须固定否则每次运行结果都不一样后续验证没法复现。2.2 waveModelInit 参数初始化的常见做法waveModelInit.m在我的使用习惯里是入口函数负责把海况参数、船体几何参数、仿真时间步长统一封装成一个结构体。这个设计比散落一堆全局变量干净得多后面所有函数都从param结构体取数改参数不用动业务代码。function param waveModelInit(spectrumType, Hs, Tp, varargin) % 输入: spectrumType 谱类型 PM 或 JONSWAP % Hs 有义波高 (m), Tp 谱峰周期 (s) % 输出: param 结构体, 包含波浪、船体、仿真参数 p inputParser; addRequired(p, spectrumType, ischar); addRequired(p, Hs, (x) x 0); addRequired(p, Tp, (x) x 0); addParameter(p, h, 50); % 水深 (m)默认深水 addParameter(p, gamma, 3.3); % JONSWAP 谱峰因子 addParameter(p, Nfreq, 200); % 频率离散个数 addParameter(p, seed, 42); % 随机种子保证可复现 parse(p, spectrumType, Hs, Tp, varargin{:}); param.spectrumType p.Results.spectrumType; param.Hs p.Results.Hs; param.Tp p.Results.Tp; param.h p.Results.h; param.gamma p.Results.gamma; param.Nfreq p.Results.Nfreq; rng(p.Results.seed); % 频率范围0.3*omega_p 到 3*omega_p wp 2 * pi / Tp; wmin 0.3 * wp; wmax 3.0 * wp; param.w linspace(wmin, wmax, Nfreq); param.dw param.w(2) - param.w(1); param.phase rand(1, Nfreq) * 2 * pi; end这段代码里最值得说明的是inputParser的使用。很多新手会用 nargin 手动判断缺省值参数一多就乱套inputParser把必填项、可选项、类型校验集中处理后续扩展船体参数时只需要往addParameter里追加即可。频率范围取 (0.3\omega_p) 到 (3\omega_p) 不是随便定的PM 谱和 JONSWAP 谱 99% 的能量都落在这个区间低于 0.3 倍谱峰频率的分量对波面贡献极小高于 3 倍的成分能量也很小但会显著拖慢仿真。需要特别注意的是种子固定rng(42)不写的话第二次运行同一个脚本得到的是完全不同的波面序列这对标定waveForce里的系数很不友好。2.3 linearWaveSimulation 与 spectrumPM 的实现要点2.3.1 PM 谱离散化与振幅计算specturmPM.m文件名里的拼写不影响调用实现的是 Pierson-Moskowitz 谱它是国际船舶与海洋工程界描述充分发展风浪的标准谱模型。单参数形式适合给定风速的场景双参数形式(H_s) 和 (T_p)在耐波性分析中更常用[ S(\omega) \frac{5}{16} \cdot \frac{H_s^2}{\omega_p} \cdot \left(\frac{\omega_p}{\omega}\right)^5 \cdot e^{-\frac{5}{4}(\omega_p/\omega)^4} ]我在用的时候建议用双参数版本因为实际设计工况一般直接给有义波高和谱峰周期。振幅和谱密度的关系是[ a_i \sqrt{2 S(\omega_i) \Delta\omega} ]这里的 (2) 是单边谱到双边谱的转换因子很多人漏掉它导致合成的波面有效波高只有设定值的 (1/\sqrt{2}) 左右。我开头说的 4 米变成 2.6 米一半原因就在这。function S spectrumPM(w, Hs, Tp) % 双参数 PM 谱 % w: 角频率数组 (rad/s), Hs: 有义波高 (m), Tp: 谱峰周期 (s) wp 2 * pi / Tp; w w(:); % 转列向量保证外部广播正常 S (5/16) * (Hs^2 / wp) .* (wp ./ w).^5 .* exp(-(5/4) * (wp ./ w).^4); S(w 0) 0; % 频率为负无物理意义 endPM 谱的优点是形式简单、参数少缺点是无法刻画有限风区和不完全发展的海况。如果做港湾或近岸工程建议换成 JONSWAP 谱它多了峰增强因子 (\gamma)在谱峰附近把能量集中起来也比 PM 谱更容易出现高频尾部的数值问题。waveModelInit里预留了gamma参数切换到 JONSWAP 时对应代码如下function S spectrumJONSWAP(w, Hs, Tp, gamma) wp 2 * pi / Tp; sigma ones(size(w)); sigma(w wp) 0.09; sigma(w wp) 0.07; A exp(-(w - wp).^2 ./ (2 * sigma.^2 .* wp.^2)); S spectrumPM(w, Hs, Tp) .* gamma .^ A; end2.3.2 波面合成与 bomian.m 绘图bomian.m的功能是波面可视化但它不只是画一条曲线还能画出在不同时刻沿船长方向分布的波面高度这是后面算波浪力的输入。我在实际项目中习惯把波面输出成eta(t, x)矩阵横轴时间纵轴船体纵向位置这样waveForce可以直接对矩阵做切片计算。参数推荐值影响调参方向Nfreq150300波面平滑度与计算量折中低于 80 出现周期性频率上限 (3\omega_p)固定覆盖 99% 谱能量缩短仿真步长时需同步提高随机种子固定整数结果可复现调试验证时务必固定水深 (h)实际水深色散关系求解精度浅水时换用显式近似公式3. 波浪力计算waveForce 的切片积分与傅里叶处理3.1 波浪力为什么不能直接乘个系数不少初学者拿到波面eta后直接乘一个经验系数当波浪力这在概念上是错的。波浪对船舶的作用力分两类一是 Froude-Krylov 力由未受扰动的波浪压力场产生压力分布沿船体型线积分二是绕射力是船体存在导致波浪场改变产生的附加力。线性理论的习惯做法是先算 F-K 力再用修正系数或者通过边界元求解计及绕射效应。waveForce.m在实际工程里常用切片法strip theory近似沿船长把船体切成若干二维剖面每段计算单位长度上的波浪力再沿船长积分。这样做计算成本低适合初步设计阶段快速扫参缺点是忽略三维流动效应对船长与波长比小于 1 的短船误差较大。切片法的基础公式是[ dF(x,t) \left[ \rho g A(x) \eta(x,t) - C_a \rho A(x) \ddot{\eta}(x,t) \right] dx ]第一项是静水恢复力与 F-K 压力的合并效果第二项是附加质量力。实际上eta对时间的二阶导数要从波面求解比较稳妥的做法是在频域算完响应后转时域而不是对合成波面直接数值微分数值微分会把高频噪声无限放大。3.2 切片积分实现与参数说明function F waveForce(eta, param) % 切片法计算船舶纵向波浪力 % eta: 1xN 波面高度分布 (m) % param: waveModelInit 生成的结构体 % 包含 Bm (各切片船宽), dx (切片长度), Cd (阻力系数) Nsec length(param.Bm); F 0; for i 1:Nsec % 该切片瞬时淹没截面积简化线性化处理 A_sub param.Bm(i) * (param.draft eta(i)); A_sub max(A_sub, 0); % 防止出水面积负值 f_pressure param.rho * param.g * A_sub; f_viscous 0.5 * param.rho * param.Cd(i) * param.Bm(i) * ... abs(param.vertVel(i)) * param.vertVel(i); dF (f_pressure f_viscous) * param.dx; F F dF; end end注意代码里的param.rho、param.g、param.vertVel必须在waveModelInit里初始化vertVel是波面垂向速度理想情况下应该由线性波理论解析给出(v_z a \omega \cos(kx - \omega t))而不是对波面差分数值求导。我在实际跑数时发现Cd黏性阻力系数对结果影响最大它和海况、船体形状都相关建议范围是 0.61.2没有实验数据时取 0.8 起步然后观察船舶垂荡位移的幅值是否在合理区间内。另外param.draft是设计吃水不能用满载吃水代替否则算出来的波浪力偏大。3.3 时域力序列的傅里叶分析拿到F(t)时间序列后先不要急着往里喂运动方程先看看它的频谱是否合理。用fft做频谱分析时有个常见误区是直接对原始信号做变换而不处理均值趋势导致零频成分淹没低频能量。我一般会先把均值去掉再变换Fs 1 / param.dt; % 采样频率 (Hz) L length(F); Y fft(F - mean(F)); % 去均值消除零频泄漏 P2 abs(Y / L); % 双侧幅值谱 P1 P2(1:L/21); % 单侧取前半段 P1(2:end-1) 2 * P1(2:end-1); % 能量加倍 f_axis (0:L/2) * Fs / L; % 频率轴 (Hz) plot(f_axis, P1) set(gca, XScale, log)这段代码里有两个细节值得展开。第一abs(Y/L)得到的是幅值谱不是功率谱密度两者量纲不同如果要和 PM 谱对比需要计算功率谱密度 (S_f |Y|^2 / (Fs \cdot L))幅值谱的峰值对应波频分量的振幅而功率谱密度下的面积等于信号方差。第二单侧谱的幅值要乘 2除了直流分量和奈奎斯特频率点这是做单边化时最容易漏的一步。看到谱峰频率落在 0.61.2 rad/s 区间对应谱峰周期 510 秒并且和输入的Tp一致就说明波浪力计算链路基本正确。4. 船舶运动响应仿真Simulink 与状态空间对接4.1 垂荡运动的二阶微分方程建模波浪力算出后船舶运动仿真的核心是求解刚体运动方程。以垂荡heave为例最常用的线性化模型是[ (m A_{33}) \ddot{z} B_{33} \dot{z} C_{33} z F_{wave}(t) ]其中 (m) 是船舶质量(A_{33}) 是垂荡附加质量(B_{33}) 是辐射阻尼系数(C_{33} \rho g A_{wp}) 是静水恢复力系数等于水线面积乘以水密度和重力加速度。这些系数在真正的工程设计中要由势流软件如 WAMIT、AQWA或模型试验给出但课程设计和概念验证阶段可以先用经验公式估算附加质量取船舶排水量的 0.30.8 倍阻尼比取临界阻尼的 5%15%。把二阶方程改写成状态空间形式是接 Simulink 最干净的方式A [0 1; -C33/(mA33) -B33/(mA33)]; B [0; 1/(mA33)]; C eye(2); % 输出位移和速度 D [0; 0]; sys ss(A, B, C, D);状态空间的好处是回头在 Simulink 里只需要一个State-Space模块不用手动搭积分器和增益模块出错的概率低很多。而且sys可以直接用线性系统分析工具箱做传递函数分析比如看垂荡的共振频率和阻尼比这些参数在调waveForce的黏性系数时有指导作用。4.2 Simulink 模型搭建与波浪力输入在 MATLAB 脚本里建模型的自动化写法如下。这样不依赖手动拖拽鼠标模型可以通过脚本反复重建参数修改后一键重新生成。mdl ship_motion_model; if bdIsLoaded(mdl), close_system(mdl, 0); end new_system(mdl); open_system(mdl); % 添加状态空间模块 add_block(simulink/Continuous/State-Space, [mdl /ShipMotion], ... Parameters, sys); % 添加波形输入: 从工作区读取波浪力时间序列 t 0:param.dt:param.Tsim; Fwave waveForceTimeSeries(t, param); % 预先算好整段时间的力 assignin(base, Fwave, Fwave); assignin(base, twave, t); add_block(simulink/Sources/From Workspace, [mdl /F_wave], ... VariableName, Fwave, SampleTime, param.dt); add_block(simulink/Sinks/To Workspace, [mdl /z_out], ... VariableName, z_out, SaveFormat, Timeseries); add_block(simulink/Sinks/Scope, [mdl /Scope]); add_line(mdl, F_wave/1, ShipMotion/1, autorouting, on); add_line(mdl, ShipMotion/1, z_out/1, autorouting, on); add_line(mdl, ShipMotion/1, Scope/1, autorouting, on); save_system(mdl); simOut sim(mdl, Solver, ode45, StopTime, param.Tsim); z simOut.get(z_out);这段脚本里有几个关键点。From Workspace模块默认要求输入是[t, u]两列或timeseries对象我这里传入的是两个单独数组SampleTime必须和param.dt保持一致否则 Simulink 会按默认步长插值导致波浪力高频成分被平滑掉。autorouting参数让连线自动布线适合代码方式建模型。仿真器选择ode45对四阶运动方程精度足够但如果电磁学的问题浮点发散可以换ode15s看是否是刚性系统。注意sim命令里的StopTime是字符串传字符串形式而不是数字不然老版本会报类型错误。4.3 脚本与 Simulink 联合仿真的参数表模块关键参数推荐设置说明State-SpaceA/B/C/D 矩阵由ss函数生成不要手敲矩阵权重在脚本里改From WorkspaceSampleTime等于param.dt不一致时波形被重采样Solver类型与步长ode45自动步长收敛慢时换 ode15s输入幅值Fwave先归一化再接入检查峰值是否达到 (10^6) 量级初始条件z(0), zdot(0)一般是 0若瞬态过长增加仿真预热段param.Tsim的选择直接决定波面序列长度。300 秒的仿真时长在典型 Tp8 秒海况下包含约 37 个波浪周期足以让垂荡响应进入稳态。waveForceTimeSeries是对第 3 章waveForce的封装内部循环时间步并调用切片法如果每步都算切片法 300 秒 3000 步就是 90 万次循环MATLAB 里跑起来有点慢建议先用parfor并行或者预计算压力分布表。5. main.m 串联运行与仿真发散排查的实用技巧main.m是整个项目的指挥中心。它的核心任务是把waveModelInit、linearWaveSimulation、waveForce、Simulink 模型按顺序串起来并且输出便于分析的位移、速度、加速度时序。我推荐在main.m里加一个runtime计时和断言检查能从源头拦截一半的问题clc; clear; close all; tStart tic; % 1. 初始化参数 param waveModelInit(PM, 4.0, 8.5); param.rho 1025; param.g 9.81; param.L 120; param.B 18; param.draft 6; param.Nslice 40; param.dx param.L / param.Nslice; param.Bm ones(1, param.Nslice) * param.B * 0.9; param.vertVel zeros(1, param.Nslice); % 先占位后面按线性波填充 % 2. 生成波面并预计算波浪力时间序列 Nt round(param.Tsim / param.dt); FwaveAll zeros(Nt, 1); for k 1:Nt tNow k * param.dt; etaNow linearWaveSimulation(param, tNow); % 沿船长分布 param.vertVel zeros(1, param.Nslice); % 用解析速度重新填充 FwaveAll(k) waveForce(etaNow, param); end % 3. 运行 Simulink 模型假设模型已存在 simOut sim(ship_motion_model, Solver, ode45, ... StopTime, num2str(param.Tsim)); z simOut.get(z_out); elapsed toc(tStart); fprintf(仿真完成, 耗时 %.1f s\n, elapsed);跑通之后先核对三件事。第一max(abs(z.Data))是否在合理区间4 米有义波高对应的垂荡单幅值不应超过 2.5 米超过说明恢复力系数或阻尼系数设置偏小。第二绘位移时域曲线观察起始阶段是否有衰退率过大的瞬态振荡正常辐射阻尼下前 1 到 2 个周期会有衰减之后进入与波频一致的稳态。第三按下 Simulink 的 Scope 波形和 MATLAB 的z时间序列对比看数据是否对齐防止采样率不一致的问题。如果发散按效率排序检查以下位置。最先看仿真步长是否覆盖了波频最大分量最高波浪频率 (3\omega_p) 对应周期约 2.8 秒固定步长必须小于周期的五十分之一否则仿真数值上天然不稳定。其次检查Fwave是否有异常尖峰波浪力在船首和船尾的切片可能出现淹没截面突变需要在切片区段做光滑处理。再次看A_{33}和B_{33}的经验系数是否与船舶尺度匹配120 米船吃水 6 米的垂荡固有周期大约在 69 秒如果算出来只有 2 秒就要怀疑刚度系数的水线面积量级少乘了船宽。最后检查波形输出模块是否默认限制了最大数据点改To Workspace的SaveFormat为Timeseries可以保存全精度结果。最后一个实用技巧验证时把随机相位种子param.seed轮流设成 0、1、2、3 跑四组然后把四组垂荡有效值做平均这样统计值更接近理论值这也是评审老师或项目验收比较认可的做法。本文还有配套的精品资源点击获取