0 引　言

21世纪以来，以“蛟龙”号、“勇士”号、“奋斗者”号载人潜水器和“潜龙一号”、“潜龙二号”等无人潜水器为代表的水下高端装备使我国在深海空间的探索与开发领域取得了领先优势^[1-2]。

在潜水器设计中，安全性是一项至关重要的技术指标，也是设计人员的技术红线。潜水器设计与水面船舶设计的显著不同是其储备浮力小，通过注排水实现浮力调节，在紧急情况时通常抛弃部分载荷应急上浮。随着潜水器平台的大型化和多功能化要求，部分平台提出了水下铺缆、水下回收等多种作业新需求，对平台的水下起吊、拖曳回收能力有所要求。我国首个实验型35吨级深海移动工作站目前已进行了水池试验，开展了多项试验，初步对水下起吊和拖曳回收能力进行了探索性研究，其未来型将具备水下装备运载投送能力，并配备多种深海机器人，实现长期的水下科学考察与开发研究。

未来在大型多用途潜水器方案设计中，具备一定任务载荷的拖曳和起吊能力将成为一项值得关注的重要指标。拖曳和起吊作业时潜水器的航速低、舵效差，拖曳和起吊过程中缆索传递给潜水器的作用力变化可能对其安全性和作业性能有较大影响，必须对潜水器-缆索-载荷三者间相互耦合作用进行严谨的动力学仿真研究，分析缆索收放过程对缆索张力和潜水器运动相应的影响规律，从而更好为多用途潜水器的储备浮力设计策略和作业流程设计服务。

潜水器-缆索-载荷三者间相互耦合作用的动力学仿真研究问题的难点在于缆索的力学建模。目前针对缆索的仿真方法一般分为两类：一类是基于悬链线方程或缆索微元平衡方程对水下缆索力学问题进行静力学分析，典型用例为系泊系统设计中基于悬链线方程对缆索张力和曲率问题进行求解^[3-6]；一类是动力学分析，考虑其动态响应过程，按缆索微元的处理方式又分为集中质量法、有限差分法、有限段法、有限元法等。静力学分析不适用于潜水器-缆索-载荷三者间相互耦合作用的仿真分析，后者有一定的适用性，并已在潜艇拖曳线阵动态运动仿真运用中得到较好的仿真精度和效率。

目前上述方法基本未涉及收放缆索过程中变缆索长度、变边界问题的研究，或即使涉及到根据文献调研显示，也未取得较好的仿真效果（效率、精度、稳定性等）。

Ablow等^[7]首次将缆索动力学平衡方程与运动学方程投影变换至缆索局部运动坐标系，建立了缆索的偏微分控制方程组，开始基于有限差分法对动态缆索问题进行研究；Burgess^[8]和李英辉^[9]基于上述方法对深海拖曳系统在二维下拖缆与拖曳体的耦合运动进行了数值模拟。

付薇^[10]和王飞^[11]等采用集中质量法对拖缆系统的控制方法进行了研究；并且陆肇康等^[12]通过增加或减少凝聚质量点的数目，对非定长水下缆索动态运动（缆索收放）进行了初步研究。

Huston^[13]和Kamman^[14]等基于有限段法对水下缆索进行了动力学仿真；李晓平^[15]基于变拓扑多体系统理论对缆索收放过程中的变长度问题进行了初步研究。

有限元分析中一般将缆索采用梁单元或LINK单元近似，即使在Ansys2020正式加入了绳索单元CABLE280，但据了解其适用范围仍仅限于模拟线缆拉压、变形等工况，不适用缆索的大范围运动仿真。任革学团队开发了ALE缆索单元^[16]，并在缆索驱动多体系统仿真研究方面取得了很大进展，已在部分航天用大伸缩结构中得到运用，尚未清楚该方法是否适用于水下缆索收放仿真。

综上所述，采用集中质量法增加或减少凝聚质量点的数目、采用有限段法中增加或减少缆索段数目均能在一定程度上解决水下缆索收放的仿真问题，但凝聚质量点数目或缆索微元段数目的变化，将导致求解过程中质量矩阵、刚度矩阵等结构的不断改变，严重影响了求解效率和求解稳定性。

本文采用有限差分法、并结合缆索微元长度更新策略，开展潜水器-缆索-载荷三者间相互耦合作用研究中涉及的缆索变长度问题，并通过拖航变深度航行、拖航过程收缆索进行起吊作业2种工况的数值仿真对提出的方法进行验证。

1 仿真方法 1.1 坐标系定义

