引 言

我国要实现向航运强国的目标，必须更加注重船舶操纵性是否能满足安全上的要求，这对国内学者关于船舶航行和操纵安全的研究提出了更严格的需求。近年来，全回转尾推进型拖轮（Azimuth Stern Drive Tug，ASD 拖轮），由于配备双导管桨，能在水平面内绕垂直轴 360° 旋转，操纵灵活、马力强劲、安全可靠，所以其功能任务具有多样性，在港口辅助操作、海洋平台拖航、安全救助等任务中，起到了不可或缺的主力军作用。国内外学者对其操纵运动特性进行了大量研究，尚未有理想的数学模型能适用于四象限全工况操纵运动的计算与仿真。

对培训拖轮驾驶员来说，仅依靠经验理论掌握操纵性还不够，全回转拖轮操纵模拟器^[1] 对航海人才的教育培训、科学研究具有重要价值。在拖轮操纵模拟器环境中，存在大量静态或动态的障碍物，同样要求全回转导管桨工作在四象限，必要时须将桨偏角改变 180° 实现紧急停船。完整的数学模型能有效提高全回转拖轮模拟器的行为真实感^[2]。仿真质量的提高，提高培训效率及成果的同时亦能节约资金投入。本文着重介绍一种新的导管桨四象限推力模型，同时提出的阻力模型能满足于全速域要求，最后通过 Visual Studio 2013 建立测试平台，开展一系列操纵性仿真试验^[3] 研究，包括速度试航试验、螺旋线试验、Z 形试验^[4]、回转试验和停船试验^[5]，覆盖了 ASD 拖轮运动从小桨偏角到大桨偏角，常速域到低速域的运动全程。

1 ASD 拖轮操纵运动数学模型

假设 ASD 拖轮为刚体，在深、广、静、稳的水域中航行。ASD 拖轮运动采用定系 ${O_{\rm{0}}}{x_{\rm{0}}}{y_{\rm{0}}}{z_{\rm{0}}}$ 与动系Oxyz两个坐标系统。如图 1 所示，动系原点为 ASD 拖轮重心。按照 MMG（日本船舶运动数学模型研讨小组 Ship Maneuvering Mathematical Model Group）建模思想，模型如下：

