简介这是一份面向海洋工程与船舶仿真学习者的MATLAB源码包聚焦船舶在波浪中的运动模拟涵盖线性波与随机波生成、波浪力计算及船体运动响应等核心环节。资源共7个文件其中6个为MATLAB的.m脚本分别承担波浪模型初始化、线性波模拟、PM谱随机波浪生成、波力计算与主程序运行等任务另有1个txt说明文档可辅助理解代码结构与参数设置。压缩包整体仅2KB轻量紧凑适合初学者对照源码逐行学习波浪理论在MATLAB中的实现。目前已有1555人学习下载。通过该资源读者可以掌握基于傅里叶级数的线性波构建方法了解JONSWAP或Pierson-Moskowitz谱的随机波浪模拟思路并理解波浪与船体相互作用的力学处理流程为后续船舶耐波性分析、Simulink联合仿真等研究提供基础。1. 船舶在波浪中的仿真matlab从运动方程到可运行程序把一艘船放进不规则波里用 MATLAB 看它怎么颠簸很多人第一反应是先画波面再算浮力。实际做下来会发现船舶在波浪中的运动响应大头在水动力记忆效应船体运动扰动流场后辐射力要延迟一段时间才衰减完不能把阻尼当常数直接乘速度。船舶在波浪中的仿真 matlab 常见路线是先用切片法或面元法求频域附加质量、阻尼和波浪力传递函数再用 Cummins 方程搬进时域状态方程求解。同一套频域数据可同时给出规则波 RAO 和不规则波时历计算量远小于 CFD。这篇文章写的是如何用 MATLAB 把船舶在波浪中的运动方程落到可运行脚本并把发散排查和参数边界讲透。适合正在做耐波性分析、想自己控制求解过程而不是只拖 Simulink 模块的人。2. 船舶在波浪中的运动模型怎么搭从六自由度到 Cummins 状态空间2.1 为什么先做垂荡和纵摇2 自由度运动方程船舶在波浪中的运动完整描述需要六个自由度纵荡、横荡、垂荡、横摇、纵摇、艏摇。如果把六个自由度都放进矩阵每个矩阵都是 6×6耦合项有几十项而且垂荡与纵摇在迎浪工况下基本主导垂向加速度、砰击和甲板上浪。所以第一版仿真从这两个自由度开始既能讲清楚模型又不用过早陷入横摇非线性和漂移力的泥潭。时域运动方程常用 Cummins 形式(M A(∞)) ξ″(t) ∫₀ᵗ K(τ) ξ′(t−τ)dτ Cξ(t) F_wave(t)其中 ξ [z; θ]z 是垂荡位移θ 是纵摇角M 是船体质量与纵摇转动惯量A(∞) 是频率趋于无穷时的附加质量矩阵C 是回复力矩阵F_wave(t) 是波浪激励力。与弹簧质量系统最大的差别是那个卷积积分它把过去所有时刻的速度都计入当前辐射力这就是波浪记忆效应。C 矩阵对垂荡和纵摇分别取C diag(ρgA_wp, ρg∇GM_L)。其中 A_wp 是水线面面积∇ 是排水体积GM_L 是纵稳心高。注意垂荡刚度单位是 N/m纵摇刚度单位是 N·m/rad两者差着量级写状态方程时最容易丢的就是这一处。符号物理意义常见获取方式M船体质量与转动惯量总布置图、装载手册A(∞)无穷频率附加质量切片法外推、BEM 高频点K(τ)脉冲响应函数频域阻尼的余弦变换C静水回复刚度ρgA_wp、ρg∇GM_LF_wave(t)波浪激励力Froude-Krylov 加绕射力这个表格是后续代码里每个矩阵的直接来源。对有经验的工程师值得注意的不只是 A(∞) 本身而是 A(∞) 与频域附加质量 A(ω) 的关系BEM 软件只算到有限频率A(∞) 需要做高频渐近外推外推误差会在卷积积分前段造成尖峰后面第四章会专门讲这个现象的排查。2.2 用 MATLAB 把频域阻尼 B(ω) 转成时域记忆函数 K(τ)频域阻尼 B(ω) 在每个频率点是一个 2×2 矩阵。按 Cummins 理论K(τ) 是 B(ω) 的余弦逆变换K(τ) 2/π ∫₀^∞ B(ω) cos(ωτ) dω实际数据里 ω 只能给到某个最大值积分必须截断。经验约束是 ω_max 至少取到遭遇频率上限的 2 到 3 倍τ 的采样长度要覆盖 K(τ) 衰减到零的时间。船体越大衰减越慢100 米船长量级上我一般取 τ_max 60 秒步长 0.05 秒。function K computeRetardation(omega, B, tau) % omega: 1xNw 频率向量rad/s % B: Nw x 2 x 2 频域阻尼矩阵 % tau: 1xNt 延迟时间向量s Nw numel(omega); Nt numel(tau); K zeros(Nt, 2, 2); B_mat squeeze(B); % 转为 Nw x 2 x 2 for it 1:Nt % R2016b 之后自动隐式扩展Nw x1 余弦乘到每个矩阵元素 integrand B_mat .* cos(omega(:) * tau(it)); K(it,:,:) (2/pi) * trapz(omega(:), integrand, 1); end end这段代码里最关键的是trapz(omega(:), integrand, 1)它沿频率轴对 2×2 里每个元素分别积分。omega(:) * tau(it)得到 Nw×1 向量会和 Nw×2×2 的B_mat自动广播。如果你还在用 R2016a 之前的版本需要先用repmat把余弦扩展到同样尺寸。拿到 K(τ) 后先画曲线正常的 K(τ) 在 τ0 处有有限峰值随后带阻尼振荡并趋近于零。如果曲线出现锯齿说明频域阻尼点数不够或 ω 采样不均匀。另一个常见问题是切片程序导出的 B(ω) 频率点不等距建议在频域线性插值后再做变换。K(τ) 直接参与时域积分它的数值质量决定仿真能不能不发散。2.3 把方程写成 MATLAB 状态空间形式时域积分需要把二阶方程降成一阶。取状态向量 x [z; θ; ż; θ̇]方程写成x′ [0 I; −M_total⁻¹C, −M_total⁻¹D] x [0; M_total⁻¹] F_wave其中 D 是常数阻尼项。严格做法是保留卷积项并把 K(τ) 拟合成状态空间用一组指数函数近似然后增加内部状态 μ这样每次计算不必从 0 积分到当前时刻效率高也方便放进 Simulink。M_total M_rigid A_inf; detM det(M_total); if detM 0 error(A_inf 导致质量矩阵非正定检查高频外推值); end A_state [zeros(2), eye(2); - (M_total \ C), - (M_total \ B_lin)]; B_state [zeros(2); M_total \ eye(2)];这里先说矩阵构造再用反斜杠求解不要显式写inv乘。虽然 2×2 差别不大但六自由度时数值差异会非常明显。B_lin是从频域阻尼 B(ω) 里取遭遇频率附近的对角项如果完全省略这一项系统变成无阻尼共振仿真结果会持续振荡而不是收敛这是新手最常调不出来的原因。3. 用 MATLAB 写船舶在波浪中的仿真最小脚本与参数设置3.1 输入参数表与波谱生成写可复现脚本之前先把输入参数定下来。下面是一组教学用假想船参数量级相当于 100 米级散货船不是任何真实船型的精确数据。参数数值单位说明船长 L105m垂线间长排水量 Δ7500t对应质量 m 7.5e6 kg纵摇惯性半径 kyy26mIyy m·kyy²水线面面积 Awp1580m²垂荡回复刚度来源纵稳心高 GM_L128m纵摇回复刚度来源航速 U6m/s约 12 kn浪向角 β180deg迎浪cosβ −1这里最容易带偏的参数是纵摇惯性半径。Iyy 不要用船长平方直接乘质量一般取 (0.22~0.28)L 范围散货船满载偏大、压载偏小。GM_L 从完整稳性计算书拿不要自己从型线图估。量级对了后续调参才有意义。不规则波力按 JONSWAP 谱离散每个频率分量幅值 aᵢ √(2S(ωᵢ)Δω)相位随机。迎浪时遭遇频率 ω_e ω − ω²U/g·cosβ所有力、力矩传递函数都要插值到遭遇频率点上。相位随机意味着每次运行结果不同做统计时至少跑 20 组随机种子。3.2 可直接运行的规则波仿真脚本RK4 最小实现下面代码先用常数阻尼代替卷积记忆项把规则波下的垂荡和纵摇跑通。这样做是为了先建立可运行的基线再在这个基线上替换记忆效应。% ship_sim_wave_linear.m % 2自由度垂荡纵摇规则波常数阻尼RK4 主循环 clear; clc; close all; % 船体参数教学用假想船 m 7.5e6; % 排水质量kg Iyy m * 26^2; % 纵摇转动惯量kg*m^2 rho 1025; g 9.81; Awp 1580; % 水线面面积m^2 GM_L 128; % 纵稳心高m M_rigid diag([m, Iyy]); A_inf diag([0.15*m, 0.08*Iyy]); % 高频附加质量示意值 M_total M_rigid A_inf; C diag([rho*g*Awp, rho*g*Awp*GM_L]); omega_e 0.63; % 遭遇频率rad/s约 10s 周期 B_lin diag([8e4, 2.5e7]); % 频域阻尼在遭遇频率处的折算 D_lin diag([1e4, 2e6]); % 附加线性粘性阻尼 A_state [zeros(2), eye(2); - (M_total \ C), - (M_total \ (B_lin D_lin))]; B_state [zeros(2); M_total \ eye(2)]; % 规则波激励波幅 1m力幅由切片法给出 A_wave 1.0; F3_amp 1.2e6; % 垂荡力幅N Q5_amp 3.0e7; % 纵摇力矩幅N*m dt 0.1; tMax 100; t 0:dt:tMax; N numel(t); Fwave [F3_amp * A_wave * sin(omega_e*t); Q5_amp * A_wave * sin(omega_e*t)]; % RK4 主循环 x zeros(4,1); his zeros(N,4); for i 1:N-1 u Fwave(:,i); k1 A_state * x B_state * u; k2 A_state * (x 0.5*dt*k1) B_state * u; k3 A_state * (x 0.5*dt*k2) B_state * u; k4 A_state * (x dt*k3) B_state * u; x x dt/6*(k1 2*k2 2*k3 k4); his(i1,:) x; end plot(t, his(:,1), t, his(:,2)*10); legend(heave/m, pitch*10/rad); xlabel(t/s); ylabel(response);这段代码可以直接保存为脚本运行。A_inf用的是高频附加质量的示意值不是波频处的 A(ω)这一点不要混用。C中纵摇刚度是 ρgAwp·GM_L量级约 2e9 N·m/rad垂荡刚度约 1.6e7 N/m差两个量级但在状态矩阵里通过左除已经做了归一化不会导致数值问题。RK4 四个 stage 共用同一个u即当前时间步的波浪力。对于线性系统且波浪力是缓变信号时这个写法就是标准 RK4如果后面把辐射力记忆项加进来就不能再四个 stage 共用同一个卷积外力那时要么用状态空间近似要么把记忆项冻结在主步开始下一节给替换位置。3.3 从常数阻尼切换到卷积记忆项要严格保留记忆效应需要把 A_state 里的 B_lin 去掉改为每步计算辐射力 F_rad。替换的核心表达式是% 预计算 K33、K35、K53、K55方法见 2.2 的 computeRetardation % 主循环内需要保存过去每个时刻的速度 vel_hist(1:i, :) % Frad(1) trapz(tau, K33.*vel_z_interp K35.*vel_p_interp); % Frad(2) trapz(tau, K53.*vel_z_interp K55.*vel_p_interp); % F Fwave(:,i) - Frad;这里vel_z_interp和vel_p_interp是在时间轴 τ 上插值出的滞后速度v(t−τ)。每次积分都要从 0 到当前时刻重新算复杂度是 O(N²)仿真 2000 步还能接受再长就会明显变慢。工程上更可靠的做法是把 K(τ) 拟合成四五阶传递函数用附加状态 μ 代替卷积历史这属于进阶优化。提示在把卷积项写进主循环之前先在零波浪力工况跑一次自由衰减。K(τ) 首个峰值应为正若整体偏负回头检查 B(ω) 矩阵是否传反了转置方向。4. 仿真发散排查船舶在波浪中的运动数值稳定性与控制4.1 发散现象与第一排查顺序船舶在波浪中的运动仿真发散最常见的表现是位移在几十步内跳到 1e15、出现 NaN或者位移一直增长但数值没有溢出。前者说明动态不稳定后者更像静力失稳或缺回复刚度。拿到发散曲线不要急着调小 dt先按顺序查状态矩阵特征值、质量矩阵正定性、阻尼是否有正贡献、波浪力幅值是否过大。确定线性部分是否稳定用一行代码就能定位eigA eig(A_state); if any(real(eigA) 1e-6) fprintf(A_state 不稳定最大实部 %.3e\n, max(real(eigA))); end如果 A_state 本身有不稳定特征值任何步长下都会发散。这时重点检查 C 矩阵是否有漏项纵摇刚度写 ρgAwp·GM_L 时 GM_L 必须用米垂荡刚度如果顺手写成 ρg∇虽然量级接近但物理不对会改变特征频率。4.2 步长、阻尼和卷积截断带来的振荡排除结构失稳后还发散就是数值原因。RK4 的稳定域比欧拉大但遇到波浪高频分量时过大的 dt 会让高频成分进入不稳定区时历上表现为明显锯齿。有效做法是看 FFT 频谱是否在奈奎斯特频率附近出现折叠峰。卷积截断造成的发散和步长发散不同它通常在仿真开始后几十秒才出现。K(τ) 只截到 30 秒而水动力记忆还没衰减完相当于每个波浪周期都给系统注入一点能量。排查方法是把 K(τ) 画出来看末端是否归零如果没有延长 τ_max 并检查 ω_max 是否覆盖阻尼峰值。发散现象可能原因优先处理前 20 步直接 NaNA_inf 导致质量矩阵奇异检查 det(M_total)位移指数增长C 矩阵缺失或符号错检查 GM_L 单位和 C 正定性高频振荡叠加dt 超过波频限制dt 降到 π/ω_max 以下延迟到 50s 后增长卷积截断过短τ_max 加倍纵摇发散但垂荡正常辐射阻尼非对角项被忽略补 B35/B534.3 抑制自由衰减与初始化技巧把初始速度设为零时船体会产生一组自由衰减响应叠加在波浪响应上如果波浪力从 t0 突然加载还会产生虚假冲击。常见做法是给波浪力前 5 秒加平滑窗ramp sin(pi/2 * min(t/5, 1)).^2; Fwave(:, i) Fwave_full(:, i) * ramp;ramp 在 0 到 5 秒内从 0 平滑到 1保留波浪力频率成分但避免阶跃冲击。另一个值得养成的习惯是给阻尼矩阵加一个下限防止某些频率点上 B(ω) 接近零造成伪共振峰。下限取最大阻尼的 1%不会明显改变耐波性统计结果却能避免数值上出现不合理的窄带共振。5. 船舶在波浪中的运动仿真验证从垂荡纵摇扩展到六自由度5.1 用规则波 RAO 校验仿真脚本时域脚本跑通后先验证再算统计值。最省事的验证是输入一个频率固定、波幅为 a 的规则波运行 30 个波浪周期取后 20 个周期的垂荡或纵摇幅值除以波幅得到 RAO然后与频域解直接比较。% 规则波验证取稳态段幅值 z_steady his(t 40, 1); RAO_sim (max(z_steady) - min(z_steady)) / (2 * A_wave); RAO_fd abs(F3_amp / (C(1,1) - (m A_inf(1,1)) * omega_e^2)); % 示意更精确的做法是用 findpeaks 对后 20 个完整周期的峰值取平均再除以波幅。若 RAO 误差超过 5%不要先调阻尼回头检查遭遇频率映射很多时历脚本的幅值谱和相位谱都对但把遭遇频率当成了自然频率导致共振峰偏移。5.2 从 2 自由度扩展到六自由度时的工程取舍扩展时矩阵从 2×2 变成 6×6A(∞)、B(ω)、C、F_wave 全部变大。最需要注意的是横摇回复刚度 ρg∇GM_T 与垂荡纵摇刚度量级差异横摇刚度通常只有垂荡刚度的百分之一不到耦合进状态方程后容易产生刚性问题。固定步长 RK4 此时效率明显下降建议改用 ode45 或 Simulink 自适应求解如果系统同时存在大刚度和波浪高频则要上隐式求解器避免时步被刚性限制拖到很小。另一个细节是横荡、艏摇与垂荡纵摇的弱耦合。斜浪里横荡与艏摇通过 B24/B42 与横摇耦合忽略这些项会让横摇 RAO 低估正横浪下零航速时耦合仍然存在。也就是说六自由度的耦合项选择依赖浪向不是为了全而全直接对角化会漏掉共振频率。5.3 用 Simulink 做联合仿真时如何保留记忆函数Simulink 里搭船舶在波浪中的运动模型时卷积记忆项不能直接用连续积分模块因为积分上限是当前仿真时间变步长求解器每步都可能变化。我一般先把 K(τ) 拟合成 4 到 8 阶传递函数放入ss对象作为线性系统与六自由度运动方程并联。拟合可以用tfest或对 K(τ) 做 Prony 近似误差控制在频域阻尼重构的 2% 以内耐波性结果差异通常可忽略。如果仿真只用于评估垂向加速度和砰击输出端加一个低通滤波器截止频率取 1.5 倍最大波浪频率。后续把这段运动时历接到结构载荷计算里时垂向加速度输出也要用同一截止频率滤波否则砰击载荷峰值统计会比实际偏高。本文还有配套的精品资源点击获取