大地坐标系为 $E{{ - }}\xi \eta \zeta $ 、随体坐标系为 $ b{{ - }}xyz $ ，两坐标系原点重合， $E{{ - }}\xi \eta \zeta $ 坐标系依次绕z-y-x轴转动 $ \psi , $ $ \theta ,\phi $ 后与 $ b{ - }xyz $ 坐标轴重合。缆索局部坐标系为 $c{ - }tnb$ ， $t$ 为缆索的切线方向，规定沿缆索长度s增长的方向为正；b与t垂直并位于水平面内（ ${\psi _c} = \text{π} /2$ ）。 $ {\varphi _c} $ 为t轴与 $ij$ 平面夹角， $ {\theta _c} $ 为t轴在 $ij$ 平面的投影与 $i$ 轴的夹角。由大地坐标系 $E{{ - }}\xi \eta \zeta $ 依次绕z-x-z转动 $ - {\theta _c}, - {\psi _c},{\varphi _c} $ （欧拉角）获得。转换关系如下：

$ \left( {\begin{array}{*{20}{c}} {\boldsymbol{t}} \\ {\boldsymbol{n}} \\ {\boldsymbol{b}} \end{array}} \right) = \left[ {\begin{array}{*{20}{c}} {cos{\varphi _c}cos{\theta _c}}&{cos{\varphi _c}sin{\theta _c}}&{sin{\varphi _c}} \\ { - sin{\varphi _c}cos{\theta _c}}&{ - sin{\varphi _c}sin{\theta _c}}&{cos{\varphi _c}} \\ {sin{\theta _c}}&{ - cos{\theta _c}}&0 \end{array}} \right]\left( {\begin{array}{*{20}{c}} {\boldsymbol{i}} \\ {\boldsymbol{j}} \\ {\boldsymbol{k}} \end{array}} \right)。$

(1)

1.2 潜水器模型定义

潜水器的运动模型采用Fossen^[17]使用的符号表示： ${\boldsymbol {\eta }} = {\left[ {n,e,d,\phi ,\theta ,\psi } \right]^{\rm{T}}}$ , ${\boldsymbol{v}} = {\left[ {u,v,w,p,q,r} \right]^{\rm{T}}}$ 。其中 $ n,e,d $ 表示在 $E{\text{ - }}\xi \eta \zeta $ 坐标系下沿北、东、向下指向的位置， $ \phi ,\theta ,\psi $ 则表示相应的欧拉角（滚转、俯仰和偏航）； $ u,v,w, p, $ $ q,r $ 表示潜水器在随体坐标系 $ b{{ - }}xyz $ 下的广义速度。

潜水器的运动模型采用美国海军研究院的NPS AUV Ⅱ^[18-20]，该模型采用尾部2台正反转螺旋桨推进，在首部、尾部有多个姿态调节舵面，与我国首个实验型35t级深海移动工作站具备一定相似性。

NPS AUVⅡ 潜水器的运动方程为：

$ {{M}}\dot v + {{C}}\left( v \right)v + {{D}}\left( v \right)v + {{g}}\left( \eta \right) = {{Q}} 。$

(2)

其中： ${{M}}$ 为系统的惯性矩阵（包含附加质量）； ${{C}}\left( v \right)$ 为科里奥利向心力矩阵； ${{D}}\left( v \right)$ 为阻尼矩阵； ${{g}}\left( \eta \right)$ 为重力、浮力的力和力矩向量； ${\boldsymbol{Q}}$ 为外力和外力矩的合向量。其具体的表达形式见文献[19]。

1.3 缆索收放过程力学建模

Ablow等将缆索视为圆截面理想柔性弹性体，张力始终作用于缆索切向，忽略弯曲、扭转和剪切力的影响，且缆索处于张紧状态。重力加速度为 $g, \; {\text{km/}}{{\text{s}}^2}$ ；水流密度为 $\rho {\text{，}} {\text{kg/}}{{\text{m}}^{\text{3}}}$ ；缆索杨氏模量为 $ E $ ，缆索的等效直径为 $d \;{\text{m}}$ ，单位长度质量为 $m\; {\text{kg/m}}$ ，其法向和切向拖曳力系数分别为 $ {C_n} $ 和 $ {C_\iota } $ 。缆索系统的建模采用符号： ${{\boldsymbol{X}}_c} = {\left[ {x,y,z} \right]^{\rm{T}}}$ ， ${{\boldsymbol{Y}}_c} = {\left[ {T,{v_t},{v_n},{v_b},{\theta _c},{\varphi _c}} \right]^{\rm{T}}}$ 表示。其中 $ x,y,z $ 表示缆索上各点在大地坐标系 $E{{ - }}\xi \eta \zeta $ 的位置， $ {v_t},{v_n},{v_b} $ 为缆索上各点在随体坐标下的各个速度分量， $ T $ 为缆索上各点的缆索张力， $ {\theta _c},{\varphi _c} $ 为缆索上各点的姿态角。

