0 引　言

锚是船舶最重要的设备之一，除了保证船舶安全停靠之外，锚设备还可以辅助操纵船舶，因此，分析拖锚运动的动态变化对于预测和控制船舶的行为具有重要意义。目前，国内外对船舶运动模型在锚泊操纵方面的研究较少，对于船舶拖锚运动建模，更是少有人进行全面系统的研究。

锚链作为船舶拖锚运动中密不可分的一部分，其形式多样，受力情况复杂。近年来，已有多种数值计算模型用来求解锚链的张力与构型，其中包括动态计算方法（集中质量法^{[1 - 2]}、有限差分法^[3]）、稳态计算方法（悬链线法^{[4 - 5]}、有限元法^[6]、直接积分法^{[7 - 8]}等）。考虑到船舶获取锚链张力计算所需条件的局限性，本文采用直接积分法对锚链力进行求解。Chiou等^[7]、Sun等^[8]采用了直接积分法构建动力学方程，前者率先应用此方法深入探究了缆绳不规则运动的特性，后者在此基础上进行了进一步的改良和提升，其引入缆绳速度和姿态角作为关键变量，研究缆绳在收放过程中的水动力响应以及低应力条件下的动态特性。

本文利用直接积分法建立锚链力计算模型，将上述作用力与分离型模型（Manoeuvring Model Group，MMG）相耦合，建立船舶三自由度拖锚运动数学模型，使其能够获取拖锚运动过程中船舶的运动参数以及锚链的受力情况。

1 建立船舶三自由度数学模型

本文在描述船舶平面运动时采用如图1所示的坐标系统。该系统包含2个坐标系^[9]，一个是空间坐标系$ {O}_{0}-{x}_{0}{y}_{0} $，即惯性坐标系，其中$ x_0y_0 $平面与水面重合；另一个是随船运动坐标系$ O-xy $，即附体坐标系，$ O $取船的中心位置，规定$ x $轴指向船首，$ y $轴指向右舷；航向角$ \varPsi $定义为$ {x}_{0} $轴和$ x $轴之间的夹角；$ \beta $为漂角；$ r $为艏摇角速度；$ u、v $分别为$ x、y $方向上的速度分量，$ \theta $为锚链方位。

基于MMG理论，在船舶纵荡、横荡以及艏摇运动之间的互相作用下，耦合锚链的作用力和力矩，以船体中心位置作为附体坐标系的原点，建立船舶三自由度拖锚运动方程：

