简介船舶操纵性研究者和相关专业学生可使用这套MATLAB仿真代码包对回转试验与Z形操舵试验进行数值建模与仿真。压缩包内共3个文件包含turningtest.m与zigzagtest.m两个主程序脚本以及一份用于对比验证的船模试验数据Excel表格。代码依据姚建喜《船舶操纵性理论基础》中的数值求解流程编写通过仿真结果可直观看到回转圈、Z形操舵响应曲线并与实际试验数据进行对照。资源整体约15KB结构紧凑便于直接运行和修改参数目前已有1971人学习下载。读者可基于该代码理解船舶运动方程离散化求解过程掌握典型操纵性仿真试验的实现方法也可将自带数据替换为自测数据用于课程设计、论文验证或操纵性预报初步评估。1. 船舶回转试验与Z形操舵仿真让MATLAB仿真代码先跑起来压缩包里同时放着 turningtest.m、zigzagtest.m 和一份数据.xlsx这比只给一段孤立代码有价值得多。两个脚本分别仿真回转试验与 Z 形操舵试验Excel 数据是采集端用 getdate 获取时间戳后导出的船模试验记录专门用来和仿真结果逐点对比。数值求解以姚建喜《船舶操纵性理论基础》为蓝本对三自由度运动方程做四阶龙格-库塔积分能够预报回转直径、超越角、操纵周期这些标准操纵性指标。这份资源适合做船舶操纵性课程设计、需要验证水动力导数取值或者想把仿真航迹与船模试验数据对齐做误差分析的从业者。水动力导数标定到位后用数值求解预报船舶操纵性是稳定可行的下面从建模、仿真到数据校核依次拆开讲。2. 三自由度运动方程搭建与MMG分离建模2.1 为什么不用一二阶Nomoto模型Nomoto 一阶模型用 K、T 两个参数就能描述航向对舵的响应在自动舵设计里非常好用但它只输出航向给不出轨迹坐标和横荡速度画不出回转圈也算不了定常回转直径。回转试验和 Z 形试验的仿真都需要轨迹级输出所以这里采用状态空间形式的三自由度方程把纵荡、横荡、艏摇三个自由度耦合在一起求解。标准的三自由度操纵方程写成下面的形式(m mx)·u̇ - (m my)·v·r XH XP XR(m my)·v̇ (m mx)·u·r YH YR(Iz Jz)·ṙ NH NR左侧是惯性力和附连水质量项右侧是船体力、螺旋桨力、舵力。代码按照姚建喜《船舶操纵性理论基础》里的数值求解思路组织力和力矩按 MMG 分离模型分别计算再叠加而不是用一个整体多项式去拟合。这样做的直接好处是想单独调舵力或者修改某个水动力导数时改动局限在一个子函数里不会牵连其他项。2.2 状态向量与参数结构体仿真里的状态有六个纵荡速度 u、横荡速度 v、艏摇角速度 r、艏向角 psi以及大地固定系的坐标 x、y。它们的关系是psi 的导数是 rx 的导数是 u·cos(psi) - v·sin(psi)y 的导数是 u·sin(psi) v·cos(psi)。后面两个坐标变换式决定了回转圈的形状写的时候别把符号搞反。参数建议全部收进一个结构体 p 里函数签名统一写成 f(s, p)调参时只改一处。状态含义初始值更新依据u纵荡速度 m/sp.u0纵向力/质量v横荡速度 m/s0横向力/质量r艏摇角速度 rad/s0艏摇力矩/惯量psi艏向角 rad0rx大地坐标东向 m0u·cos(psi)-v·sin(psi)y大地坐标北向 m0u·sin(psi)v·cos(psi)水动力导数的数量级可以先用下面的示意值起步再用数据.xlsx 里的试验结果校核p.u0 2.0; p.Lpp 3.0; p.m 250; p.Izz 280; p.Xuu -2.5; p.Yv -450; p.Yr -30; p.Nv -80; p.Nr -160; p.mx 0.05*p.m; p.my 0.75*p.m; p.Jz 0.25*p.Izz;这里的 mx、my 是附连水质量Jz 是附加艏摇惯量。my 一般在 0.7~0.9 倍船体质量之间Jz 约在 0.2~0.3 倍转动惯量之间具体取值要和模型试验的无因次化结果对应直接抄别人论文的数值经常对不上。2.3 RK4积分与仿真步长选择操纵性仿真的时间尺度是几十秒到几分钟定步长四阶龙格-库塔是性价比最高的选择稳定性好核心实现只要四行代码。步长取 0.05~0.2 s 都行我一般取 0.1 s步长超过 0.5 s 后舵效和回转圈半径会出现肉眼可见的离散误差。function s1 rk4_step(f, s, dt, p) k1 f(s, p); k2 f(s 0.5*dt*k1, p); k3 f(s 0.5*dt*k2, p); k4 f(s dt*k3, p); s1 s (dt/6) * (k1 2*k2 2*k3 k4); end函数句柄 f 每次根据当前状态返回六个状态的导数向量。注意 f 里要把 u、v、r 的导数从方程左侧解出来再返回也就是先算合力再做 (mmx) 之类的质量除法顺序反了会出现量纲错误现象是速度莫名其妙发散。3. 回转试验仿真turningtest.m的舵令执行与回转圈提取3.1 舵令处理一阶惯性加舵速限制回转试验的标准做法是船以直航速度稳定航行舵迅速打到满舵 35°保持到航向转过 360° 以上。仿真里很多人直接给舵角一个阶跃这会让初始段舵力突变算出来的战术直径偏小。真实的舵机有响应时间我在 turningtest.m 里用一阶惯性环节近似舵机动态时间常数 1.5~2 s同时限制舵速。这样既贴近物理又避免数值突变还能在输出里清楚看到舵令和实际舵角的相位差。3.2 主循环与合力计算turningtest.m 的主循环结构如下每一步先更新舵角再算合力最后做一步 RK4dt 0.1; t_end 180; cmd 35 * pi/180; % 满舵35° T_d 1.5; % 舵机时间常数 s [p.u0; 0; 0; 0; 0; 0]; % [u; v; r; psi; x; y] delta 0; for k 1:round(t_end/dt) delta delta (cmd - delta)/T_d * dt; % 一阶惯性 [X, Y, N] total_force(p, s, delta); % 合力/合力矩 f (s) state_eq(s, X, Y, N, p); s rk4_step(f, s, dt, p); out(k, :) [s; delta]; % 状态舵角一起存 endstate_eq 里先解出三个速度导数再拼上 psi_dot r 和 x_dot、y_dot 的坐标转换式。舵角 delta 存进 out 的最后一列是为了画图时能对照舵令与航向的相位关系。total_force 内部的舵力估算要取船尾处的合成来流速度而不是简单用 u否则大舵角下舵效会被明显高估。舵力这部分我一般直接在 total_force 里嵌下面的形式ua sqrt((u - 0.1*r)^2 v^2); % 舵处来流速度近似 L_r 0.5 * 1025 * p.Ar * ua^2 * p.CLa * delta; % 舵升力 Y Y L_r; N N L_r * p.x_R; % x_R 为舵力作用点纵向坐标来流速度用 u、v、r 的合成而不是单独用 u是因为回转过程中漂角和船尾横向速度会显著改变舵效这也是初版仿真和试验数据对不齐时最常见的误差来源。3.3 回转直径与战术直径的提取仿真结束后画 x-y 轨迹取航向转过 360° 之后的稳定段用三点定外接圆求半径半径的两倍就是定常回转直径。战术直径则是航向改变 180° 时船相对初始航线的横向位移直接取该时刻的 y 坐标绝对值即可不需要额外乘系数。psi_deg out(:,4) * 180/pi; % 累计航向RK4积分不取模 i180 find(psi_deg 180, 1, first); D_tac abs(out(i180, 6)); % 战术直径 i360 find(psi_deg 360, 1, first); seg out(i360:end, 5:6); % 稳定段 x,y P1 seg(1,:); P2 seg(round(end/2),:); P3 seg(end,:); a norm(P2-P3); b norm(P3-P1); c norm(P1-P2); s (abc)/2; R a*b*c / (4*sqrt(s*(s-a)*(s-b)*(s-c))); % 海伦公式定外接圆 D_std 2*R; % 定常回转直径注意 psi_deg 是积分得到的累计值超过 360° 后不要取模否则 i360 永远找不到。海伦公式要求三点不共线稳定段的轨迹近似圆弧三点间距拉开到几十个采样点以上数值上是稳定的。拿到这两个直径后再和试验数据里对应的量做对比。4. Z形操舵试验仿真zigzagtest.m的判舵逻辑与超越角提取4.1 判舵规则与超越角定义Z 形试验用来评价船舶对舵的响应快慢和阻尼特性。以最常用的 10°/10° 试验为例直航稳定后把舵打到右 10°当航向相对起点偏转到 10° 时把舵反向压到左 10°当航向相对新的起点反向偏转 10° 时再回右舵。如此反复得到一条锯齿状的航向曲线。超越角是每次反向操舵后航向继续冲过目标角的余量。比如第一次反向时航向实际冲到 12.3°则第一超越角为 2.3°。超越角小说明船舶阻尼大、惯性小反过来则说明舵效相对惯性偏弱。这个指标对水动力导数 Yv、Nr 特别敏感是校核参数时的首选目标。4.2 航向偏差检测与舵角翻转实现zigzagtest.m 的判舵核心是记录每次操舵翻转时刻的航向作为 psi_start然后判断当前航向相对 psi_start 的偏差是否达到设定试验角psi_cmd 10 * pi/180; % 10°/10° Z形试验 delta psi_cmd; % 先压右舵 psi_start 0; % 本段操舵起点航向 flipped false; % false: 右舵段; true: 左舵段 for k 1:round(t_end/dt) dev psi - psi_start; if ~flipped dev psi_cmd delta -psi_cmd; psi_start psi; flipped true; elseif flipped dev -psi_cmd delta psi_cmd; psi_start psi; flipped false; end [X, Y, N] total_force(p, s, delta); s rk4_step((s) state_eq(s, X, Y, N, p), s, dt, p); end判断里用 dev 而不是直接用 psi 和 0 比较是因为每次翻转后基准航向变了。若写成 psi 超过 ±10° 就翻转第二次翻转会提前发生超越角会被严重低估。psi_start psi 这一行必须放在更新状态之前取当前时刻的航向作为新基准不能用 RK4 更新后的值。4.3 仿真与试验数据的指标对比表跑完 10°/10° Z 形试验后从航向曲线里提取第一超越角、第二超越角、操纵周期和每个转向段的时间。以一组船模参数为例仿真输出与数据.xlsx 中的试验记录对比如下指标仿真值试验值偏差第一超越角 (deg)2.62.30.3第二超越角 (deg)3.12.80.3操纵周期 (s)34.233.80.4平均舵效时滞 (s)1.21.4-0.2表格里这组偏差说明模型整体阻尼略偏小。如果第一超越角偏大优先增大 |Nv| 或 |Yv|如果周期偏长说明惯性项或者附连水质量取大了。注意先看超越角偏差因为它同时反映线性和非线性阻尼的平衡比单独看周期更能定位问题。5. 数据.xlsx试验值与仿真结果的时间对齐与误差校核5.1 试验数据读入与航向对齐数据.xlsx 是采集端用 getdate 获取时间戳后导出的船模试验记录通常包含时间、指令舵角、实际舵角、航向等通道。读入后先确认时间轴是否从 0 开始采集系统经常从操作时刻开始记录导致和仿真起点有固定偏移。处理方式是做一次线性插值把试验数据映射到仿真时间轴上raw readmatrix(数据.xlsx); % 老版本MATLAB用 xlsread t_exp raw(:,1) - raw(1,1); % 时间归零 psi_exp raw(:,3); t_sim (0:numel(psi_sim)-1) * dt; psi_exp_aligned interp1(t_exp, psi_exp, t_sim, linear);interp1 之前必须做 t_exp - t_exp(1) 归零否则插值区间起点不对结果要么直接报错要么外推出大量 NaN。读列号之前先打开 Excel 确认每列物理量航向若是度的单位要记得转弧度这步错了后面误差指标全是废的。5.2 误差指标与参数校核技巧对齐后计算航向的均方根误差和最大绝对误差err psi_sim - psi_exp_aligned; RMSE sqrt(mean(err.^2)); maxErr max(abs(err));对 Z 形试验RMSE 能压到 1° 以内就算相当好。若 RMSE 偏大先看误差曲线是整体偏移还是单个超越角位置凹陷整体偏移是时间对齐问题超越角处的局部偏差才需要动水动力导数。另一个实用技巧是分别统计右舵段和左舵段的 RMSE船模左右不完全对称时两侧误差会明显不同这时就要检查 Y 类导数里有没有遗漏不对称项。提示每改一次参数就重跑完整仿真太慢。可以先跑短时长的 Z 形试验前两个周期调阻尼类导数定下数量级后再跑完整回转试验验证回转直径整个调参过程能省一半以上的时间。本文还有配套的精品资源点击获取