在海洋资源探测、地形地貌测量及军事反潜作战中，需要UUV在某一特定海域进行探测或监测工作，为保证UUV维持其空间位置并保持良好的姿态，必须时刻通过推进器对UUV进行控制，从而消耗掉大量能源，无法达到长时间监测或探测的目的。为了解决UUV续航能力差、工作时间短的局限性，通过锚链将UUV与锚块连接起来，组成驻留UUV锚泊系统，使UUV处于水下一定深度并维持其静平衡状态^[1]，减少能源消耗。为便于UUV水下探测或监测工作的开展及其二次启动发射的实现，对UUV锚泊状态时的姿态有一定要求，必须研究系统在海流作用下的动态运动响应。

水下锚链、缆索广泛应用于拖曳及锚泊系统中，在非均匀海流作用下其动态运动解析求解十分困难，只能采用数值方法得到时域内的动态运动响应，有限差分方法则是其中最为普遍采用的方法。文献^[2]对二维状态下水面船舶的锚泊状态动态响应进行了求解分析；文献^{[3, 4]}分别研究了水下锚系导弹发射系统在海流作用下的三维运动情况以及通讯缆索对水下航行器运动的影响，但都没有考虑弯矩对锚链姿态的影响；文献^{[5, 6]}将缆索弯曲刚度考虑在内，建立了三维状态下拖曳系统的耦合运动模型，并采用有限差分方法求解。

与拖曳系统数学模型相比，锚泊系统的数学模型与之类似却有不同之处。两者不仅具有截然不同的边界条件，而且更值得注意的是，锚泊系统与海流之间的相对速度较低，海流非均匀性造成的锚链流体动力非线性更加显著，锚链易呈现低应力状态，锚链弯曲刚度的影响就变得不可忽视。另外，由于锚泊系统的下端点固定，UUV由于具有较大的运动惯性，其运动状态的改变往往滞后于锚链，这使得某些特殊情况下锚链姿态大幅度变化、个别欧拉角具有不确定性，会导致的运动方程奇异。

本文采用有限差分方法建立包含弯矩在内的锚链运动数学模型，使用四元数代替欧拉角，避免数值求解过程产生奇异，并通过耦合条件把驻留UUV、锚链与锚块系统综合考虑，建立适合驻留UUV锚泊系统的动态运动方程，对锚泊系统不同条件下的动态运动响应进行预报分析。

1 锚泊系统三维运动数学模型 1.1 坐标系选择

为便于分析锚泊系统在不规则海流作用下的动态运动响应，需要建立两个直角坐标系：地面坐标系S_E(O₀,x₀,y₀,z₀)、航行器体坐标系S_B(O,x,y,z)。如图1所示，地面坐标系原点选取在锚点O₀处，O₀x₀轴位于水平面内，O₀y₀轴位于竖直面内，铅直向上为正，O₀z₀轴的指向参照右手系规则确定；航行器体坐标系原点位于UUV浮心O处，Ox轴沿航行器纵轴且向后为正，Oy轴垂直于Ox轴并指向上方，Oz轴垂直于Oxy平面。另外，引入局部坐标系S_T(t,n,b)，其中轴t表示锚链切向，方向为锚链长度s的收缩方向；n为锚链的法向；b为锚链的副法线方向且位于水平面内。锚块沉于海底并与锚链下端固联，锚链上端在锚系点处与UUV连接，构成完整的锚泊系统。

坐标系之间可通过相应的姿态角相互转换，欧拉角(ε,γ)为锚链微元相对于地面坐标系的姿态角。其中ε为方位角，即锚链偏离X₀轴的角度，而γ为抬升角，如图2所示。它们可以定义为