$ \left\{ {\begin{aligned} &{(m + {m_x})\dot u - (m + {m_y})vr - m{x_G}{r^2} = } \\ &{{X_H} + {X_P} + {X_R} + {X_W} + {X_A} + {X_C} + {X_T}}，\\ &{(m + {m_y})\dot v + (m + {m_x})ur + m{x_G}\dot r = } \\ &{{Y_H} + {Y_R} + {Y_W} + {Y_A} + {Y_C} + {Y_T}} ，\\ &{({I_{zz}}{\text{ + }}{J_{zz}} + mx_G^2)\dot r + m{x_G}(\dot v + ur) = } \\ &{{N_H} + {N_R} + {N_W} + {N_A} + {N_C} + {N_T}} 。\end{aligned}} \right. $

(1)

式中：$ m $为船舶质量；$ {m}_{x} $、$ {m}_{y} $为船舶分别在$ x $、$ y $轴上的附加质量；$ {I}_{zz} $为船舶绕过船体重心轴的转动惯量；$ {J}_{zz} $为船舶附加质量惯性矩；$ {x}_{G} $为船舶重心坐标；$ \dot{u} $、$ \dot{v}、\dot{r} $分别为前进方向、横移方向和转动方向的加速度；$ \left({X}_{H},{Y}_{H},{N}_{H}\right) $为作用于船体上的力和力矩；$ {X}_{P} $为作用于螺旋桨上的力；$ \left({X}_{R},{Y}_{R},{N}_{R}\right) $为作用于舵上的力和力矩；下标$ A $、$ W、C、T $分别为风、浪、流以及锚链。其中，

$ {\left\{\begin{gathered}X_H=\left(\dfrac{\rho}{2}\right)Ld_m\left\{X'_0uU+\left(m_y+X'_{vr}+\Delta X'_{vr}\right)Lvr\right\}，\\ Y_H=\left(\dfrac{\rho}{2}\right)Ld_m\left[ \begin{array}{*{20}{c}}\left\{Y'_vv\left|u\right|+\left(Y'_r-m_x\right)Lru\right\}-\dfrac{C_D}{\boldsymbol{L}}\cdot \\ \int_{-\boldsymbol{L}/2}^{\boldsymbol{L}/2}\left(1-\tau'\dfrac{x}{\boldsymbol{L}}\right)\left|v+C_{rY}rx\right| \\ \left(v+C_{rY}rx\right){\rm{d}}x\end{array} \right]，\\ N_H=\left(\dfrac{\rho}{2}\right)L^2d_m\left[\begin{array}{*{20}{c}}\left\{N'_vvu+N'_rLr\left|u\right|\right\}-\dfrac{C_D}{L^2}\int_{-\boldsymbol{L}/2}^{\boldsymbol{L}/2}\cdot \\ \left(1-\tau'\dfrac{x}{\boldsymbol{L}}\right)\left|v+C_{rN}rx\right| \\ \left(v+C_{rN}rx\right)x{\rm{d}}x\end{array}\right]。\\\end{gathered} \right. }$

(2)

式中：$ \rho $为海水密度，取1025 kg/m³；$ L $为船舶质量；$ {d}_{m} $为船舶标准吃水；$ {X}_{0}' $为直航阻力系数；$ {Y}_{v}' $、$ {Y}_{r}' $、$ {N}_{v}' $、$ {N}_{r}' $分为横向和艏摇水动力导数；$ {X}_{vr}' $为纵向粘性水动力导数；$ {C}_{rY}、{C}_{rN}、{C}_{D} $均为模型修正系数；式(1)与式(2)中各参数及力的计算可参考文献[9 - 10]。

2 锚链与锚力的计算 2.1 锚链微分方程

为了更为方便地求解锚链张力及构型，本文建立了如图2所示的船舶单锚泊状态下的坐标系。图中：$ {O}_{0}-{x}_{0}{z}_{0} $为惯性坐标系；$ {x}_{0} $指向锚链延展方向；$ {z}_{0} $垂直指向水面；设锚位为坐标原点$ o $；$ \varphi $为锚链在$ o $处的欧拉角；$ W $为锚链水中的重力；$ T $为$ o $处锚链的张力；$ {L}_{c} $为锚链释放总长；$ {l}_{c} $为卧底锚链长度；$ {d}_{a} $为锚链孔距锚位水平距离；$ {h}_{a} $为锚链孔距锚位垂直距离；$ {T}_{s} $为锚链孔处锚链的张力；$ {\varphi }_{s} $为锚链孔处的欧拉角。

设锚链是细长、柔性的，不考虑扭矩的影响。锚链微元$ \mathrm{d}s $所受张力和水中重力分别为$ \overrightarrow{T} $和$ \overrightarrow{W} $。$ s $为未拉伸的锚链长度，$ S $为拉伸后的锚链长度。一般锚链材料的拉伸特性为$ {S}'=\partial S/\partial s=1+eT $。其中，$ e=1/EA $，$ E $为拖缆弹性模量；$ A $为拖缆横截面积。根据牛顿第二定律，锚链微元的静力平衡方程为：

$ \frac{{\partial \mathop T\limits^ \to }}{{\partial s}} = S'\mathop W\limits^ \to。$

(3)

将式(3)在锚链附体坐标系下沿各个坐标轴方向展开，同时根据惯性坐标系与附体坐标系的转换关系可得平衡方程的标量形式：

$ \left\{ \begin{gathered} {{{\text{d}}T} \mathord{\left/ {\vphantom {{{\text{d}}T} {{\text{d}}s}}} \right. } {{\text{d}}s}} = {w_c}g\sin \phi ，\\ {{{{{\text{d}}\phi } \mathord{\left/ {\vphantom {{{\text{d}}\phi } {{\text{d}}s}}} \right. } {{\text{d}}s}} = {w_c}g\cos \phi } \mathord{\left/ {\vphantom {{{{{\text{d}}\phi } \mathord{\left/ {\vphantom {{{\text{d}}\phi } {{\text{d}}s}}} \right. } {{\text{d}}s}} = {w_c}g\cos \phi } T}} \right. } T} ，\\ {{{\text{d}}x} \mathord{\left/ {\vphantom {{{\text{d}}x} {{\text{d}}s}}} \right. } {{\text{d}}s}} = \left( {1 + Te} \right)\cos \phi ，\\ {{{\text{d}}z} \mathord{\left/ {\vphantom {{{\text{d}}z} {{\text{d}}s}}} \right. } {{\text{d}}s}} = \left( {1 + Te} \right)\sin \phi 。\\ \end{gathered} \right. $

(4)

式中，$ {w}_{c} $为水中每米锚链重量。

2.2 锚链微分方程求解

锚链的构型随锚链孔与锚位间距离$ {d}_{a} $的不同而不同，如图3所示。$ {d}_{a1} $为锚链水平拉力$ {T}_{H}=0 $时的最大船-锚间距；$ {d}_{a2} $为卧底锚链$ {l}_{c}=0 $时的最小船-锚间距；$ {d}_{a3} $为最大船-锚间距。

当$ {d}_{a}\leqslant {d}_{a1} $时，锚链垂直于水平面，部分锚链可能随意铺设或堆积在海底，对船体不产生水平拉力；当$ {d}_{a1}\leqslant {d}_{a} < {d}_{a2} $时，存在部分卧底锚链，锚链在水底处的欧拉角$ \varphi =0 $，在锚链孔处$ {\varphi }_{s} > 0 $，对船体产生水平拉力；$ {d}_{a2}\leqslant {d}_{a} < {d}_{a3} $时，不存在卧底锚链，锚链在水底处的欧拉角$ \varphi $和在锚链孔处的欧拉角$ {\varphi }_{s} $均大于$ 0 $，对船体产生水平拉力。

因此，根据上述3种情况分别确定锚链两端的边界条件，然后利用二分法和龙哥库塔方法对锚链微分方程积分求解。

1）$ {d}_{a}\leqslant {d}_{a1} $时锚链的张力及构型

如释放锚链总长为$ {L}_{c} $，则在不考虑弹性形变的情况下可认为$ {d}_{a1}={L}_{c}-{h}_{a} $。在考虑锚链弹性形变的情况下，采用如下计算方法。根据锚链海底端的边界条件可知：

$ \left\{ \begin{gathered} {T_0} = 0 ，\\ {\phi _0} = 0 ，\\ {x_0} = {L_c} - {h_a} ，\\ {z_0} = 0。\\ \end{gathered} \right. $

(5)

将式(5)作为初始值代入式(4)中，采用二分法迭代求解$ x $，即$ {d}_{a1} $，直至满足：

$ \left| {\int_0^{{L_c} - x} {(1 + eT)\sin \phi {\text{ds}} - {h_a}} } \right| < \varepsilon。$

(6)

式中，$ {x}_{0}\in \left[{L}_{c}-{h}_{a},{L}_{c}-\left(1-{h}_{a}{w}_{c}\mathrm{g}e\right){h}_{a}\right] $。式(3)中$ T= 0 $为一奇点，可简单地将其设置为一个很小的数来取代0；$ \varepsilon $为设置的迭代误差。

可以认为，在$ {d}_{a}\leqslant {d}_{a1} $的情况下，锚链的张力和构型与$ {d}_{a}={d}_{a1} $时的情况一致。此时，虽然锚链孔端的锚链张力$ {T}_{s}\ne 0 $，但由于锚链孔端的φ_s = 90°，因此，锚链对船体不产生水平拉力。

2）$ {d}_{a1} < {d}_{a}\leqslant {d}_{a2} $时锚链的张力及构型

根据锚链海底端的边界条件可确定式(4)的积分初值如下：

$ \left\{ \begin{gathered} {T_0} = 0，\\ {\phi _0} = 0，\\ {x_0} = 0 ，\\ {z_0} = 0 。\\ \end{gathered} \right. $

(7)

当$ {d}_{a}={d}_{a2} $时，不存在卧底锚链，可根据锚链海底端张力的变化范围，采用二分法迭代求解$ {d}_{a2} $，直至满足：

$ \left|\displaystyle\int_0^{L_c}(1+eT)\mathrm{sin}\varphi\text{ds}-h_a\right| < \varepsilon\text{ }\text{，}\text{ }T_0\in\left(0,T_{\mathrm{\max}}\right]。$

(8)

当$ {d}_{a1} < {d}_{a} < {d}_{a2} $时，存在卧底锚链，此时，可根据卧底锚链长度，即$ {x}_{0} $的变化范围和锚链海底端张力$ {T}_{0} $的变化范围，采用二重二分法迭代求解锚链张力和构型，直至满足：

$ \left\{\begin{gathered}\left|\int_0^{L_c-x}(1+eT)\cos\phi\text{d}s-d_a\right| < \varepsilon，\\ \left|\int_0^{L_c-x}(1+eT)\sin\phi\text{d}s-h_a\right| < \varepsilon。\\ \end{gathered}\right. $

(9)

式中：$ {T}_{0}\in (0,{T}_{\max}] $, $ {x}_{0}\in [0,{d}_{a1}) $，奇点$ T=0 $的处理方法同式(1)；$ {T}_{\max} $可取锚链的最大破断力。

3）$ {d}_{a2} < {d}_{a}\leqslant {d}_{a3} $时锚链的张力及构型

此时，不存在卧底锚链，式(4)的积分初始值为式(7)。可根据锚位处$ {\varphi }_{0} $的变化范围和张力的变化范围，采用二重二分法迭代求解锚链张力和构型，直至满足下式：

$ \left\{ \begin{gathered} \left| {\int_0^{{L_c}} {(1 + eT)\cos \phi {\text{d}}s - {d_a}} } \right| < \varepsilon，\\ \left| {\int_0^{{L_c}} {(1 + eT)\sin \phi {\text{d}}s - {h_a}} } \right| < \varepsilon。\\ \end{gathered} \right. $

(10)

式中：$ {T}_{0}\in (0,{T}_{\max}] $，$ {\varphi }_{0}\in [0,a\mathrm{s}\mathrm{i}\mathrm{n}({h}_{a}/{L}_{c}\left)\right] $。令$ {T}_{0}={T}_{\max} $，联立式(4)、式(7)和式(10)即可迭代求解出$ {d}_{a3} $。

利用上述公式求解出锚链孔处的锚链张力$ {T}_{s} $和欧拉角$ {\varphi }_{s} $后可按下式计算作用在船体上的锚链力和力矩：

$ \left\{ \begin{gathered} {X_C} = {T_s}\cos {\phi _s}\cos \theta，\\ {Y_C} = {T_s}\cos {\phi _s}\sin \theta，\\ {N_C} = - {X_C} \cdot {y_C} + {Y_C} \cdot {x_C}。\\ \end{gathered} \right. $

(11)

式中：$({x_C},{y_C})$为锚链孔相对于船舶中心的位置。

2.3 船舶锚抓力的计算

锚的抓力主要由2个因素构成，即锚的自身抓力和锚链与海底的摩擦力。

锚的总抓力计算式为：

$ P = {P_a} + {P_c} = {\lambda _a} \cdot {w_a} + {\lambda _c} \cdot {w_c} \cdot {l_c} .$

(12)

式中：$ {P}_{a} $为锚的总抓力；$ {P}_{c} $为锚链的总抓力；$ {\lambda }_{a} $为锚抓力系数，可按表1和表2取值；$ {\lambda }_{c} $为锚链抓力系数，泥沙底质通常取0.6；$ {w}_{a} $为锚在空气中的重量，kN；$ {w}_{c} $为单位每米锚链在空气中的重量，N/m；锚的静抓力系数由多个因素决定，包括锚的种类、底质特性等。表1为运输船舶常用的几种锚型的静抓力系数^[11]。

在船失去锚固的情况下，抓力系数与拖锚时的动态锚抓力系数具有相同数值。当海底底质为泥沙时，锚的动抓力与锚链的出链长度及水深之间存在一定关系（见表2）。

锚链在锚链孔处的张力$ {T}_{s} $可由2.2节计算，当锚链的水平拉力$ {T}_{s}\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_{s} > P $时，按表2进行锚动抓力的计算。

3 数值仿真

为了验证模型的正确性和适用性，本文以“北斗”轮风中拖锚漂航实船试验为例进行风中拖锚漂航运动仿真。船舶主尺度及其锚链参数如表3和表4所示^[11]。

3.1 锚链构型与张力仿真计算

图4为200 m锚链随船锚间距的增加而变化的构型，图5为锚链孔端的锚链张力和水平交角以及卧底锚链长度随船锚间距的增加而变化的过程。可知，锚链作用于船体的张力随船锚距离的增加而增加，而水平交角和卧底锚链的长度则随之减小。在已知锚位、船位和出链长度的情况下，利用该算法，不仅可以预报锚链的张力和构型，而且还可计算出卧底锚链的长度，即可较准确地计算出锚和锚链的抓力，用于预报船舶是否走锚。在实际工作中，锚位、船位和出链长度较易获取，因此，该算法具有较强的实用性。

3.2 船舶操纵性仿真试验

为了验证所建船舶运动模型的准确性，本文以“北斗”轮为例进行旋回与停船仿真试验，并与实船试验结果^[12]进行对比验证。

1）旋回仿真试验验证

船舶旋回仿真试验的初始条件为：航速11.2 kn、航向0° 、舵角35° 、南风5.6 m/s。当仿真时长为400 s时，仿真结果如图6和图7所示。比对结果如表5所示。

2）停船仿真试验

船舶停船仿真试验的初始条件为：航速12 kn、航向0° 、主机转速0 、风速4.0 m/s，风向70°。仿真时长400 s，仿真结果如图8和图9所示。比对结果如表6所示。

由船舶旋回和停船仿真试验可知，所建模型操纵性符合度良好，误差值均小于10%，仿真结果符合实际，仿真精度满足工程需要。

3.3 风中拖锚漂航仿真试验

本文以“北斗”轮风中拖锚漂航实船试验为例^[13]进行风中拖锚漂航运动仿真，初始条件为：航速0 kn、横风风速15 m/s、水深22 m，泥底，锚链长2节。仿真试验时长8 min，仿真结果如图10和图11所示。比对结果如表7所示。

可知，由于船舶正横受风，导致船体横向漂移速度逐渐增加，最终由于锚链力及力矩和横流水动力及力矩的作用，横向漂移速度、锚链方位角和锚链张力逐渐趋于稳定；同时由于锚链力及力矩的作用使船舶最终的稳定状态既没有表现出固定锚泊船的偏荡运动，也没有表现出自由船舶风中的偏转规律，而是保持小于90°的风弦角向下风漂移，符合风中拖锚漂航的实际运动规律。

通过表7的相对误差可以看出，船舶拖锚漂航仿真试验结果与实船试验结果的符合度较好，与实船试验的最大误差为5.7%，能够满足工程及仿真的精度需求。通过对仿真结果的分析，更好理解了拖锚漂航的机制和影响因素，也进一步证明了本文拖锚运动数学模型的准确性。

4 结　语

本文基于直接积分法，提出了一种在已知锚位、船位及出链长度即可计算锚链张力和构型方法建立锚链与锚力的计算模型，并将此锚链力计算模型与船舶运动数学模型耦合，建立全面的三自由度船舶拖锚运动数学模型。同时，以“北斗”轮为实船研究对象，根据实船实验条件及数据，对本文提出的拖锚运动数学模型进行仿真验证，开展了拖锚漂航仿真试验以及不同风速下的漂航运动分析，仿真结果与实船数据的相对误差均在6%以下，表明了本文拖锚运动数学模型的正确性和适用性。