大地坐标系 $E{{ - }}\xi \eta \zeta $ 到缆索局部坐标系 $ c - tnb $ 的变换矩阵位 ${\boldsymbol{R}}_n^c\left( {{\Theta _c}} \right)$ ，则潜水器收放缆索处的速度 ${{\boldsymbol{v}}_E}$ （在大地坐标系下），转换到局部坐标系下，该点的速度为 ${\boldsymbol{v}} = {\boldsymbol{R}}_n^c\left( {{\Theta _c}} \right){{\boldsymbol{v}}_E}$ 。海流速度在大地坐标系下的速度为 $ {{\boldsymbol{J}}_E}\; m/s $ ，缆索绞车装置的收放速度为 $ {{\boldsymbol{v}}_p} $ ，则缆索微元相对洋流的速度为：

$ {\boldsymbol{u}} = {\boldsymbol{R}}_n^c\left( {{\Theta _c}} \right)({{\boldsymbol{v}}_E} - {{\boldsymbol{J}}_E}) + {{\boldsymbol{v}}_p} ，$

(3)

结合文献[7]可得到考虑缆索收放速度变化的缆索偏微分方程组控制方程组：

$ {\boldsymbol{M}}\left( {s,t} \right)\frac{{\partial {\boldsymbol{Y}}\left( {s,t} \right)}}{{\partial s}} = {\boldsymbol{N}}\left( {s,t} \right)\frac{{\partial {\boldsymbol{Y}}\left( {s,t} \right)}}{{\partial t}} + {\boldsymbol{Q}}\left( {s,t} \right) 。$

(4)

式中：

$ {\boldsymbol{M}}\left( {s,t} \right) = \left[ {\begin{array}{*{20}{c}} 1&0&0&0&0&0 \\ 0&1&0&0&{ - {v_b}\cos \varphi }&{ - {v_n}} \\ 0&0&1&0&{ - {v_b}\sin \varphi }&{{v_t}} \\ 0&0&0&1&{{v_t}\cos \varphi + {v_n}\sin \varphi }&0 \\ 0&0&0&0&{T\cos \varphi }&0 \\ 0&0&0&0&0&T \end{array}} \right]，$

$ \begin{gathered} {\boldsymbol{Q}}\left( {s,t} \right) = \hfill \\ \left[ {\begin{array}{*{20}{c}} {w\sin \varphi + 0.5\sqrt {1 + eT} \rho \text{π} d{C_\iota }{u_t}\left| {{u_t}} \right|} \\ 0 \\ 0 \\ 0 \\ {0.5\rho d\sqrt {1 + eT} {C_n}{u_b}\sqrt {u_n^2 + u_b^2} - \rho A{{\dot J}_b}} \\ {wcos\varphi + 0.5\rho d{C_n}\sqrt {1 + eT} {u_n}\sqrt {u_n^2 + u_b^2} - \rho A{{\dot J}_n}} \end{array}} \right] \hfill \\ \end{gathered} ，$

$ \begin{split}{\boldsymbol{N}}\left( {s,t} \right) =& \left[ {\begin{array}{*{20}{c}} {\dfrac{{ - me{v_t}}}{{1 + eT}}}&m&0&0 \\ e&0&0&0\\ 0&0&1&0\\ 0&0&0&1\\ {\dfrac{{ - e({m_1}{v_b} - \rho {\rm A}{J_b})}}{{1 + eT}}}&0&0&{{m_1}}\\ {\dfrac{{ - e({m_1}{v_n} - \rho {\rm A}{J_n})}}{{1 + eT}}}&0&{{m_1}}&0 \end{array}} \right.\\ & \left. {\begin{array}{*{20}{c}} { - ({m_1}{v_b} - \rho {\rm A}{J_b})\cos \varphi }&{({m_1}{v_n} - \rho {\rm A}{J_n})} \\ 0&0 \\ 0&{1 + eT} \\ {(1 + eT)\cos \varphi }&0 \\ {m{v_t}\cos \varphi + ({m_1}{v_n} - \rho {\rm A}{J_n})\sin \varphi }&0 \\ { - ({m_1}{v_b} - \rho {\rm A}{J_b})\sin \varphi }&{m{v_t}} \end{array}} \right]。\end{split} $

