简介面向海洋工程与船舶设计研究者的 MATLAB 仿真资源包聚焦船舶在风浪流作用下的六自由度动态响应与运动性能分析。压缩包共 8 个文件包含 5 个 Simulink 模型、2 个 MATLAB 脚本及 1 个数据说明文档总大小约 27KB结构紧凑。模型文件覆盖船舶动力学方程、神经网络预测模块和控制系统可搭建从外部环境激励、船舶运动响应到控制器调节的完整仿真链路脚本分别用于神经网络训练与 BFGS 参数优化配套文本提供风浪数据和初始工况输入。借助这套资源读者能学习船舶运动数学模型的 Simulink 实现方法、参数辨识与调优思路并将神经网络、优化算法应用于运动预测与控制器设计。已有 1528 人浏览学习适合具备一定 MATLAB 基础、希望系统开展船舶海上运动仿真研究的海洋工程与船舶工程技术人员。1. 船舶海上运动 matlab仿真别急着搭六自由度先看清这个仿真在解决什么船舶海上运动 matlab仿真是船舶水动力与耐波性方向最常见的落地场景之一它把船体在波浪、风、流作用下的摇荡运动用微分方程描述再交给 MATLAB 做时间域求解。你最终得到的是横摇、纵摇、垂荡这些运动量在一条条海况下的时历曲线以及对应的统计幅值。它带给你的直接价值是“这条船在这个海况下到底晃多大”“改船型参数有没有效果”。这个方向我建议新人先做垂荡—纵摇二维模型而不是一上来就六自由度参数少、物理意义直观、还方便和理论解核对把坐标系、运动方程和数值积分这三件事练扎实再往三维扩展就顺了。2. 船舶运动建模坐标系、六自由度基础与垂荡-纵摇降阶2.1 惯性系与船体坐标系角度符号定死后面才不乱做船舶仿真第一件事不是写方程而是把坐标系写清楚。MATLAB 里随意定义坐标轴短期没问题一旦开始加风浪方向、做控制闭环符号冲突的坑就全冒出来了。我将默认采用行业惯例惯性系大地坐标系原点取静水面X 轴指向船艏方向Z 轴向上为正船体坐标系固连在船体上原点取重心或水线面中心附近三个轴随船转动。六个自由度的位置量用惯性系描述速度量平移速度和转动角速度默认放在船体坐标系里写。两个坐标系之间用欧拉角旋转矩阵转换常见的顺序是 Z-Y-X也就是先艏摇、再纵摇、最后横摇。注意这里的“Z-Y-X 顺序”是旋转先后顺序不是坐标轴位置顺序写矩阵时别和三维图形学里的习惯混用。初学者最常犯的错是把横摇角和纵摇角对调导致仿真结果里横摇幅值跑到纵摇曲线上我见过不止一次查了半天最后发现是矩阵列写反了。提示把角度正方向写进注释里例如“横摇角 φ 右舷下沉为正”。如果你是左舷下沉为正整套作用力符号都要跟着翻转宁可在代码开头定清楚不要在 300 秒仿真之后对着结果猜。2.2 六自由度运动方程骨架惯性、附加质量、阻尼、恢复与激励刚体在波浪中的六自由度运动简化后可以写成如下形式(M A(ω)) ·v_dot τ_rest τ_damp τ_wave τ_wind/current其中 M 是船体质量惯性矩阵A(ω) 是附加质量矩阵它代表船体推动周围水一起运动时等效出来的质量。在实海域仿真中A(ω) 随波浪频率变化在大多数入门实现里用一个在其固有频率附近取值的常数矩阵就可以起步。右侧分别是恢复力矩、阻尼力矩、波浪激励力和风/流干扰力。恢复力来自浮力和重力不平衡。垂荡的恢复力系数在静水条件下直接等于 ρ·g·A_wp其中 A_wp 是水线面面积。纵摇的恢复力矩来自水线面面积在纵向的重新分布严格讲要用水线面面积的纵向二次矩 I_wp_yy 乘以 ρ·g 得到系数。阻尼力则复杂得多包含兴波阻尼、涡流阻尼和摩擦阻尼。简化建模时常用线性阻尼系数 b33、b55并调到一个能稳定收敛的量级。真实的阻尼特性还要附加上辐射阻尼随频率的变化那属于势流理论的内容初版仿真不需要。初版仿真想让结果有意义核心是把频率特性对准恢复力系数决定固有周期阻尼系数决定共振峰高度附加质量修正固有周期数值。三者配合不对仿真画出来的时历曲线长得很像回事但共振点偏移一两秒你怎么调都无法和试验值匹配。所以建模之初不要盲目加复杂项先把这三类参数的物理意义和单位同步标定好。2.3 垂荡-纵摇降阶纵向平面能解决 80% 的耐波性评估问题把六自由度砍成垂荡-纵摇二维模型的理由很简单多数排水型船在迎浪和斜浪工况下纵向平面运动最容易造成甲板上浪、螺旋桨飞车和垂向加速度超标横摇虽然危及安全但它主要靠横摇周期附近的阻尼和减摇鳍来解决物理机制相对独立。初版仿真把纵向平面做对再横向扩展横摇每一步的结果都能单独验证。垂荡-纵摇耦合方程可以写成矩阵形式[[mA33, A35], [A53, IyyA55]] · [ẍ₁; ẍ₂] [[b33, b35], [b53, b55]] · [ẋ₁; ẋ₂] [[c33, c35], [c53, c55]] · [x₁; x₂] [F_wave_z; F_waveθ]这里的 x₁ 是垂荡位移x₂ 是纵摇角。很多船型的 A35、c35 耦合项不大初版可取零但阻尼耦合项 b35 在某些高速船上相对明显如果需要算垂荡加速度和纵摇角谁影响谁建议保留。取耦合项之前先跑一个不带耦合的版本用垂荡固有周期和纵摇固有周期分别校验再决定加不加。耦合项一开始就打开出了问题很难定位。2.4 把二阶方程改写成 ode45 可解的一阶状态方程MATLAB 的 ode45 只能处理一阶常微分方程组。二阶运动方程必须引入状态变量才进得了求解器。我惯用的状态定义是x [z; theta; z_dot; theta_dot]垂荡位移、纵摇角、垂荡速度、纵摇角速度那么状态方程写作x_dot(1) x(3)x_dot(2) x(4)x_dot(3) inv(有效质量矩阵) · (恢复项 阻尼项 波浪力) 的第 1 行x_dot(4) 同样的第 2 行这段逻辑是整篇文章的核心。后面所有加减波浪谱、风载荷甚至控制闭环都是在扩展这里右侧的力项而不动数值求解框架。我在第 3 章给出可直接运行的 MATLAB 示例。3. 用 MATLAB 脚本跑通最小垂荡-纵摇仿真完整代码与参数3.1 船型参数与初始条件先把物理量设对再跑数值先给出一个虚拟的小型巡逻船尺度参数适合做方法验证数值不必对应真实船但量纲必须自洽排水量 120 吨垂荡固有周期约 7.2 秒纵摇固有周期约 6.5 秒水线面面积 180 平方米。附加质量系数垂荡取 0.8 倍排水质量纵摇取 0.4 倍转动惯量。阻尼比垂荡取临界阻尼的 10%纵摇取 12%这是实海域数据常见的量级。波浪激励先用一个圆频率等于 0.85 rad/s 的规则正弦力幅值垂荡力 60 kN纵摇力矩 220 kN·m。% 最小垂荡-纵摇仿真: 主脚本 % 状态 x [z, theta, zdot, theta_dot] % 单位制: SI (m, kg, s, rad) rho 1025; % 海水密度 kg/m^3 g 9.81; % 重力加速度 m/s^2 m 120e3; % 排水量 kg, 对应120吨 Awp 180; % 水线面面积 m^2 Iyy 3.2e7; % 纵摇转动惯量 kg*m^2, 示意值 % 附加质量 A33 0.8 * m; A55 0.4 * Iyy; % 恢复力系数 c33 rho * g * Awp; % N/m c55 rho * g * (Awp * 30); % 用纵向水线二次矩的简单近似, 单位 N*m/rad % 阻尼系数: 按临界阻尼比近似 omega_z sqrt(c33 / (m A33)); omega_t sqrt(c55 / (Iyy A55)); b33 2 * 0.10 * (m A33) * omega_z; b55 2 * 0.12 * (Iyy A55) * omega_t; % 波浪激励: 初始用规则波正弦近似 wave_freq 0.85; % rad/s Fz0 60e3; % 垂荡力幅值 N Ftheta 220e3; % 纵摇力矩幅值 N*m参数说明阻尼系数 b33、b55 在这里按“临界阻尼比”回推这样比直接拍脑袋给一个数要可靠。临界阻尼比 0.100.15 是排水型船纵向运动的常见范围。恢复系数 c55 我用的是水线面面积的纵向二次矩简化式严格算要积分水线面几何初版这样近似已经能让固有周期落在合理区间。3.2 动力学函数定义给 ode45 一个干净的右端函数function xdot ship_longitudinal_ode(t, x, Mmat, Bmat, Cmat, F) % 提取状态 z x(1); theta x(2); zdot x(3); q x(4); % 纵摇角速度 % 波浪力按时间变化 wave_freq F.wave_freq; Fz F.Fz0 * sin(wave_freq * t); Ft F.Ftheta * sin(wave_freq * t); % 恢复力和阻尼力 tau_rest -Cmat * [z; theta]; tau_damp -Bmat * [zdot; q]; tau_wave [Fz; Ft]; % 求解加速度: acc inv(Mmat) * (tau_rest tau_damp tau_wave) acc Mmat \ (tau_rest tau_damp tau_wave); xdot [zdot; q; acc(1); acc(2)]; end逻辑说明右端函数返回四维列向量 [速度加速度]。恢复力项和阻尼项我特意写成负号因为它们是反抗运动的力。用“\”运算符做矩阵求逆MATLAB 会采用数值稳定的求解方式比 inv(Mmat)*force 更稳也更快。波浪力这里强制只由时间决定和当前船体位置无关这就是所谓的“弗劳德-克雷洛夫力 入射波力”最简化近似等船体运动幅值大起来以后这种近似会误差大但入门阶段够用。3.3 主程序里调 ode45 并画图看时历曲线验证合理性% 组装质量/阻尼/恢复矩阵 Mmat [m A33, 0; 0, Iyy A55]; Bmat [b33, 0; 0, b55]; Cmat [c33, 0; 0, c55]; % 输入结构体 F.wave_freq 0.85; F.Fz0 60e3; F.Ftheta 220e3; t_span [0, 300]; x0 [0; 0; 0; 0]; % 从静水面静止出发 % 求解 [t, x] ode45((t, x) ship_longitudinal_ode(t, x, Mmat, Bmat, Cmat, F), ... t_span, x0); % 画图 figure; subplot(2,1,1); plot(t, x(:,1) * 100); % 垂荡, 以cm显示更直观 ylabel(垂荡 (cm)); grid on; subplot(2,1,2); plot(t, x(:,2) * 180 / pi); % 纵摇角, 转换为度 ylabel(纵摇 (deg)); xlabel(时间 (s)); grid on;参数说明x0 从零开始前几十秒会有瞬态过程这属于初始激励的响应分析稳态幅值时记得跳过前 100 秒再取极值。我通常把仿真时长设到 300 秒以上至少覆盖 20 个固有周期否则统计幅值没有代表性。垂荡位移单位换算成厘米显示纵摇角换算成度显示是为了人工读图时有一个直观的运动幅度概念。3.4 调参规律把激励频率扫过固有频率看共振是否出现这个最小模型跑通以后第一件值得做的事是扫频。把 wave_freq 从 0.4 到 1.4 rad/s 变化每个频率跑 300 秒记录稳态垂荡幅值和纵摇幅值画成幅频响应曲线。你会看到在 0.87 rad/s 左右幅值明显抬升这就是垂荡固有频率理论值可以由 sqrt(c33/(mA33)) 算出来做对比。纵摇固有频率略低于垂荡频两条峰值错开是正常现象。若峰值对应的频率和理论固有周期差得超过 15%先查附加质量和恢复系数有没有写错。这里我要强调一个初学者常误判的点波浪激励频率并不等于船体运动频率。船体运动是强迫振动在激励频率附近振荡但当激励接近固有频率时运动幅值显著放大所以“共振”是运动对激励的响应不是说船的摇荡周期会变成波浪周期。搞清楚这一层后面看 P-M 谱、JONSWAP 谱出来的大振幅时间区段就不会误以为仿真出 bug。4. 把不规则波装进模型海况谱选择与波浪力加载4.1 规则波只能做验证做设计要上不规则波第 3 章的正弦波叫规则波实验室里造波机常用来标定船舶模型。但真正海面上的波浪是随机过程波高和周期都在变化所以实海域响应用谱分析方法把海浪看作无数频率规则波的线性叠加每个频率给一个幅值幅值大小由波谱密度函数决定。行业内最常见的两个谱是 Pierson-Moskowitz 谱P-M 谱成熟外海限风区长浪和 JONSWAP 谱有限风区、成长中的波浪。P-M 谱只有一个参数适用于规范设计海况JONSWAP 谱多了峰值增强因子 γ能描述波浪能量更集中的海况。选谱的标准很直接你手里的环境资料能给出什么参数。如果只知道有效波高 H_s 和平均跨零周期 T_z用 P-M 谱如果是针对某海域实测数据拟合的波浪谱多半是 JONSWAP 改进型。下面的代码用 P-M 谱生成波面时间序列和对应的垂荡激励力。4.2 用 P-M 谱生成波面时间序列谐波叠加法代码function [t_wave, wave_elev, Fz_wave] generate_pm_wave(Hs, Tz, T_total, dt) % P-M 谱: Hs 有效波高(m), Tz 平均跨零周期(s) % 输出波面时历(m)和简化垂荡力时历(N) g 9.81; omega0 2 * pi / Tz; % 谱峰圆频率近似, rad/s domega 0.005; % 频率分辨率 % omega 从 0.1 到 3 * omega0, 超出这个范围的能量极小可忽略 omega 0.1:domega:3 * omega0; % P-M 谱密度, 单位 m^2*s/rad % 常数 A 0.0081*g^2, B 0.74*(g/(2*pi*Tz))^4 A_pm 0.0081 * g^2; B_pm 0.74 * (g / (2 * pi * Tz))^4; S_omega A_pm ./ omega.^5 .* exp(-B_pm ./ omega.^4); % 随机相位 rng(42); % 固定随机种子, 结果可复现 phase 2 * pi * rand(size(omega)); % 谐波叠加合成波面 t_wave (0:dt:T_total); wave_elev zeros(size(t_wave)); for k 1:length(omega) ak sqrt(2 * S_omega(k) * domega); wave_elev wave_elev ak * cos(omega(k) * t_wave phase(k)); end end参数说明domega 是频率离散间隔取 0.005 能保证高频段和谱峰附近采样够密更小的 domega 会让波面更光滑但循环次数线性增加不需要刻意追求极小值。色谱峰 omega0 按 2π/Tz 近似P-M 谱的实际峰频率比这个值略低但用在力谱叠加里差异不大。固定 rng(42) 是为了让每位读者跑出来结果一致方便对比如果你要做蒙特卡洛统计把这个种子删掉多跑几十次再统计响应极值。4.3 把波面转成船体激励力从静水压力项到简化力不规则波产生的船体波浪力不能用波面直接乘一个系数。真正做法是计算入射波在船体湿表面上的压力积分这需要边界元或切片理论。MATLAB 里做初版工程评估常见做法是先用切片理论近似把船体沿船长切成若干切片每个切片当一段二维浮体求解单位幅波下的垂向力和纵向力矩响应再和波面时历做频域叠加。% 把波面激励转入运动方程: 垂荡力近似幅值响应 % 简化处理: 垂荡力幅值 rho * g * Awp * eta_a * 载荷系数 % eta_a 是切片段处的波面幅值, 载荷系数取 0.9 (典型排水型巡逻船) Hs 2.5; Tz 6.8; T_total 600; dt 0.1; [t_wave, wave_elev, ~] generate_pm_wave(Hs, Tz, T_total, dt); % 与船体吃水附近波高关联: 垂荡力近似 Awp 180; rho 1025; g 9.81; coef_load 0.9; Fz_irr rho * g * Awp * coef_load * wave_elev; % 更新动力学函数里 wave 的输入, 改为查表插值 % 在 ship_longitudinal_ode 中改为: % Fz_now interp1(t_wave, Fz_irr, t, linear, 0);逻辑说明interp1 做线性插值是因为 ode45 的时间步 t 不是固定 dt离散的力表需要插值才能和求解器对齐。注意外插选项设为 0防止仿真时间超过力表末尾时出现 NaN 或边界跳跃。实际情况里波浪力随频率的传递函数变化明显这个 coef_load 常数只会产生尺度偏差不会改变共振位置因此做“定性选型”足够做“定量极值”还得上水动力软件算 RAO幅值响应算子。我通常在初版仿真里把这一步的误差心理预期定在 ±25%。4.4 有义波高与周期的选择设计海况不是拍脑袋船舶耐波性评估常用的有义波高 H_s 等级H_s ≈ 2 米算轻浪4 米算中浪6 米以上算恶劣海况。设计时取哪一档要看你的船跑什么航线、作业窗口限制在哪。一个容易忽略的点是 H_s 与 T_z 强相关不能只调 H_s 不调周期。经验关系里 T_z ≈ sqrt(H_s) 的系数在 3.54.5 之间短宽型船取低值长瘦型船取高值。仿真时如果把二者配错谱峰频率会撞到纵向固有频率上算出来的极值比真实情况严重得多那就是自己吓自己。提示做谱分析时波面对应的是能量谱船体运动响应对应的是响应谱。把激励谱峰值和固有频率错开运动幅值就上不去刻意让谱峰对准固有频率算出来的是该海况最坏情况。两者都有用但写结论时一定要说清“我这个仿真取的是哪种工况”。5. 船舶 matlab 仿真避坑发散、步长、单位与坐标系5.1 仿真发散曲线在 50 秒附近突然飞上天现象垂荡位移或纵摇角按时历推进前半段看起来正常到某一时刻突然变成单调递增的抛物线状数值迅速到 1e10 量级图已经没法看。原因最常见的不是物理建模错误而是数值积分引入的刚性问题。你的质量矩阵里加了附加质量之后垂荡和纵摇的特征频率变化不大但耦合项、阻尼项的量级差异可能跨越三个数量级MATLAB 会在一台电脑上自动缩小步长……美中不足的是步长缩小到一定程度浮点误差累积反而又导致大跳变。解决先试 odeset 里 MaxStep 限制步长我一般在固有周期 T_n 的 1/50 以内例如 T_n 7.2sMaxStep 取 0.1s。把这个参数写进 opts odeset(MaxStep, 0.1, RelTol, 1e-6, AbsTol, 1e-7)。改完基本能压制到 10 分钟仿真不发散如果再发散再考虑换 ode15s 刚性求解器。每次都性能优化的诀窍是先固定随机种子舵偏转时要留出用于区分初值困难的余量。5.2 数值步长与固有周期的匹配MaxStep 到底设多大现象仿真能跑但垂荡时历看起来是锯齿状尤其是在波浪峰顶出现明显尖点。原因步长相对于固有周期太大了。典型值是 0.5 秒一个步长而垂荡固有周期才 7 秒一个周期只有 14 个采样点对正弦类响应来说分辨率明显不够。要认真讨论峰值的幅度再现一个周期至少要对齐 3050 个点。解决我一般以最小的固有周期为基准取 T_min / 40 做初始 MaxStep。如果不确定 T_min 是多少直接看恢复系数和等效质量算出的固有周期更保险的是用 ode45 默认自动步长然后对比 MaxStep 从 0.5 改成 0.1 时稳态幅值变化在 2% 以内这时视为收敛。这个“两步法”是我每次上仿真必做的体检。5.3 单位不统一用 kN 的力配了 kg 的质量现象垂荡位移漂移纵摇角看起来偏小但时历形状又有规律。原因力从外部资料里复制写着 60 kN代码里写成 60000本来是 N而质量用了吨一来二去差 1000 倍。另一个常见是周期单位外部数据是秒谱分析里圆频率用了 Hz整个激励频率偏出共振区。解决一手代码只选一套单位制。我习惯全程 SI质量用 kg力用 N长度用 m角度用 rad。凡是外部公式降下来的一次都在代码开头统一换算再进入后续计算。尤其注意角速度显示时转成度每秒计算时保留弧度每秒输出图层再做即时转换不要在状态量里存角度制。5.4 角度单位混用纵摇角一会儿 rad 一会儿 degree现象纵摇角时历看起来像是被压缩了几十倍还带明显的相位偏移。原因解析表达式里用了 sin(omega*t 相位)而 phase 看是红窗来自 deg2rad 的转换没做或者谱生成代码里 omega 是 rad/s而输入 T_z 给了秒结果相位和幅值全错位。解决写 ODE 函数第一行加一行长注释说明所有角度单位是 rad。然后在状态初值、波浪力相位、绘图前每一处转换都明确写成 180/pi 或 pi/180 系数。不要相信“我记得刚才已经转换过”把转换集中在两个位置做避免散布在代码各处。我曾在几个不同版本的脚本里三次转换相位结果偏偏又转回了原值这种错误浪费一下午排查。5.5 初始瞬态误判为共振现象300 秒仿真的前 30 秒垂荡幅值猛增到峰值之后慢慢回落初学者误以为数值发散。原因初值从静止出发相当于给系统一个阶跃激励响应包含自由衰减振动和强迫振动两部分。初始几秒的自由振动幅值由初值决定可能比稳态幅值大好几倍它是物理现象不是 bug。解决稳态分析时只取仿真后 50% 的时间段或者先设置一段衰减期。我通常在统计幅值时用代码忽略前 100 秒并且额外检查一个更长的仿真看稳态幅值是否与初始段无关。如果是两种不同初始条件的仿真稳态结果不能显著依赖初始条件否则说明阻尼太弱或者数值有奇异解。5.6 中文注释乱码与脚本编码现象MATLAB 2026b 打开别人发来的 .m 文件注释里全是乱码甚至报错。原因脚本编码不是 UTF-8MATLAB 新版本默认按 UTF-8 解析旧版 GBK 编码的中文注释在中文版 Windows 上有历史遗留问题。解决统一把脚本保存为 UTF-8。用 2026b 建议直接在“预设项-编辑器-语言”里把文件编码改成 UTF-8。代码正文本身不用中文也可以但注释里如果有错别字提示解释不清排查时极容易误导人。另外注意如果你的仿真脚本要在 Linux 服务器上用 MATLAB 跑批UTF-8 几乎是唯一不踩坑的选择。6. 验证结果再上闭环我的检验习惯和进阶方向仿真跑通只是起点真正可信的判断来自验证。我最常用的验证方法是对齐共振频率把激励谱在 0.31.5 rad/s 频段扫描记录垂荡稳态幅值找出峰值对应的频率再和理论固有频率 sqrt(c33/(mA33)) 对比。偏差在 10% 以内说明恢复力和附加质量建模基本正确偏差很大请回头先查附加质量或恢复系数的量纲这个动作能筛掉至少一半的隐藏错误。第二层验证是能量检查给一根裸船不加波设定初始垂荡位移 0.5 m观察 500 秒自由衰减看振动周期和阻尼衰减率是否符合设计值时同时看系统有没有额外能量进来。如果自由衰减曲线出现“越衰减越大”那肯定是符号写反了振荡只能减弱不能增强。往下走的路径我一般建议两条。一是把平面运动模型扩展到 3DOF加入横摇并接上舵鳍控制然后在 Simulink 里搭闭环船舶动力学模块用 S-Function 包住这个 ode 函数控制律用 PID波浪力 From Workspace 接入第 4 章的谱生成器这样就能看减摇鳍控制下的横摇幅度衰减效果。二是换掉初版简化用切片理论或势流边界元求解 RAO替掉那个人为的 0.9 载荷系数把激励力谱更新成频域传递函数仿真精度可以从 ±25% 收窄到 ±5% 以内这时候就能扛住船级社的耐波性校核了。做船舶 matlab 仿真这几年我最大的体会是先别急着追求六自由度大而全的模型把那两三个自由度的物理机制、量纲、频率关系磨清楚每个结果都能落到理论或实验数据上去核对。数值仿真赚不到“黑匣子”的便宜它所有输出都是方程和参数的镜像你的判断力决定它的价值。遇到诡异结果先怀疑单位、后怀疑符号、再怀疑步长这顺序我踩了无数个坑才总结出来。希望帮到你。本文还有配套的精品资源点击获取