简介：本资源是一份面向船舶运动建模初学者的Matlab基础教学实践包，适用于本科及硕士阶段的船舶与海洋工程、控制工程等相关专业师生开展三自由度MMG标准模型仿真学习。资源基于Matlab 2019a开发，完整实现船舶纵荡、横荡、首摇三自由度非线性运动方程求解与可视化，涵盖主控脚本、MMG参数计算模块、实例工况仿真及运行结果图示，助读者理解MMG模型结构、系数选取逻辑与数值求解流程。压缩包共7个文件（6个.m函数/脚本文件承担建模、求解与绘图功能，1张.jpg直观展示仿真输出曲线），总大小仅31KB，轻量易部署，适合课堂演示、课程设计与自主复现。目前已有2870人学习下载，内容聚焦核心算法实现，无冗余依赖，附带清晰调用关系与注释，可直接运行并支持参数修改拓展。
1. 为什么用 Matlab 模拟船舶三自由度 MMG 模型？——不是为了炫技，而是让操纵性预报真正“可算、可调、可验证”
你手头有一份新船型的线型图和主尺度参数，船东催着要给出“满舵回转直径”和“Z 形操纵响应曲线”，但还没到水池试验阶段，CFD 又太贵、周期太长。这时候，一个能跑通、能调参、能复现经典 MMG 标准试验工况（如纯横荡、纯首摇、耦合回转）的三自由度数学模型，就是你桌上最硬的“数字舵轮”。它不替代实船或水池，但能把设计早期的操纵性风险从“凭经验猜”变成“用方程推”——横荡（surge）、纵荡（sway）、首摇（yaw）三个自由度联动，把非线性水动力导数、螺旋桨推力、舵力、风浪干扰全塞进一组常微分方程里，再用 Matlab 的 ode45 或 Simscape 跑起来。这不是玩具模型：IMO《船舶操纵性标准》附录里的推荐测试流程（如 10°/10° Z 形试验），MMG（Manoeuvring Modelling Group）发布的标准导数表（如 S1720、S1721 系列），以及国内《船舶操纵性预报指南》（CB/T 3903-2021）都明确要求用此类模型做预评估。新手靠它理解水动力耦合逻辑，老手拿它快速迭代舵角/舵面积/螺旋桨布置方案。关键在于：它必须是“可解释的黑匣子”——每个导数项有物理意义，每条曲线能对标实测数据，每次参数修改都能看到响应变化。否则，再漂亮的动画也只是幻灯片。
2. 从 MMG 标准方程到 Matlab 实现：三自由度动力学建模的四步落地法
MMG 三自由度模型不是一堆符号堆砌，而是一套经过数十年水池验证、被全球主流船型数据库（如 KRISO、MARIN、CSSRC）反复校准的工程化方程组。它的核心是把船舶运动分解为三个平动/转动自由度，并对每个自由度建立牛顿-欧拉方程。Matlab 的优势在于：符号计算（Symbolic Math Toolbox）能帮你推导雅可比矩阵，ode 求解器能稳定处理刚性项，Simulink 提供可视化调试界面，而最关键的——你能把 MMG 推荐的无量纲导数（如 $X_{\dot{u}}$、$Y_{v}$、$N_{r}$）直接映射成可调参数，而不是写死在 C 代码里。下面这四步，是我带新人上手时必走的最小闭环路径。
2.1 搭建状态空间与运动学关系：先定义“船在哪、往哪去”
船舶三自由度的状态向量通常定义为：
$$ \mathbf{x} = [u,\ v,\ r,\ x,\ y,\ \psi]^T $$
其中 $u$（纵向速度）、$v$（横向速度）、$r$（首摇角速度）是动力学变量；$x$、$y$、$\psi$（航向角）是运动学变量。注意：MMG 标准中 $u$ 是船体坐标系下沿艏艉轴的速度，不是地速；$v$ 和 $r$ 同理。运动学方程由坐标系转换决定：
% 运动学微分方程（地心坐标系 -> 船体坐标系） dxdt = u*cos(psi) - v*sin(psi); dydt = u*sin(psi) + v*cos(psi); dpsidt = r;提示：这里
cos(psi)和sin(psi)的符号极易写反——记住口诀：“x 方向投影是 u 在 x 轴分量减 v 在 x 轴分量”，画个坐标系草图比背公式管用。很多翻车案例都是因为 psi 定义方向（顺时针 vs 逆时针）和坐标系（北东地 vs 船体前右上）没对齐。
2.2 填入 MMG 标准动力学方程：把“水怎么推船”翻译成代码
MMG 三自由度动力学方程（以 S1720 模型为例）本质是三个力/力矩平衡式：
$$ \begin{cases} m(\dot{u} - vr) = X_H + X_P + X_R \ m(\dot{v} + ur) = Y_H + Y_P + Y_R \ I_z \dot{r} = N_H + N_P + N_R \end{cases} $$
其中 $X_H/Y_H/N_H$ 是船体水动力（含非线性项），$X_P/Y_P/N_P$ 是螺旋桨推力（含伴流、推力系数），$X_R/Y_R/N_R$ 是舵力（含舵角、舵效因子）。Matlab 实现时，绝不能把整个方程当黑箱函数扔给 ode45。我习惯拆成三个子函数，每个对应一类力源，便于单独调试：
function [XH, YH, NH] = hull_forces(u, v, r, delta, L, B, d, beta, t, w, Xuu, Yvv, Yvr, Nvr, Nrr, ...) % 输入：当前速度、舵角、船型参数（L=船长,B=船宽,d=吃水）、MMG 导数 % 输出：船体水动力 XH, YH, NH（单位：N, N·m） % 1. 计算无量纲速度（MMG 标准要求） U = sqrt(u^2 + v^2); beta = atan2(-v, u); % 漂角，注意负号！MMG 定义漂角为正时船向右漂 % 2. 船体水动力（简化版，仅含主要非线性项） XH = m * Xuu * u * abs(u); % 纵向阻力平方项 YH = m * (Yvv * v * abs(v) + Yvr * v * r + Yrr * r * abs(r)); NH = Iz * (Nvr * v * r + Nrr * r * abs(r)); % 3. 加入漂角修正项（如 Ybeta, Nbeta）——此处省略具体系数，实际需查 MMG 表 end参数说明：
Xuu、Yvv等是 MMG 推荐的无量纲导数（如 $X'{uu} = X{uu}/(0.5\rho L^2 d)$），必须乘以 $0.5\rho L^2 d$ 才得有量纲力。ρ 取 1025 kg/m³（海水），L、B、d 单位必须统一为米。新手最大误区：直接把文献里的 $X'{uu}$ 当成 $X{uu}$ 用，导致力大 100 倍——船一启动就飞出屏幕。
2.3 集成螺旋桨与舵力模型：让“人控”介入系统
螺旋桨推力 $X_P$ 和舵力 $Y_R$、$N_R$ 是操纵性的主动控制输入。MMG 对这两者有成熟建模规范：
- 螺旋桨推力：$X_P = (1-t) \rho n^2 D^4 K_T$，其中 $K_T$ 是推力系数，查图谱或用经验公式（如 Keller 公式）；$t$ 是推力减额分数，典型值 0.1~0.25。
- 舵力：$Y_R = \frac{1}{2}\rho V_R^2 A_R a'_H \sin(\delta)$，$N_R = Y_R \cdot x_H$，其中 $V_R$ 是舵处来流速度（含伴流），$A_R$ 是舵面积，$a'_H$ 是舵效因子（MMG 推荐 2.8~3.2），$x_H$ 是舵压力中心距重心距离。
Matlab 中，我把舵角delta和螺旋桨转速n设为外部输入（可随时间变化），这样就能模拟 Z 形试验（δ=±10°阶跃）或回转试验（δ=35°恒定）：
% 在主 ODE 函数中调用 [XP, NP] = propeller_force(u, n, t, KT_func, rho, D); % KT_func 可是查表函数或多项式 [YR, NR] = rudder_force(u, v, r, delta, AR, aH, xH, rho, beta, w); % w 是伴流分数注意：
delta单位必须是弧度！Matlab 三角函数全用 rad。若你从 Excel 导入的是度数，务必delta_rad = deg2rad(delta_deg)。曾见同事因这一步没转，舵力始终为零，调试三天才发现 log 里sin(35)算的是 35 弧度——约等于 sin(2000°)，结果接近 0。
2.4 构建 ODE 求解主干：用 ode45 跑通第一个完整轨迹
把上述所有力汇总，写出完整的六维状态微分方程：
function dxdt = ship_3dof_ode(t, x, params) % x = [u; v; r; x; y; psi] u = x(1); v = x(2); r = x(3); x_pos = x(4); y_pos = x(5); psi = x(6); % 获取输入（此处设为常量，实际可用 interp1 查时间序列） delta = params.delta_func(t); % 舵角函数，如 @(t) (t>10)*pi/18; % 10s 后打 10° 舵 n = params.n_func(t); % 螺旋桨转速 % 计算各力 [XH, YH, NH] = hull_forces(u, v, r, delta, params.L, params.B, params.d, ...); [XP, NP] = propeller_force(u, n, params.t, params.KT, params.rho, params.D); [YR, NR] = rudder_force(u, v, r, delta, params.AR, params.aH, params.xH, params.rho); % 总力 X_total = XH + XP; Y_total = YH + YR; N_total = NH + NP + NR; % 动力学方程（注意质量项） m = params.m; % 船舶排水量（kg） Iz = params.Iz; % 惯性矩（kg·m²） du = (X_total)/m + v*r; % 注意：MMG 方程中 du/dt 含 +vr 项（科氏力） dv = (Y_total)/m - u*r; dr = N_total/Iz; % 运动学方程 dx = u*cos(psi) - v*sin(psi); dy = u*sin(psi) + v*cos(psi); dpsi = r; dxdt = [du; dv; dr; dx; dy; dpsi]; end然后调用求解器：
% 初始条件：静止，艏向 0° x0 = [0; 0; 0; 0; 0; 0]; tspan = [0 120]; % 仿真 120 秒 options = odeset('RelTol',1e-6,'AbsTol',1e-8); % 高精度，尤其对 r 的小变化敏感 [t, x] = ode45(@(t,x) ship_3dof_ode(t,x,params), tspan, x0, options);关键参数
params必须包含所有船型与导数：L,B,d,m,Iz,rho,D,AR,t,w,aH,xH, 以及全部 MMG 导数（Xuu,Yvv,Yvr,Nvr,Nrr等）。我习惯用结构体传参，避免函数签名过长。切记：m是船舶质量（ρ×排水体积），不是排水量吨数！10000 吨排水量 ≈ 1.025e7 kg。
3. MMG 导数从哪来？——三种来源的实操对比与参数校准铁律
MMG 模型的精度，90% 取决于导数的可靠性。没有准确导数，再完美的代码也只是空中楼阁。现实中，导数来源只有三条路：查标准表、CFD 计算、水池试验拟合。每条路都有坑，也都有绕不开的硬规则。
3.1 查 MMG 标准导数表：最快但最危险的起点
MMG 发布过多个标准船型（如 S1720、S1721、S1722）的导数表，覆盖货船、油船、集装箱船。这些表是水池试验的结晶，但绝不是“拿来即用”的万能钥匙。例如 S1720（5万吨散货船）的 $Y_v$ 是 -0.92，而某 20 万吨 VLCC 可能是 -1.35——差 45%，回转直径误差超 30%。我的做法是：
- 第一步：选最接近母型。用船长/吃水比（L/d）、方形系数（Cb）、棱形系数（Cp）三项指标匹配。比如你的船 Cb=0.82，S1720 是 0.81，S1721 是 0.76，则优先取 S1720。
- 第二步：按尺度缩放。MMG 导数是无量纲的，但需还原为有量纲力。缩放公式为： $$ X_{uu}^{real} = X'_{uu} \times \frac{1}{2}\rho L^2 d $$ 注意：$L$、$d$ 必须用你船的实际值，不能用母型值。
- 第三步：重点校验漂角相关项。$Y_\beta$、$N_\beta$ 对 Z 形试验响应影响极大，标准表常偏保守。我一般把 $Y_\beta$ 初值设为标准值 ×1.2，后续用实测数据反调。
提示：MMG 官网（mmg.org）已归档旧表，但最新版（如 2023 年更新的 S1723）需订阅。国内高校图书馆常购有《MMG Standard Test Cases and Derivatives》纸质版，第 47 页起有完整表格。别信网上流传的 PDF——我见过三个版本，$N_{rr}$ 数值相差 0.15，足够让首摇发散。
3.2 CFD 计算导数：可控但昂贵，必须守住三个收敛底线
用 Star-CCM+ 或 OpenFOAM 做强迫振荡（Forced Oscillation）或斜拖（Oblique Towing）计算，是获取定制导数的主流方法。但 CFD 不是点鼠标就出结果的黑盒：
- 网格必须满足 y+ < 1（壁面第一层网格），否则粘性底层解析失败，$Y_v$、$N_r$ 严重偏低；
- 时间步长需满足 CFL < 0.5，尤其对首摇运动，否则 $N_{rr}$ 振荡失真；
- 湍流模型必须用 SST k-ω，RANS 模型中它对分离流和舵后涡最稳健；k-ε 模型在大漂角下会高估阻尼。
我曾用 CFD 算一艘 3000TEU 集装箱船，初始设置用 k-ε，得到 $Y_v = -0.78$，实船 Z 形试验拟合后发现应为 -1.02。换 SST k-ω 重算，结果变为 -0.99，误差 <3%。CFD 报告里写的“网格独立性验证”不是摆设——必须做至少 3 套网格（粗/中/细），看 $Y_v$ 变化 <2% 才算过关。
3.3 水池试验数据反演：最准但最慢，教你用最小二乘法“扒”导数
当你有实船或船模的 Z 形试验数据（δ(t) 和 r(t)、v(t) 曲线），就可以用最小二乘法反演关键导数。这是我在船厂现场最常用的救急手段。核心思想：把动力学方程写成线性组合形式，如：
$$ \dot{r} = a_1 v + a_2 r + a_3 v|r| + a_4 r|r| + a_5 \delta + a_6 \delta|v| $$
其中 $a_i$ 是待求系数，对应 $N_v/m$, $N_r/I_z$ 等。用试验数据构造矩阵 $A$ 和向量 $b$，解 $A \mathbf{a} = b$ 即可。
% 假设已有试验数据：time_vec, r_vec, v_vec, delta_vec % 构造设计矩阵 A（每行一个时刻） A = [v_vec, r_vec, v_vec.*abs(r_vec), r_vec.*abs(r_vec), delta_vec, delta_vec.*abs(v_vec)]; b = diff(r_vec)./diff(time_vec); % dr/dt，用梯形法微分更稳 a = A \ b; % 最小二乘解 % a(1) 就是 N_v / Iz，乘以 Iz 得 N_v Nv_est = a(1) * params.Iz;注意：微分噪声极大！
diff直接用会放大测量误差。血泪经验：先用sgolayfilt(v_vec, 3, 21)（Savitzky-Golay 滤波，窗长 21 点）平滑 v、r、δ，再微分。否则反演结果全是毛刺，根本没法用。
4. 常见问题排查：三自由度 MMG 模型跑不稳、不对、不收敛的 5 个真实翻车现场
模型跑出来轨迹乱飞、速度爆炸、角度突变？别急着重写代码——90% 的问题出在参数、坐标系或数值设置上。以下是我在项目中记录的 5 个高频翻车点，每个都附带现场日志特征和秒级定位法。
4.1 现象：u或v在 1 秒内飙升至 1000 m/s（远超船速）
- 原因：质量
m或惯性矩Iz设为 0 或极小值（如 1），导致加速度 $\dot{u} = X/m$ 爆炸。常见于忘记赋值params.m = 1.025e7;，Matlab 默认m = []，除零报错被忽略后返回Inf。 - 解决：在 ODE 函数开头加断言：
assert(isfinite(params.m) && params.m > 1e3, 'Mass m is invalid!');
4.2 现象：psi（航向角）持续线性增长，不收敛也不振荡
- 原因：运动学方程
dpsi = r正确，但r（首摇角速度）本身为常数非零值。根源是N_total计算错误——最常见是舵力NR = YR * xH中xH符号反了（压力中心在重心前还是后？），导致净力矩恒正。 - 解决：单独打印
N_total随时间变化。若为平稳非零值，立刻检查xH符号（标准定义：重心为原点，向艏为正，舵通常在重心后，xH应为负值）。
4.3 现象：Z 形试验中，r响应滞后严重，峰值比实测晚 5 秒以上
- 原因：
Y_v、N_v导数过小，或Y_{vr}、N_{vr}缺失。MMG 模型中，v*r交叉项对首摇启动至关重要。若只填了Y_v、N_r，漏掉Y_{vr}，则初期响应迟钝。 - 解决：查 MMG 表，确认是否填全所有
*v*r项。S1720 表中Y_{vr}= 0.12，N_{vr}= -0.18，缺一不可。
4.4 现象：ode45报错Failure at t=xxx. Unable to meet integration tolerances，且t值很小（<0.1s）
- 原因：方程存在刚性（stiffness），如
X_{uu}项（$u|u|$）在u≈0附近导数无穷大，导致 ode45 步长崩坏。 - 解决：换刚性求解器
ode15s，或对非线性项做平滑处理：% 替代 u*abs(u)，避免 u=0 处不可导 u_smooth = u * sqrt(u^2 + 1e-6);
4.5 现象：同一舵角下，不同初速（u0=2vsu0=5）的回转直径几乎相同
- 原因：导数未按速度无量纲化。MMG 导数是基于无量纲速度 $u/U_0$ 定义的，但代码中直接用了
u。正确做法是：所有导数项乘以参考速度 $U_0$（如服务航速），或在力计算中显式引入 $U_0$。 - 解决：检查
hull_forces函数，确认XH = m * Xuu * (u/U0) * abs(u/U0) * U0^2类似结构。若没U0，模型就失去速度相似性。
5. 进阶技巧：用 Simulink 实现交互式操纵性试验台，让船东当场改参数看效果
写完脚本只是开始。真正让模型产生工程价值，是把它变成一个“所见即所得”的试验台——船东坐在你电脑前，调个舵角滑块，3 秒后屏幕上就跑出回转轨迹；设计师拖动螺旋桨转速条，Z 形试验的超调量实时刷新。这靠纯 M 文件很难实现，但 Simulink + MATLAB Function Block + Dashboard 控件，5 分钟就能搭出来。
5.1 构建 Simulink 主框架：信号流比代码更直观
核心是把六维 ODE 封装成一个自定义模块。我习惯用MATLAB Function Block（而非 S-Function），因为无需编译，支持所有 MATLAB 语法：
- 新建 Simulink 模型，拖入
MATLAB Function模块； - 双击进入编辑器，粘贴
ship_3dof_ode的核心计算逻辑（去掉ode45调用，只留dxdt计算）； - 设置输入端口：
u_in,v_in,r_in,delta_in,n_in（来自 Dashboard 滑块）； - 设置输出端口：
u_out,v_out,r_out,x_out,y_out,psi_out； - 在模型中添加
Integrator模块链，构成状态反馈闭环（x_dot→Integrator→x→MATLAB Function输入）。
关键细节：
Integrator模块的初始条件必须设为[0;0;0;0;0;0]，且勾选 “Enable zero-crossing detection” ——这对捕捉舵角阶跃时刻的瞬态响应至关重要。
5.2 添加 Dashboard 控件：让参数调整像调音台一样简单
Simulink Dashboard 提供滑块（Slider）、旋钮（Knob）、开关（Toggle Switch）等控件，直接绑定模型参数：
- 拖入
Slider，右键 →Block Parameters→Signal Name设为delta_cmd； - 在
MATLAB Function内部，用coder.extrinsic('evalin')读取工作区变量，或更稳妥地：将delta_cmd作为额外输入端口接入； - 同样添加
Knob控制n_cmd（螺旋桨转速），Edit Box输入Xuu等导数，实时生效。
% 在 MATLAB Function Block 中 function [u_out, v_out, r_out, x_out, y_out, psi_out] = fcn(u_in, v_in, r_in, ... delta_in, n_in, Xuu, Yvv, Yvr, Nvr, Nrr, m, Iz, L, B, d, rho, D, AR, aH, xH, t, w) % 所有参数均通过输入端口传入，无需全局变量 % ... 计算逻辑同前 ...5.3 配置 Scope 与 XY Graph：一目了然看操纵性
Scope显示r(t)、v(t)曲线，用于分析 Z 形试验的响应时间、超调量；XY Graph绘制(x_out, y_out)，实时画出回转轨迹，标出直径（max(y)-min(y)）；- 添加
Dashboard Gauge显示当前psi（航向角），Dashboard Display显示u（船速）。
实战技巧：用
Simulation Data Inspector（运行后点击图标）对比多组参数下的r(t)曲线。选中两条线，右键 →Compare，它会自动计算 ITAE（Integral of Time-weighted Absolute Error），量化响应差异——这比肉眼判断“哪个更快”专业十倍。
5.4 导出为 Standalone App：告别 Matlab License 依赖
客户没装 Matlab？用Application Compiler打包成独立 exe：
APPS→Application Compiler；- 主文件选你的 Simulink 模型
.slx； - 勾选
Include MATLAB Runtime（自动下载）； - 编译后生成
setup.exe，客户双击安装即可运行，无需 Matlab。
我的血泪教训：第一次打包时忘了勾选
Simulink依赖项，exe 启动报错 “Simulink not found”。后来发现必须在Library Compiler中手动添加simulink和simscape两个 toolbox。现在我的编译清单固定为：matlab,simulink,simscape,control,optimization。
最后说句实在话：这个模型不会让你一夜之间成为水动力专家，但它会把你从“等试验报告”的被动角色，变成“主动推演方案”的技术枢纽。我见过太多项目，因为没在设计早期跑通这个三自由度模型，导致实船回转半径超标，后期只能加装侧推——成本增加百万。而用这套流程，你在 PPT 里展示的每一条曲线，背后都是可追溯、可修改、可验证的方程。希望帮到你。
本文还有配套的精品资源，点击获取