其中： $A = 0.25\;\text{π} {d^2},w = (m - \rho A)g, {m_1} = m + \rho A,e = 1/EA$ 。

1.4 边界条件及求解算法

根据潜水器拖曳、起吊作业需求，可知缆索上端点为速度约束，下端点根据是否拖曳载荷，分为自由边界或平衡约束。由于有限差分法不能处理缆索张力为零的情况，自由边界约束简化为微小张力的平衡约束。

自由边界：缆索末端未系重物，缆索尾部处于无张力状态，末端缆索单元的姿态角沿缆索长度变化率为0。

平衡约束：拖曳或吊放载荷端在拖曳或起吊过程中，缆索下端点受到的外力为拖曳载荷受到的重力、浮力、流体阻力的合力。

速度约束：缆索上端点在局部坐标系下的速度为潜水器航速在缆索局部坐标下的投影与缆索收放速度的合速度。

求解算法：在偏微分方程的求解中，采用具有2阶精度的中心差分格式。收放缆索的模拟采用在每一步迭代更新中，在求解非线性方程组时增加一次只更新缆索微元的长度的求解过程。

2 数值仿真验证 2.1 稳态过程仿真验证

美国海军泰勒舰船研究发展中心在20世纪80年代对水下缆索进行了一系列实船拖曳试验，获得并公开了大量试验数据^[21]。采用本文方法对其直航拖曳试验工况进行仿真，通过仿真数据与试验数据对比，对本文方法进行验证。实船拖曳直航试验的试验缆索分为6节，从下拖点编号1～6，各节缆索物理参数如表1所示。缆索杨氏模量为 $ 2.0\;{\rm{E}}\;10{\text{ Pa}} $ ，流体密度取1024 kg/m³。缆索上设有A和B两个监测点，分别位于第2和第4节缆索末端。

实船拖曳系统以18.5 kn航速直航，缆索上拖点与拖曳系统速度耦合，下拖点为自由端。数值仿真其监测点A和B深度，在该工况下缆索形态如图2所示。仿真得到的监测点A和B的深度与试验结果基本相符（见表2）。

2.2 动态过程仿真验证 2.2.1 验证算例及参数规定

潜水器在拖曳、起吊作业工况下，一般航速较低，舵效差。大型潜水器通常具备浮力调节水箱、辅助推进器等设备，在拖曳和起吊作业下提供较大的浮力和姿态调节能力。避免分析因素过多，影响缆索收放影响的分析，在算例中不考虑其浮力变化或辅助推进器推力变化。

本文算例中采用NPS AUV Ⅱ潜水器模型，该模型只有尾部正反转推进器的推进动力及首、尾的水平和垂直操纵舵面。计算潜水器拖曳一定任务载荷调整航行深度后转直航、转直航后进行边拖曳、边收缆的回收任务载荷过程2个工况。通过上述2种工况的分析和比较对缆索收放过程对潜水器运动及缆索最大张力的影响。

缆索和拖曳载荷的基本参数如表3和表4所示。缆索上端点位置位于潜水器的浮心处，潜水器的模型参数、水动力系数见文献[19]。

表中的参数为默认参数，未明确说明不做改变。

2.2.2 潜水器拖曳载荷变深度转直航工况仿真研究

场景演示：潜水器的缆索下端点与任务载荷相连，释放的缆索长度50 m，任务载荷已完成作业任务，准备进行回收任务载荷作业。为保障任务载荷与海底可能的触碰导致载荷损坏情况的发生，首先进行爬升航行高度流程。潜水器模型见文献[19]。

执行步骤：1）设定航行高度并准备；2）推进装置大功率推进，并操上浮舵；3）达到预设航行高度后，调整进入定深直航工况。

仿真工况：螺旋桨推进转速提高至1500 r/min，航行高度调整100 m，自动操舵，仿真时间1000 s。仿真结果如图3～图5所示。

图3为潜水器航行轨迹与缆索形态曲线变化过程。图4为潜水器航行高度、升降舵舵角、航行速度、机身纵倾角随时间的变化历程图。图5为缆索上端点张力与缆索长度的变化时间历程图。

受力分析：开始阶段（稳态过程）潜水器无速度，缆索近乎垂直，缆索上端点张力约等于拖曳载荷及释放的缆索的重力与浮力之差：

$ {{{T}}^{(1)}}{\text{ = (}}{G_b} - {B_b}{\text{) + }}(m - 0.25\;\rho \;\text{π} {d^2})g \approx 189\; {\text{N}}。$

与图5仿真结果的缆索上端点初始张力一致。通过速度及受力分析可以得出以下结论：