$\left\{ {\begin{aligned} & {(m+{m_x})\dot u - (m+{m_y})vr = {X_H}{\rm{+}}{X_P}}\text{，}\\ & {(m+{m_y})\dot v+(m+{m_x})ur = {Y_H}+{Y_P}}\text{，}\\ & {({I_{zz}}+{J_{zz}})\dot r = {N_H}+{N_P}}\text{。} \end{aligned}} \right.$

(1)

式中：m，m_x，m_y 为 ASD 拖轮质量在x 轴和y 轴的附加质量；I_zz 为绕z 轴的质量惯性矩；J_zz 为附加质量惯性矩；u，v，r 为纵向、横向速度和回转角速度；X_H，Y_H，N_H 为船体水动力和首摇力矩；X_P，Y_P，N_P 为导管桨水动力和力矩。

ASD 拖轮操纵运动的轨迹和姿态方程如下，运用四阶龙格-库塔法解算：

$\left\{ {\begin{aligned} & {{{\dot x}_0} = u\cos \psi - v\sin \psi }\text{，}\\ & {{{\dot y}_0} = u\sin \psi+v\cos \psi }\text{，}\\ & {\dot \psi = r\cos \psi }\text{。} \end{aligned}} \right.$

(2)

其中：x₀、y₀ 为位置坐标；Ψ 为航向。

2 拖船阻力模型 2.1 惯性力和粘性力

流体惯性水动力的计算，参考周昭明关于附加质量m_x、m_y 和附加质量惯性矩J_zz 的回归公式。井上模型在计算小漂角运动时的粘性水动力较理想：

$\left\{ {\begin{aligned} & {{X_H} = X(u)+{X_{vv}}{v^2}+{X_{vr}}vr+{X_{rr}}{r^2}}\text{，}\\ & {{Y_H} = {Y_v}v+{Y_r}r+{Y_{|v|v}}|v|v+{Y_{|v|r}}|v|r+{Y_{|r|r}}|r|r}\text{，}\\ & {{N_H} = {N_v}v+{N_r}r+{N_{{\rm{|}}r|r}}|r|r+{N_{vvr}}{v^2}r+{N_{vrr}}v{r^2}}\text{。} \end{aligned}} \right.$

(3)

其中：“′”为 MMG 模型提出的一撇系统进行的无量纲化。关于水动力导数，采用井上公式计算线性项 $Y_v'$，$Y_r'$，$N_v'$，$N_r'$；采用周昭明公式计算非线性项 $Y_vv'$，$Y_vr'$，$Y_rr'$，$N_rr'$，$N_vvr'$，$N_vrr'$。

其中纵向水动力导数的模型如下：

${X'_{{H}}} = X'(u')+{X'_{vv}}{v'^2}+{X'_{vr}}v'r'+{X'_{rr}}{r'^2}\text{。}$

(4)

其中：${X'_{vv}}$，${X'_{vr}}$ 与 ${X'_{rr}}$ 各项计算参考文献[6]。

2.2 直航阻力的计算

ASD 拖轮直航阻力的计算采用二因次法：

$X'(u') \!=\! {X'_{uu}}{u'^2} \!=\! - \frac{S}{{Ld}}{C_{\rm{t}}}{u'^2} \!=\! - \frac{S}{{Ld}}\left( {{C_{\rm{f}}} \!+\! {C_{\rm{r}}} \!+\! \Delta {C_{\rm{f}}}} \right){u'^2}\text{。}$

(5)

式中：S 为湿表面积；L 为船长；d 为吃水；C_t 为总阻力系数；C_f 为平板摩擦阻力系数；C_r 为剩余阻力系数；ΔC_f 为粗糙度补贴系数。

直航阻力计算的关键在于剩余阻力的确定。图 2 为塔格尔特拖轮剩余阻力图谱，专门用于拖轮剩余阻力系数的计算。

由图 2 可知，图谱在傅汝德数Fr 为 0～0.18（船速 0～5.6 kn）的范围内未给出曲线，C_r 与Fr 在Fr 为 0.18～0.25 之间已近似成线性关系，在 0～0.18 之间曲线曲率更小，为避免拟合时C_r 出现负值，对该段曲线作线性处理；在Fr 大于 0.18 时，采用 5 阶多项式拟合。最终的 ASD 拖轮剩余阻力系数图谱如图 3 所示，分段拟合后可用于计算机仿真，能适用于全速域条件下的 ASD 拖轮操纵运动要求。

3 导管桨推力模型 3.1 联合推力模型

ASD 拖轮操纵运动时，2 个导管桨操作方向相同或相反，横向力及转矩大小相同，方向相同或相反。其联合推力模型如下：

$\left\{ {\begin{aligned} & {{X_P} = {T_{(p)}}cos({\delta _{(p)}})+{T_{(s)}}cos({\delta _{(s)}})}\text{，}\\ & {{Y_P} = - {T_{(p)}}sin({\delta _{(p)}}) - {T_{(s)}}sin({\delta _{(s)}})}\text{，}\\ & {{N_P} \!=\! ({T_{(p)}}cos({\delta _{(p)}}) \!-\! {T_{(s)}}sin({\delta _{(s)}})) \cdot {L_{(ps)}}/2\! - \!\!{Y_P} \cdot {x_p}}\text{。} \end{aligned}} \right.$

(6)

式中：p 和s 分别为左右导管桨；${L_{\left( {{{ps}}} \right)}}$ 为两导管桨之间的距离；$\delta $ 为导管桨的偏转角；$T$ 为导管桨的推力；${x_{{p}}}$ 为桨轴在动系中的纵向坐标值。

3.2 导管桨四象限推力模型 3.2.1 导管桨第一象限推力系数计算模型

导管桨有导流管和螺旋桨 2 部分组成，性能优良。导管桨总推力系数 ${K_{{T}}}$，为导流管推力系数 ${K_{{{TN}}}}$ 与螺旋桨推力系数 ${K_{{P}}}$ 的有效合成。但图谱资料只有第一象限曲线，由式（7）可看出，当紧急停船时，导管桨转速n 趋于 0，进速系数J 会趋于无穷大，超出图谱范围，仿真程序无法运行。

$\left\{ {\begin{aligned} & {T = (1 - {t_{\rm{p}}})\rho {n^2}{D^4}{K_{\rm{T}}}(J)}\text{，}\\ & {J = (1 - {w_{\rm{p}}})u/(nD)}\text{，}\\ & {{K_T} = \sum {\sum {{A_{ij}}{{(H/D)}^i}{{(J)}^j}} } }\text{，}\\ & {{K_P} = \sum {\sum {{B_{ij}}{{(H/D)}^i}{{(J)}^j}} } }\text{，}\\ & {{K_{TN}} = \sum {\sum {{C_{ij}}{{(H/D)}^i}{{(J)}^j}} } }\text{。} \end{aligned}} \right.$

(7)

式中：ρ 为海水密度；n 为桨转速；D 为桨直径；J 为进速系数；${t_{{p}}}$ 为推力减额；${w_{{p}}}$ 为伴流系数；$H/D$ 为螺距比；${A_{ij}}$、${B_{ij}}$、${C_{ij}}$ 分别为导管桨总推力、螺旋桨推力、导流管推力多项式系数。

在简易导管桨图谱资料中，上海交大设计的 JD75+Ka 导管桨系列曲线形式直观，性能优良。本文以图 4 导管桨为原型（4 为桨叶数，70 为盘面比）。图 4 中x₀，y₀ 为记下的特殊点坐标值。

3.2.2 四象限螺旋桨推力系数计算模型

有界形式的四象限螺旋桨工作特性图谱，突破了实际船舶四象限全工况表达上的困难，但未考虑带有导流管的情况。

$\left\{ \begin{array}{l} P = (1 - {t_{\rm{p}}})\rho {D^2}({v_{\rm{P}}}^2+{n^2}{D^2}){K_{\rm{p}}}^\prime (J')\text{，}\\[5pt] {v_{{P}}} = (1 - {w_{\rm{P}}})u\text{，}\\[5pt] J' = {v_{\rm{P}}}/\sqrt {{v_{\rm{P}}}^2+{n^2}{D^2}} = \left\{ \begin{array}{l} J\sqrt {1+{J^2}} ;\;\;\;\;n > 0{\rm{ }}\text{，}\\ - J/\sqrt {1+{J^2}} ;\;n < 0{\rm{ }}\text{，} \end{array} \right.\\[5pt] {K_{{p}}}^\prime (J') = P/\rho {D^2}({v_{{P}}}^2+{n^2}{D^2}) = {K_{{p}}}/1+{J^2}\text{。} \end{array} \right.$

(8)

式中： ${K'_{{p}}}\left( {J'} \right)$ 为螺旋桨通用推力系数；J′ 为通用进速系数；$v_p$ 为螺旋桨相对水的速度。

在四象限螺旋桨推力图谱^[7] 资料中，哥德堡水池试验得出的诺尔特斯特洛姆^[8]Nordstrom 系列应用广泛。该图谱适合螺距比 $H/D$ 从 0～1.6，桨叶数为 4。图 5是对其部分曲线进行切比雪夫 Chebyshev 多项式拟合^[9]所得，图 5中x₁，y₁ 为记下的特殊点坐标值。

3.2.3 导管桨四象限推力系数计算模型

对图 4 用到的的推力曲线 ${K_{{T}}} - J$ 处理如下：

${x^\prime } = x/\sqrt {1+{x^2}} ;{y^\prime } = y/1+{x^2}\text{，}$

(9)

则处理后得出如图 6 中的 ${K'_{{T}}} - J'$ 推力曲线，作为导管桨四象限中的第一象限推力曲线，相应特殊点的坐标值变为 ${x_0}^\prime $ 和 ${y_0}^\prime $。

对图 5 用到的 ${K'_{\rm{P}}} - J'$ 曲线伸缩变换处理如下：

$\left\{ \begin{array}{l} \!\!\!\! a = {{y'}_0}/{y_1}\text{，}\\ \!\!\!\! {b_1} = {{x'}_0}/{x_1};\;\;\;\;\;\;\;\;\;\;\;{b_2} = \left( {1 - {{x'}_0}} \right)/\left( {1 - {x_1}} \right)\text{，}\\[5pt] \!\!\!\! n > 0,{\rm{ }}0 \leqslant J' < {x_1}{\rm{;}}\;\;x' = {b_1}x,{\rm{ }}y' = ay\text{，}\\[5pt] \!\!\!\! n \!\!> \!\!0,{\rm{ }}{x_1} \leqslant J'\!\! \leqslant \!\!1;\;\;\;x' \!\!=\!\! {{x'}_0}+{b_2}(x - {x_1}),{\rm{ }}y' \!\!=\!\! ay\text{，}\\[5pt] \!\!\!\! n > 0, - 1 \leqslant J' < 0;\;\;x' = x,{\rm{ }}y' = ay\text{，}\\[5pt] \!\!\!\! n < 0, - 1 \leqslant J' \leqslant 1;\;\;x' = x,{\rm{ }}y' = ay\text{。} \end{array} \right.$

(10)

则处理后得出图 7 中的第二、三和四象限推力曲线 ${K'_{{T}}} - J'$，综合后的全方位导管桨四象限推力曲线即如图 7 所示 。

式（11）为最终的导管桨四象限推力系数计算模型。该模型既考虑了导流管的作用，又符合四象限全工况实际工作特性；既避免了仿真严重失真，又突出了定常正倒车之间、正车到倒车与倒车到正车过程的不对称性。

$T = (1 - {t_{{p}}})\rho {D^2}({v_{{P}}}^2+{n^2}{D^2}){K_{{T}}}^\prime (J')\text{。}$

(11)

4 操纵性试验及分析

表 1 为船桨参数。仿真对比资料是从丹麦 Force Technology 公司购买的全回转拖轮运动数学模型时配套提供的试验数据（记为 DMI 数据），但模型内部结构和参数看不到，原型船为 ASD 拖轮 Delta Linda。

图 8为使用 Visual Studio 2013 软件构建的 Modeler 仿真平台，其既能实时显示船舶的位置、航向及速度信息，又能对指定数据在指定时间段内进行存取。

4.1 速度试航试验

表 2 显示误差值最大为 4.7%，仿真效果良好，证明该模型在直航状态下较为理想。

4.2 螺线试验

从螺线试验可看出各桨偏角下的定常回转角速度r，船舶是否直线稳定以及稳定性的好坏，不具有直线稳定性船舶的不稳定程度等操纵性能。如图 8其试验运动轨迹形如螺线。试验步骤：先让拖船保持匀速直航，调整桨偏角 $\delta $ 作回转运动，等r 变为常数时，记录r、$\delta $；再次调整 $\delta $ 作回转运动，等r 变为常数时，记录r、$\delta $；试验中依次将桨偏角r 调整为：15°→10°→8°→6°→4°→2°→1°→0°→-1°→-2°→-4°→-6°→-8°→-10°→-15°。

图 9 为 DMI 数据与仿真数据得出的r-$\delta $ 曲线。可以看出，仿真数据与 DMI 数据基本吻合，它是一条单调曲线，不存在滞后环，说明 ASD 拖轮具有直线稳定性；原点处切线的斜率较其他地方大，说明该 ASD 拖轮虽具有直线稳定性，但稳定性不好；在$\delta $ 较小时，r 仍很大，说明 ASD 拖轮具有良好的回转性能，牺牲掉了部分直线稳定性。

表 3 显示了螺线试验回转速降对比情况，DMI 数据与仿真数据最大误差为 3.3%，精度较高。两者桨偏角越大，漂角越大，纵向速度降幅越大，回转半径越小。进一步表明了 ASD 拖轮回转性能的优良，当受到环境干扰时，为保持航向，需要及时调整桨偏角，调整后又过于灵敏，与拖船操作实况相符。

4.3 Z 形试验

利用 Z 形试验结果可以计算操纵性指数，易于判别船舶操纵性能，特别是中、小桨偏角下保持和改变航向的能力。10° Z 形试验步骤：先让拖船保持匀速直航，尽快将桨偏角调至左 10°；当首向角左偏初始航向 10°时，立刻将桨偏角调至右 10°，记为第 2 次操作；当首向角右偏离初始航向 10°时，再向左操作；如此来回操作 5 次结束试验。20° Z 形试验同上，将桨偏角每次改为 20°即可。

图 10 分别给出了 10°和 20° Z 形试验的桨偏角和首向角的历时曲线。可看出，2 个 Z 形试验仿真数据与 DMI 实船数据变化基本一致，误差在合理范围内。其中，10°较 20°准确，后者较前者执行桨偏角稍有延迟。另外，执行桨偏角越大，首向角向左右舷各摆动一次的周期越大，拖船运动越偏离中心线，在实际操船中 10°桨偏角也更常用。在 10°Z 形试验图中，ASD 拖轮仿真计算的第一、第二超越角分别为 14.1°和 19.0°，第一、第二超越时间分别为 4.0s 和 4.8s，周期为 27.8s；相对应的 DMI 数据为 13.8°和 15.7°，3.0s 和 4.0s，周期为 27.7s。普通船舶超越角在 10°左右，该 ASD 拖轮超越角较大，表明其回转性良好，跟从性一般。

4.4 回转试验

旋回性能代表了船舶航向改变的能力。ASD 拖轮作旋回运动时，先直航前进，改变桨偏角并稳定，拖轮重心运动的轨迹称为旋回圈，其特征参数为判别 ASD 拖轮旋回性能好坏提供了依据。

旋回运动轨迹如图 11 所示。本次试验初始航向为正北，初始航速为 13.61kn，将桨偏角均向右偏 35°操作。相应特征参数的对比结果如表 4。

从图 11 和表 4 可知，完成旋回运动纵向和横向占用水域距离均小于 2 倍拖轮船长，体现了 ASD 拖轮优良的旋回性能。由表 4 对比数据可知，旋回初径的误差值最大，仍在合理范围内。

4.5 停船试验

普通单桨船倒车停船时，螺旋桨旋转方向改变，航迹不是直线，沿曲线轨迹作偏转运动。ASD 拖轮的倒车停船操纵有所差异，其在倒车停船时，螺旋桨旋转方向不变，把桨偏角改变 180°，航迹为“1”字直线型。下令停船到 ASD 拖轮船速为降为 0 时船舶纵向前进的距离，即为停船冲程。DMI 数据给出的倒车停船试验，在桨偏角旋转 180°的同时，分为主机转速不变和主机转速减半 2 种情况。仿真时采取完全相同的初始条件，图 12 给出了纵向距离历时曲线x₀ 和纵向速度历时曲线u 的对比情况。

从图 12 可知，在主机转速不变的条件下，停船冲程 DMI 数据为 58.4 m，仿真数据为 61.1 m，误差为 4.6%；停船冲时 DMI 数据为 16.0 s，仿真数据为 14.6 s，误差为 8.8%。在主机转速由全速变为半速的条件下，停船冲程 DMI 数据为 61.1 m，仿真数据为 66.8 m，误差为 9.3%；停船冲时 DMI 数据为 18.9 s，仿真数据为 16.8 s，误差为 11.1%。图中 DMI 与仿真纵向距离历时曲线变化趋势相当一致；纵向速度历时曲线对比中，DMI 数据与仿真数据变化趋势略有不同，两导管桨相对位置不同，干扰会有变化，表示该模型在两桨干扰方面有待改进。总体来看，2 组试验中的衡准参数误差均在可接受范围内，特别是左图 ASD 拖轮紧急停船的运动过程，一致性良好，进一步表明了所建模型的有效性。

5 结 语

基于操纵性理论和分离型建模思想，本文给出了 ASD 拖轮适用于全工况操纵运动的数学模型及其仿真平台。通过对 ASD 拖轮进行速度试航试验，表明该模型在常速域直航状态下的精度较高；通过螺旋线试验、Z 形试验和回转试验，体现了该模型具有保持和改变航向的能力；通过停船试验，验证了该模型从常速域到低速域均满足推阻力匹配要求。综上，本文所建船舶运动数学模型正确有效，对全回转拖轮操纵模拟器的开发与完善具有参考意义。