$$\varepsilon = \left\{ {\matrix{ {\arcsin \left( {{d_{{z_0}}}/\left( {{d_l}\cos \gamma } \right)} \right)} & {{d_{{x_0}}} \ge 0} \cr {\pi - \arcsin \left( {{d_{{z_0}}}/\left( {{d_l}\cos \gamma } \right)} \right)} & {{d_{{x_0}}} \ge 0} \cr } } \right.$$

(1)

γ=arcsin(dy0/dl)

(2)

式中:d_l为锚链微元的长度，d_x0、d_y0、d_z0为微元在地面坐标系下的3个分量。通过坐标系旋转，可以获得局部坐标系到地面坐标系的转换矩阵： $$C_E^T = \left[ {\matrix{ {\cos \varepsilon \sin \gamma } & { - \cos \varepsilon \sin \gamma } & { - \sin \varepsilon } \cr {\sin \gamma } & {\cos \gamma } & 0 \cr {\sin \varepsilon \cos \gamma } & { - \sin \varepsilon \sin \gamma } & {\cos \varepsilon } \cr } } \right]$$

1.2 姿态四元数与欧拉角映射关系

为描述锚链微元的空间姿态，可以采用式(1)~(2)中描述的欧拉角方法或者采用四元数及其导出形式。欧拉角方法中采用空间旋转造成的人为割裂给姿态描述及姿态控制带来了便利，但其本身也存在许多不足，其过多的三角运算影响计算速度与精度，同时也不适用于大幅度的姿态运动描述，在特定情况下，运动方程出现奇异现象^[7]。为此，本文引入哈密尔顿四元数方法描述锚链微元的空间姿态。

由文献^{[7, 8]}可知，定点运动刚体由某一位置到另一位置的任意有限转动可以用四元数转动表示。从而局部坐标系到地面坐标系的转换矩阵为 $$C_E^T = \left[ {\matrix{ {q_0^2 + q_1^2 - q_2^2 - q_3^2} & {2\left( {{q_1}{q_2} + {q_0}{q_3}} \right)} & {2\left( {{q_1}{q_2} + {q_0}{q_2}} \right)} \cr {2\left( {{q_1}{q_2} + {q_0}{q_3}} \right)} & {q_0^2 - q_1^2 + q_2^2 - q_3^2} & {2\left( {{q_2}{q_3} - {q_0}{q_1}} \right)} \cr {2\left( {{q_1}{q_2} - {q_0}{q_2}} \right)} & {2\left( {{q_0}{q_1} - {q_2}{q_3}} \right)} & {q_0^2 - q_1^2 - q_2^2 + q_3^2} \cr } } \right]$$ 式中：哈密尔顿四元素q₀、q₁、q₂、q₃与姿态角ε、γ具有如下关系： $$\eqalign{ & {q_0} = \cos {\gamma \over 2}\cos {\varepsilon \over 2},{q_1} = \sin {\gamma \over 2}\sin {\varepsilon \over 2} \cr & {q_2} = \cos {\gamma \over 2}\cos {\varepsilon \over 2},{q_3} = \sin {\gamma \over 2}\sin {\varepsilon \over 2} \cr} $$

1.3 锚链运动学模型

本文建立的数学模型中，假设锚链为连续的细长圆柱状缆索，材质均匀且具有各向同性，在整个锚链长度上光滑连续。锚链微元的运动控制方程为^[9]

F′+CETW+H=Ma${\ddot r}$

(3)

Q′+r′×F=J${\dot \omega }$

(4)

式中：F为局部系下的锚链微元受力，Q为局部系下的锚链微元内部力矩，H为锚链微元的流体阻力，W为惯性系下单位长度锚链湿重，M_a为包括附加质量在内的单位长度锚链质量矩阵，r为局部系中的锚链微元位移向量，J为单位长度锚链的转动惯量，ω为锚链微元的角速度向量。控制方程中的运算符号定义如下，表示对括号内部的变量进行等式右侧的计算：

$\mathop {()}\limits^. $=∂/∂t+ω×()
()′=∂/∂s+κ×()

由于锚链沿长度s方向具有二阶光滑连续性，可以得到

$$\mathop {\left( {{{\partial r} \over {\partial s}}} \right)}\limits^. = {\left( {{{\partial r} \over {\partial s}}} \right)^\prime }$$

(5)

由坐标系之间的转换关系以及锚链微元的转动角速度、曲率的物理含义可知：

κ=κtκnκbT=Λ∂χ/∂s

(6)

ω=ωtωnωbT=Λ∂χ/∂t

(7)

式中：χ=εγ0^T；Λ为坐标系投影产生的矩阵，其表达式为 $$\Lambda = \left[ {\matrix{ {\sin \gamma } & 0 & 1 \cr {\cos \gamma } & 0 & 0 \cr 0 & 1 & 0 \cr } } \right]$$

综合式(3)~(7)，把锚链微元的运动控制方程分解到局部系t、n、b方向并写为矩阵形式^[10]：

$$M{{\partial Y} \over {\partial s}} = N{{\partial Y} \over {\partial t}} + P$$

(8)

式中：Ys,t=[F_t F_n F_b V_t V_n V_b κ_n κ_b γ ε]为锚链微元的控制变量；M、N均为10×10矩阵，矩阵M中的非零元素如下所示： M(i,i)=1,i=1，2，…10，M(1,9)=－F_n，M(1,10)=CγF_b，M(2,9)=F_t，M(2,10)=－SγF_b，M(3,10)=SγF_n－CγF_t，M(4,9)=－V_n，M(4,10)=CγV_b，M(5,9)=V_t，M(5,10)=－SγV_b，M(6,10)=SγV_n－CγV_t
矩阵N中的非零元素如下所示： N(1,4)=N(2,5)=N(3,6)=ρ_cA， N(1,9)=－ρ_cAV_n，N(1,10)=ρ_cACγV_b，N(2,9)=ρ_cAV_t，N(2,10)=－ρ_cASγV_b，N(3,10)=ρ_cA(SγV_n－CγV_t)，N(4,1)=1/EA，N(5,9)=1+F_t/EA，N(6,10)=－(1+F_t/EA)Cγ $$P = \left[ {\matrix{ {{1 \over 2}{C_t}{\rho _w}{d_c}{U_t}\left| {{U_t}} \right| - \sin \gamma \left( {{\rho _w} - {\rho _c}} \right)gA} \cr {{1 \over 2}{C_t}{\rho _w}{d_c}{U_n}/\sqrt {{U_n}^2 + {U_b}^2} - \cos \gamma \left( {{\rho _w} - {\rho _c}} \right)gA} \cr {{1 \over 2}{C_n}{\rho _w}{d_c}{U_b}/\sqrt {{U_n}^2 + {U_b}^2} } \cr 0 \cr 0 \cr 0 \cr {{\kappa _t}{\kappa _b} + {{\left( {1 + {F_t}/EA} \right){F_b}} \over {EI}}} \cr { - {\kappa _t}{\kappa _n} + {{\left( {1 + {F_t}/EA} \right){F_n}} \over {EI}}} \cr {{\kappa _b}} \cr {{\kappa _n}/\cos \gamma } \cr } } \right]$$ 式中:下标t、n、b表示局部系中各轴向,C_t、C_n为锚链的切向和法向阻力系数,A为锚链横截面积,E为锚链的杨氏弹性模量,ρ_w为海水密度,ρ_c为锚链密度,U为锚链微元相对海流速度。

1.4 边界条件

要求解方程(8)，必须给定锚链的初始形态值以及锚链上下两端点处的边界条件，形成封闭的系统运动模型。

锚链下端点处与锚块通过铰接方式相连接，在该点处弯矩为零，且锚点固定不动，可知

κn|s=0=κb|s=0=0

(9)

Vt|s=o=Vn|s=o=Vb|s=o=0

(10)

锚链上端点处与航行器铰接耦合，且速度与航行器耦合点处速度相同：

κn|s=L=κb|s=L=0

(11)

$$C_T^E\left[ {\matrix{ {{V_{t|S = L}}} \cr {{V_{n|S = L}}} \cr {{V_{b|S = L}}} \cr } } \right] = C_B^E\left[ {\matrix{ {{\upsilon _{xM}}} \cr {{\upsilon _{yM}}} \cr {{\upsilon _{zM}}} \cr } } \right] = \left( {\left[ {\matrix{ {{\upsilon _x}} \cr {{\upsilon _y}} \cr {{\upsilon _z}} \cr } } \right] + \left[ {\matrix{ {{\omega _x}} \cr {{\omega _y}} \cr {{\omega _z}} \cr } } \right] \times {r_A}} \right)$$

(12)

式中：C_B^E为UUV体坐标系到地面坐标系的转换矩阵，[v_xM v_yM v_zM]^T为体坐标系下锚链与航行器耦合点处速度，r_A为耦合点在UUV体坐标系下的坐标矢量。

2 锚泊系统运动模型数值求解方法

本文选用有限差分方法在时间和空间上对锚链微元的控制方程(8)~(12)离散化。通过n+1个节点将锚链划分为n段长度为Δs的微元，并将时间划分为一系列时间步长Δt。应用具有二阶精度的中心差分格式，将方程(8)在中间节点j+1/2和时间节点i+1/2处展开，可得到 $$\eqalign{ & \left[ {M_{j + 1}^{i + 1} + M_j^{i + 1}} \right]{{Y_{j + 1}^{i + 1} - Y_j^{i + 1}} \over {\Delta {s_j}}} + \left[ {M_{j + 1}^i + M_j^i} \right]{{Y_{j + 1}^i - Y_j^i} \over {\Delta {s_j}}} = \cr & \left[ {M_{j + 1}^{i + 1} + M_{j + 1}^i} \right]{{Y_{j + 1}^{i + 1} - Y_{j + 1}^i} \over {\Delta {t_j}}} + \left[ {N_j^{i + 1} + N_j^i} \right]{{Y_j^{i + 1} - Y_j^i} \over {\Delta {t_j}}} + \cr & P_{j + 1}^{i + 1} + P_j^{i + 1} + P_{j + 1}^i + P_j^i \cr} $$

应用Newton-Raphson方法迭代求解差分得到的非线性方程组(13)，即可得到锚泊系统的动态运动响应情况。

3 模型验证与数值计算结果 3.1 模型验证

为验证本文建立的包含弯矩作用在内的锚链数值模型精确性，选取Hopland拖曳实验^[11]与本文模拟数值进行对比。试验参数：锚链直径0.033 2 m,锚链总长360 m,锚链密度3 121 kg/m³,锚链弹性模量77.5 GPa。拖曳过程中，拖船速度历时60 s，由2.5 kn线性降至1.0 kn。不同时刻锚链姿态如图3所示。可以看到，模型数值模拟得到的锚链姿态动态值与实验值吻合良好，锚链阻力系数的不确定性造成了模拟值与实验值存在细微偏差。

3.2 数值计算结果分析

基于上述的锚泊系统运动数学模型，仿真UUV与锚链在海流作用下的运动过程。本文涉及的驻留航行器具有回转体外形，锚泊系统的主要参数如表1所示。UUV运行至距海底高度10 m处释放锚链，锚块及锚链在重力作用下下落，实现锚泊驻留功能，初始时刻锚链姿态呈竖直状态。为模拟海流的复杂缓慢变化，引入如下正弦波模型模拟海流速度J(T为海流运动周期)：

J=－1+0.5×sin(2πTt)，T=60 s

驻留系统的动态响应过程如图4~9所示。

从仿真结果可以看出，在海流作用下，锚泊系统位置发生偏移，由于锚块固联于海底，且驻留UUV相对于锚链具有较大的运动惯性，故而其运动状态变化滞后于锚链，如图4所示，锚链中间部分在锚泊运动初期呈凸起状，随着时间推移，航行器海流方向速度逐渐增加，锚链随之呈现悬链状。

由于海流波的周期性特征，由图所示，航行器的姿态偏降、速度、攻角及锚链拉力均以海流波动周期波动，并且比海流提前一定的相位角；航行器的俯仰角在海流及辅助推进器的作用下，被限制在[5°，-3°]范围内，可以保证良好的驻留姿态，便于水下工作的开展。

4 结论

1)与试验数据对比后发现，本文所建立的锚链三维运动数学模型不仅能有效避免数值求解过程中产生的奇异现象，同时还具有较高的计算精度，能准确地反映出锚链等拖曳系统与锚泊系统中常用的水下缆索张力传递及其自身形态变化情况。

2)驻留UUV的水下锚泊姿态与运动响应是整个锚泊系统设计的关键内容，对系统的运动响应进行仿真分析后发现，航行器俯仰角这一关键参数在海流作用下逐渐增大，在采用辅助推进器进行姿态调节后能明显将其限制在安全范围内，避免发生系统走锚或航行器触底现象。同时，这些系统的参数变化规律可以为航行器姿态控制和稳定性研究提供理论依据，为实现锚泊系统正常工作及驻留UUV二次启动提供有意义的参考。