1）在相同推进器转速和操舵策略条件下，由于拖曳体与缆索施加给潜水器拖曳阻力，潜水器需更长时间才能调整到预设安全航行高度，并且由于拖航体和拖缆的流体阻力其稳定后的航行速度明显减小；

2）在潜水器爬升段中，拖曳系统施加在潜水器上的缆索张力先增加后减小，潜水器爬升至预定高度后，升降舵减小、潜水器纵倾角减小，上端点缆索张力随之减小，并在潜水器航速稳定后缆索张力逐渐稳定。该过程中缆索张力的变化与潜水器航行速度在缆索局部坐标系 $t$ 轴的投影长度一致。

2.2.3 变深度转直航后收缆回收工况的仿真研究

场景演示：进行回收任务载荷作业，已完成爬升航行高度流程，执行缆索回收流程。潜水器模型见文献[19]。

执行步骤：1）完成调整航行高度并转入定深直航；2）稳定直航后收缆；3）缆索长度小于阈值后停止收缆。

仿真工况：从上述算例可知在仿真时间400 s左右，潜水器调整至预设航行高度，且速度稳定在0.9 m/s（1.8 kn）左右。可以认为在仿真时间500 s时，其已经达到稳定直航。并开始以0.1 m/s的速度收缆，缆索总长度小于0.5 m时停止收缆。螺旋桨推进转速维持1500 r/min，自动操舵维持航行高度，仿真时间1500 s。仿真结果如图6～图8所示。

图6为拖航收缆工况下潜水器航行轨迹与缆索形态曲线变化过程。图7为仿真过程中潜水器航行高度、升降舵舵角、航行速度、机身纵倾角随时间的变化历程图。图8为缆索上端点张力与缆索长度的变化时间历程图。

受力分析：潜水器停止收缆，达到稳定航速后，其缆索上端点张力为拖曳载荷受到的重力、浮力以及流体阻力的矢量合。潜水器定深直航，流阻与重浮力垂直：

${{{T}}^{(2)}}{\text{ = }}\sqrt {{{{\text{(}}{G_b} - {B_b}{\text{)}}}^2} + {{{\text{(0}}{\text{.5}}\;c_x^ + \rho {L^2}{{\text{u}}^2}{\text{)}}}^2}} \approx 160\; {\text{N}}$ 。与图8仿真所得缆索上端点稳定后张力一致。

通过速度及受力分析可以得出以下结论：

1）在相同推进器转速和操舵策略条件下，在收缆起吊过程中潜水器航行深度和纵倾角均保持在零度附近，并且随着缆索长度减小，拖曳系统施加给潜水器的缆索张力减小，潜水器航速增加。

2）在开始收缆时，拖缆点速度变化导致，缆索张力出现短暂跳变，说明缆索上端点张力对收放缆索的速度变化较敏感，应缓慢改变起吊速度以避免缆索张力剧烈变化。

3 结　语

本文针对潜水器-缆索-载荷三者间相互耦合作用研究中涉及的缆索变长度问题，在仿真办法调研过程中受到集中质量法和有限段法中改变缆索单元数目或凝聚质量点数目来实现缆索变长度问题仿真的启发，基于Ablow推导得到的缆索偏微分控制方程，在非线性方程组迭代时间接更新缆索微元段长度实现缆索收放过程的模拟，并通过稳定拖曳工况与文献对比验证。最后通过拖航变深度航行、拖航过程收缆索进行起吊作业2种算例，对本文方法在水下缆索收放场景模拟中的运用可行性进行了验证，并通过客观分析说明仿真结果的合理性和正确性，相关结论对于开展多用途潜水器拖曳和起吊能力设计具有一定指导意义。主要结论如下：

1）基于Ablow的缆索偏微分控制方程及缆索微元段长度更新策略，可以有效实现缆索收放过程的模拟。

2）潜水器拖曳任务载荷，在相同动力条件下，速度响应与无拖曳工况比较变慢、且稳定航速也降低。潜水器拖曳载荷爬升过程中由于艇体与缆索夹角小，缆索上端点张力较大。

3）在潜水器收缆作业中，缆索张力对收放缆索的速度变化较敏感，应缓慢改变起吊速度以避免缆索张力剧烈变化。同时，随着缆索长度减小，缆索上端点张力也明显减小，说明缆索自身也造成了相当的流体阻力。

本文建立的仿真办法可有效仿真潜水器拖曳航行和起吊航行工况下的运动响应和缆索张力变化，对潜水器水下缆索相关作业的研究提供了合适的仿真手段，但该方法在仿真速度和精度上有待进一步改进。