Matlab模拟船舶操纵运动模型：从运动方程到旋回试验

简介面向本科、硕士及船舶运动控制初学者基于Matlab2019a编写的船舶操纵运动仿真模型资源包。主要包含船舶平面运动建模、回转试验与Z形试验仿真脚本涉及舵角控制策略、相对距离与偏航位置调整等核心模块适合用于课程设计、毕业设计或科研预研。资源共16个文件其中9个.m脚本覆盖主程序与功能函数配合2个xlsx数据表和1个.mat数据文件可完成参数输入与结果存储另有4张jpg/png图片用于展示仿真曲线与船舶轨迹压缩包仅81KB结构精简轻量便于快速部署。目前已有1380人学习下载对于希望快速掌握船舶操纵运动仿真框架、理解回转半径与舵角关系的教研用户可直接运行并在此基础上扩展自己的控制算法。1. Matlab模拟船舶操纵运动模型从舵令到航迹的核心思路Matlab模拟船舶操纵运动模型简单说就是把一条船当作水平面内的刚体在三自由度运动方程上做数值积分把舵令→航迹这条链路完整跑通。要处理的事情包括定义坐标系与操纵性导数、组带有附加质量和重心偏置的耦合质量矩阵、用ode45积分得到旋回圈和Z形响应再从中提取战术直径、纵距、横距这些操纵性指标。适合做船舶操纵性评估、自动舵或航迹控制算法预研、课程设计与毕业设计的读者。比较反直觉的一点是运动方程本身半小时能写完真正卡住大家的往往是符号约定、参数量级和指标提取方式这三件事比解方程更值得花时间。2. 船舶操纵运动模型的运动方程、坐标系统与参数约定2.1 大地固定系与随船坐标系状态量怎么选模拟船舶操纵运动的第一步是把参考系定死。工程上几乎都采用双坐标系大地固定系 O-X₀Y₀Z₀ 固定在地球上原点取在仿真起点附近随船系 G-xyz 原点在船体重心x 轴指向船艏y 轴指向右舷z 轴向下。操纵运动只研究水平面内的纵荡、横荡和艏摇三个自由度所以状态量取 6 个状态量含义单位x_e, y_e重心在大地系中的位置mψ航向角随船系 x 轴相对大地系 X₀ 轴的夹角radu随船系纵向速度即前进速度m/sv随船系横向速度右舷方向为正m/sr转艏角速度右转为正rad/s这里最容易出问题的就是ψ的符号因为 z 轴向下右手系里ψ顺时针增大也就是说船右转时ψ变大。很多二维绘图场景习惯把 y 轴画成向上一旦把坐标轴方向改掉整个旋回图就会镜像但方程本身没错。建议从头到尾固定x 向右、y 向右手定则的 y、z 向下这一套别在半路切换。随船系速度到大地系位置的转换就是标准的二维旋转dxe u*cos(psi) - v*sin(psi); dye u*sin(psi) v*cos(psi); dpsi r;这组运动学关系没有任何近似是精确的。注意横向速度 v 直接参与位置更新不要因为船主要往前开就把它丢掉——旋回过程中 v 的量级能达到 u 的 8%10%直接影响轨迹形态。2.2 三自由度运动方程与质量矩阵求逆动力学方程写在随船系里最方便因为船体水动力导数都是在随船系下测的。带附加质量和重心偏置的常用形式如下(m mx)·u̇ (m my)·v·r X(m my)·v̇ m·xG·ṙ -(m mx)·u·r Ym·xG·v̇ (Iz Jzz)·ṙ -m·xG·u·r N其中 X、Y、N 是船体、螺旋桨、舵产生的合外力与艏摇力矩mx、my、Jzz 是附加质量与附加惯矩xG 是重心相对原点的纵向位置。两个科氏项 (mmy)·v·r 和 -(mmx)·u·r 是离心力在随船系下的分量旋回时它们和舵力、船体力共同决定稳态转艏角速度不能省略。注意第二个和第三个方程里 v̇ 和 ṙ 通过 m·xG 耦合在一起。只要重心不在随船系原点上就必须做 2×2 矩阵求逆才能解出 v̇ 和 ṙ。在Matlab中定义微分方程时最省事的做法是每次调用都做一次矩阵左除M [P.mP.my, P.m*P.xG; P.m*P.xG, P.IzP.Jz]; rhs [Y - (P.mP.mx)*u*r; N - P.m*P.xG*u*r]; acc M \ rhs; dv acc(1); dr acc(2);这里的 M 是常值矩阵放在子函数里每次左除的开销可以忽略。有不少初学代码会把 xG 设成 0 避开耦合这能跑通但换到真实船的参数时货船 xG 一般在 -3%L 到 -5%L转艏和横荡的耦合会明显改变旋回直径建议一开始就按耦合形式写。2.3 Abkowitz与MMG建模路线怎么取舍有了方程框架下一步是决定 X、Y、N 怎么表达常见两条路线。Abkowitz 模型把船体、桨、舵的合力统一写成关于 u、v、r、δ 及其交叉项的多项式系数通过平面运动机构试验整体拟合。优点是形式紧凑、系数数量少适合教学和整体参数辨识缺点是物理分界模糊换舵或者换桨之后整套系数都要重拟合。MMG 模型Maneuvering Modeling Group把力拆成船体力 XH、螺旋桨力 XP、舵力 XR 三个模块分别建模船体力常用漂角 β 和无量纲转艏角速度的多项式表示。优点是模块独立换螺旋桨、加浅水修正、加舵效修正都只动局部缺点是系数多入门门槛高。船模试验和实船操纵性预报的工程实践中MMG 体系更常见航海模拟器里的六自由度模型也基本是它的扩展。对比项Abkowitz 多项式MMG 分模块力的组织整体多项式船体/桨/舵分离系数数量少多换部件后需重拟合只换对应模块适合场景教学、整体辨识工程预报、模拟器实际做 Matlab 仿真时多数人的做法是混合写船体力用 Abkowitz 截断多项式舵力和螺旋桨推力单独列项。这样做既能控制参数规模又保留了操纵部件独立调整的能力下面章节的代码就采用这种写法。3. 在Matlab中搭建船舶操纵运动模型的可运行代码3.1 用结构体管理船型参数与操纵性导数仿真代码的参数组织直接决定后面调参是否痛苦。建议把所有船型参数和水动力导数放进一个结构体 P用一个独立 m 函数返回不要散落在脚本里。下面这份参数表是教学示例值量级参考一艘 160 m 级货船的公开操纵性数据做真船评估时需要用船模试验或 CFD 结果替换参数含义单位示例值L / U0垂线间长 / 设计航速m, m/s160.9, 7.7m / mx / my质量 / 纵荡附加质量 / 横荡附加质量kg1.704e7 / 8.52e5 / 1.36e7Iz / Jzz艏摇惯矩 / 附加惯矩kg·m²2.70e10 / 1.08e10xG重心纵向位置舯前为正m-5.96Yv, Yv2横荡线性/平方阻尼导数N/(m/s), N/(m/s)²-4.0e6, -3.0e5Yr转艏引起的横荡力导数N/(rad/s)2.2e8Nv, Nv2横荡引起的艏摇力矩导数N·m/(m/s), N·m/(m/s)²-5.0e7, -5.0e7Nr艏摇阻尼导数N·m/(rad/s)-5.0e9Yδ, Nδ舵力/舵力矩导数右舵为正N/rad, N·m/rad-4.0e6, 1.8e8Xu, Xu2纵荡线性和平方阻尼N/(m/s), N/(m/s)²-1.5e6, -2.0e5Xvv, Xrr漂角和转艏引起的阻力N/(m/s)², N/(rad/s)²-4.0e6, -2.0e8T0简化推力模型系数kg/m1.5e5对应代码function P ship_params() % 船型与操纵性参数SI单位。示例值量级参考约160m货船 % 用于跑通仿真流程真船评估时用船模试验或CFD结果替换。 P.L 160.9; % 垂线间长 Lpp, m P.rho 1025; % 海水密度, kg/m^3 P.m 1.704e7; % 排水质量, kg P.Iz 2.70e10; % 绕z轴转动惯量, kg*m^2 P.xG -5.96; % 重心纵向位置(舯前为正), m P.U0 7.7; % 设计航速, m/s P.mx 0.05*P.m; % 纵荡附加质量 P.my 0.80*P.m; % 横荡附加质量 P.Jz 0.40*P.Iz; % 艏摇附加惯矩 % 船体水动力导数Abkowitz截断形式 P.Xu -1.5e6; % N/(m/s) P.Xu2 -2.0e5; % N/(m/s)^2 P.Xvv -4.0e6; % N/(m/s)^2 P.Xrr -2.0e8; % N/(rad/s)^2 P.Yv -4.0e6; % N/(m/s) P.Yv2 -3.0e5; % N/(m/s)^2 P.Yr 2.2e8; % N/(rad/s) P.Nv -5.0e7; % N*m/(m/s) P.Nv2 -5.0e7; % N*m/(m/s)^2 P.Nr -5.0e9; % N*m/(rad/s) % 舵力导数约定 delta0 为右舵 P.Yd -4.0e6; % N/rad P.Nd 1.8e8; % N*m/rad P.Xd2 -1.0e6; % N/rad^2 舵的阻力分量 % 螺旋桨推力简化速度修正模型 Xp T0*(U0^2 - u^2) P.T0 1.5e5; % kg/m end所有导数都带符号这是有意的Yv、Nr 是阻尼项必须为负Nv 对航向稳定的船通常为负Yδ 与 Nδ 的符号由舵角定义决定。换用论文里的无量纲系数时第一步就是把符号约定对齐否则仿真跑出来船会朝反方向转。3.2 运动微分方程函数质量矩阵左除的做法在Matlab中定义微分方程推荐把状态导数写成独立 m 函数。状态向量 x [xe; ye; psi; u; v; r]函数返回 dxrudder_fun 是舵角函数句柄这样同一份方程代码可以服务直航、阶跃、蛇行各种工况function dx ship_3dof(t, x, P, rudder_fun) % 三自由度船舶操纵运动方程 psi x(3); u x(4); v x(5); r x(6); delta rudder_fun(t); % 舵角, 右舵为正, rad % 船体力Abkowitz多项式的截断形式 du u - P.U0; XH P.Xu*du P.Xu2*du*abs(du) P.Xvv*v^2 P.Xrr*r^2; YH P.Yv*v P.Yv2*v*abs(v) P.Yr*r; NH P.Nv*v P.Nv2*v*abs(v) P.Nr*r; % 舵力与螺旋桨推力 XR P.Xd2*delta^2; YR P.Yd*delta; NR P.Nd*delta; XP P.T0*(P.U0^2 - u^2); X XH XR XP; Y YH YR; N NH NR; % 质量矩阵求逆解出加速度 M [P.mP.mx 0 0; 0 P.mP.my P.m*P.xG; 0 P.m*P.xG P.IzP.Jz]; rhs [X (P.mP.my)*v*r; Y - (P.mP.mx)*u*r; N - P.m*P.xG*u*r]; acc M \ rhs; dx [u*cos(psi) - v*sin(psi); u*sin(psi) v*cos(psi); r; acc(1); acc(2); acc(3)]; end三点说明。第一船体力里的 v²、v|v| 写法是有讲究的v|v| 保证力的方向始终与速度方向相反v² 在某些文献里也常见但遇到大漂角时 v|v| 更符合物理两种写法在 v 为正时等价。第二舵角 delta 在函数开头一次性取出后续所有舵力项共用不要在每个表达式中重复调用 rudder_fun否则遇到查表型舵令会多出大量插值开销。第三X 中加了 (mmy)·v·r 科氏项而 Y、N 行里对应的是 -(mmx)·u·r 和 -m·xG·u·r这三个符号是最容易抄错的地方建议每次改完先做一次直航仿真确认 u 能保持 U0。3.3 用ode45积分并正确处理角度回绕积分器直接用 ode45。船舶操纵运动是慢变的平滑系统ode45 的默认容差通常够用但长时积分建议收紧P ship_params(); x0 [0; 0; 0; P.U0; 0; 0]; opts odeset(RelTol, 1e-6, AbsTol, 1e-7); % 第1段直航稳定5秒让速度场收敛 [t1, x1] ode45((t,x) ship_3dof(t,x,P,(t) 0), [0 5], x0, opts); % 第2段阶跃右满舵35°积分600秒 [t2, x2] ode45((t,x) ship_3dof(t,x,P,(t) deg2rad(35)), [5 600], x1(end,:), opts); t [t1; t2]; x [x1; x2];关键在中间那句 x1(end,:)。直接让 ode45 从 t0 就施加满舵船体会在第一个积分步里同时经历纵向速度调整和横向响应旋回圈起点的拉出过程会和真实试验不一致。真实操船试验是先在直航工况下稳定再执行舵令所以仿真也要分两段把第一段的末状态作为第二段初值。这段代码里还有一个隐蔽点微分方程中的 ψ 始终保持连续增长不做取模运算。如果有人在方程里写 psi wrapToPi(psi)ode45 会在 2π 跳变点处检测到不连续步长被压到极小仿真时间暴增事件检测也会出现假触发。正确做法是积分过程保持 ψ 连续只在画图或提取指标时用 wrapToPi 处理显示值。3.4 舵令输入的三种组织方式rudder_fun 用函数句柄的好处是可以随时换舵令形式。最常见三种。常数舵角给 (t) deg2rad(10)阶跃舵令给 (t) deg2rad(35)*(t5)这在前面的旋回代码里已经用了。第三种是时间序列舵令比如把试验记录的舵角 csv 导入Matlab后做插值tbl readmatrix(rudder_series.csv); % 第1列时间s, 第2列舵角deg rudder_fun (t) deg2rad(interp1(tbl(:,1), tbl(:,2), t, previous, extrap));插值方法建议选 previous因为真实舵机是保持舵角直到下一个指令到达线性插值会把舵角变化过程拖成斜坡和舵机特性不符。如果后面要把这套模型放进 Simulink 做闭环控制同样可以把 ship_3dof 的核心部分封装成 MATLAB Function 模块输入舵角、输出状态量脚本阶段调试通过的参数可以直接搬过去。4. 船舶操纵运动模型的旋回试验与Z形试验仿真4.1 满舵旋回试验纵距、横距、战术直径的提取跑完第 3 章的旋回仿真最直接的动作是画轨迹并提取标准指标。画轨迹用 plot 加 axis equal想演示动态过程可以换成 cometfigure; plot(x(:,1)/1000, x(:,2)/1000, LineWidth, 1.2); axis equal; grid on; xlabel(x_e (km)); ylabel(y_e (km)); title(右满舵35°旋回试验轨迹);指标提取要回到真实定义不能直接找离起点最远的点。以舵令开始执行的位置为基准把轨迹投影到初始航向线上。指标定义提取方式纵距 Advance航向转90°时沿初始航向的位移find(dpsi 90°) 后取 x_adv横距 Transfer航向转90°时垂直初始航向的位移同点取 y_tr战术直径航向转180°时垂直初始航向的位移find(dpsi 180°) 后取定常回转直径稳定回转段轨迹直径取最后一段曲率半径base x1(end,:); % 舵令执行瞬间的状态 psi0 base(3); dpsi x2(:,3) - psi0; i90 find(dpsi deg2rad(90), 1, first); i180 find(dpsi deg2rad(180), 1, first); dx x2(:,1) - base(1); dy x2(:,2) - base(2); x_adv dx*cos(psi0) dy*sin(psi0); % 投影到初始航向 y_tr -dx*sin(psi0) dy*cos(psi0); fprintf(纵距%.1f m (%.2fL)\n, x_adv(i90), x_adv(i90)/P.L); fprintf(横距%.1f m (%.2fL)\n, y_tr(i90), y_tr(i90)/P.L); fprintf(战术直径%.1f m (%.2fL)\n, abs(y_tr(i180)), abs(y_tr(i180))/P.L);这里用投影而不是直接用 x、y 坐标是因为如果初始航向不是 0直接读坐标会混入航向偏差。顺带看速度曲线旋回中 u 会下降这套参数大概掉 6%~8%真实货船满舵旋回一般掉 15%~25%误差主要来自螺旋桨推力模型过简只用了速度平方修正而没有用敞水特性曲线不影响方法演示。4.2 用Events事件函数做10°/10°Z形试验Z形试验和旋回试验的区别在于舵令由航向反馈触发先右舵 10°航向达到 10° 立即反向左舵 10°航向达到 -10° 再反向如此交替。Matlab 里实现这类状态到达阈值就切换的标准做法是 ode45 的 Events 功能。function [value, isterminal, direction] event_psi(t, x, thr, dir) value x(3) - thr; % 航向偏差过零即事件 isterminal 1; % 触发后停止积分 direction dir; % 只检测指定方向的穿越 enddirection 参数必须设为和当前舵角同号右舵时 ψ 上升只检测正向穿越反向操舵后 ψ 下降只检测负向穿越。如果设成 0 双向检测反向瞬间 ψ 恰好贴着阈值数值噪声可能造成连续触发。function [t_all, x_all, rud_log, ev_t] run_zigzag(P, t_final) x0 [0; 0; 0; P.U0; 0; 0]; delta deg2rad(10); % 先右舵10° thr deg2rad(10); % 目标航向偏差10° t0 0; t_segs {}; x_segs {}; rud_log []; ev_t []; for k 1:30 opts odeset(Events, (t,x) event_psi(t,x,thr,sign(delta)), ... RelTol, 1e-6, AbsTol, 1e-7); [t_seg, x_seg, te, xe] ... ode45((t,x) ship_3dof(t,x,P,(t) delta), [t0 t_final], x0, opts); t_segs{end1} t_seg; x_segs{end1} x_seg; rud_log [rud_log; delta*ones(numel(t_seg),1)]; if isempty(te) % 没有再触发事件就结束 break; end ev_t(end1) te(1); t0 te(1); x0 xe(1,:).; delta -delta; % 反向操舵 thr -thr; % 阈值反向 end t_all vertcat(t_segs{:}); x_all vertcat(x_segs{:}); end逐段积分的写法比在 ODE 函数内部用 persistent 变量存状态安全得多。persistent 变量在 ode45 的试探步里会被反复改写生成的舵令时间轴是错的而每段积分之间显式传递 x0事件时刻 te 就是精确的操舵换向时刻。画图时把航向和舵角用双 y 轴画在一起就能直接读出超越角。4.3 参数敏感性从特征值判稳到先动哪个参数换参数之前先做稳定性自检。把方程在直航状态下线性化得到横荡-艏摇二阶系统稳定性取决于 M₂⁻¹A 的特征值实部是否全为负M2 [P.mP.my, P.m*P.xG; P.m*P.xG, P.IzP.Jz]; A [P.Yv, P.Yr-(P.mP.mx)*P.U0; P.Nv, P.Nr-P.m*P.xG*P.U0]; lambda eig(M2 \ A); if all(real(lambda) 0) fprintf(方向稳定特征值实部 %s\n, mat2str(real(lambda)., 3)); else fprintf(不稳定检查 Yv、Nr 的负号与 Nv 的量级\n); end与其死记文献里的稳定性判据公式不如用这个特征值检查它不依赖具体符号约定任何参数组代入都能判断。跑通之后调参顺序有讲究。第一优先看 Nv 和 Nr这两个导数决定航向稳定性和超越角大小把 |Nr| 调小 30%Z 形试验的超越角会明显变大旋回直径变化不大。第二再看 Yv它主要影响漂角和横荡响应快慢。最后才动 Yδ 和 Nδ舵力导数直接控制旋回圈大小Nδ 增大 10% 战术直径能缩小约 7%。如果手里有实船或船模试验航迹想反推导数可以用 Optimization Toolbox 的 lsqnonlin以仿真航迹与实测航迹的残差为目标函数做参数辨识目标函数里嵌的就是这一整套仿真代码。5. 船舶操纵运动模型的无量纲化、自检与Simulink迁移5.1 无量纲系数与实船数据的换算论文和船模试验报告里的导数几乎都是无量纲的直接代入 SI 单位的方程会差出好几个数量级。常用约化方式是 SNAME 第一体系力除以 (0.5·ρ·L²·U²)力矩除以 (0.5·ρ·L³·U²)时间用 L/U。换算代码就三行scaleF 0.5*P.rho*P.L^2*P.U0^2; % 力的约化因子 scaleM 0.5*P.rho*P.L^3*P.U0^2; % 力矩约化因子 % 例论文给 Yd 0.0032转成 SI 的 Yd P.Yd 0.0032*scaleF; % 注意核对论文的舵角单位是rad还是deg最容易踩的坑有三个舵角在论文里可能用角度制约化速度取的是试验航速而不是设计航速有些文献用 L²d 而不是 L³ 做力矩约化。任何一组系数用之前先拿稳定性自检的特征值代码过一遍再跑一个 10° 小舵角旋回看方向是否正确。5.2 仿真结果自检的三个检查项每次改完参数按固定顺序检查三个现象。第一零舵角直航 60 秒u 应稳定在 U0 附近v 和 r 应衰减到接近 0任何持续增长的 v 都说明阻尼符号有误。第二右舵 35° 旋回航向 ψ 必须持续增大、轨迹向右偏如果反向检查 Yδ、Nδ 的符号约定。第三Z 形试验的航向振荡中心应回到初始航向附近超越角为正且不过大如果 Z 形试验里航向越过阈值后还在同一方向继续跑很远说明反向舵效不足优先调 Nδ。5.3 向Simulink与试验数据驱动扩展模型验证通过后常见的扩展方向有两个。一是封装成 Simulink 的 MATLAB Function 模块把舵角作为输入、状态量作为输出外面接 PID 航向控制器就能做闭环仿真脚本阶段的参数和符号约定原样保留。二是用试验数据驱动把实船记录的时间-舵角序列通过 readmatrix 导入按第 3.4 节的方式构造成 rudder_fun对比仿真航迹与实测航迹这一步往往是操纵性模型标定工作的起点。本文还有配套的精品资源点击获取