0 引　言

随着潜器技术的不断发展，潜器水下拖曳系统与相较于传统水面船舶拖曳系统，在海洋工程、海底探测、资源开发以及环境监测等领域中，具有更广泛的应用潜力。潜器水下拖曳系统的螺旋运动是一种较为复杂的三维运动模式，潜器通过控制方向舵和升降舵等舵面的角度，维持螺旋轨迹，空间螺旋运动是一种重要的机动方式，其对潜器的机动性能和稳定性具有关键的指示作用^[1]。研究潜器水下拖曳系统的空间螺旋运动不仅有助于优化潜器的操控性、稳定性和拖曳系统的性能，还为深海作业和精确水下任务提供技术保障。

KENNEDY等^[2]和ABLOW等^[3]通过海上试验和数值模拟研究了包括圆周回转、直线运动和水平正弦运动在内的水面船舶拖曳系统的操纵运动，提供了研究拖曳系统操纵性基础运动模型。WALTON等^[4]、HUANG^[5]和CHAI等^[6]建立、发展了集中质量法，该方法具有明确的物理意义、简便的算法和适应边界条件的能力，适合应用在拖船变速拖曳^[7]、平面回转操纵等拖曳系统操纵运动响应的研究中。KISHORE等^[8]改进了ABLOW^[3]的模型，对水下拖曳线阵列系统圆周回转操纵性进行分析，GROSENBAUGH^[9]重新计算了拖曳系统在船舶平面回转运动时的动态响应，拖曳缆在回转运动中的动态行为可归纳为回转周期和过渡状态的消退的平衡转换。WANG等^[10]对拖曳母船操纵导致的拖曳运动响应进行了详细的总结和比较，提出了多个无量纲系数用以描述平面回转操纵引发的拖曳系统运动响应的操纵性。王浩天等^[11]和周广礼等^[12]为得到潜器的操纵性能，以潜器在水下开展螺旋运动为研究对象，分析高速旋回运动中潜器姿态与轨迹的变化规律。在研究潜器水下拖曳系统的性能时，螺旋运动作为潜器运动的一种复杂模式，对于揭示拖曳系统在复杂水下环境中的操控性和稳定性具有重要意义。

在现有的水下拖曳系统研究中，研究对象主要集中于水面船舶在不同操纵运动下驱动的拖曳系统，已有不少学者提出了无量纲参数来描述水面船舶拖曳系统的操纵性。然而，潜器水下拖曳系统的研究相对较少，尤其是如何利用无量纲参数来分析和优化拖曳系统在复杂运动模式下的性能。本研究采用缆索动力学数值计算工具模拟了水下拖曳系统随潜水器螺旋上浮运动形成的运动响应，通过分析无量纲参数对潜器水下拖曳系统在螺旋上浮运动中的影响，探讨如何优化潜器驱动的拖曳系统的操控性与稳定性，旨在为潜器水下拖曳系统的设计与控制提供新的理论依据和方法。

1 潜器水下拖曳系统的数学模型

为适应潜器拖曳系统水下运动，本文所要建立的拖曳系统的数学模型对水下拖曳系统的相关文献 ^{[13 − 17]}给出的模型进行了改进。包括将潜器的水下六自由度操纵运动以欧拉变换与集中质量法统一。将拖曳体简化为质点模型，该模型主要受重力、浮力和水流阻力的作用。

潜器水下拖曳系统示意图如图1所示，该拖曳系统包含潜器、拖曳缆和拖曳体。将拖曳缆视为柔性缆，忽略其弯曲和扭转刚度的影响。潜器所做运动是一直保持在水面一定深度以下并且为设定静水环境下，拖曳体不上浮到水面，不考虑风、浪、流等因素对拖曳系统的影响。图1中的s为拖曳缆的坐标；L为拖曳缆的总长度；$ i = 0、s = 0 $为起点；$ (\theta ,\phi ) $为拖曳缆的姿态角，其中θ为缆与y轴的顺时针偏离角；$ \phi $为抬升角，$ (\theta ,\phi ) $与缆的空间位置相关。

1.1 系统惯性坐标系与局部坐标系

潜器拖曳时运动为螺旋上浮运动，为了更好地描述拖曳系统的运动，建立了2个坐标系，分别是描述潜器运动的惯性坐标系O – XYZ，描述拖曳缆运动的局部坐标系o - btn。其中，惯性坐标系的原点设在水平面上，即$z = 0$，Z轴选取向上为正，当z值的数据减小时，表示深度增大；反之，表示深度减小。拖曳缆的垂直坐标系中，t轴为拖曳缆的切向，n为法向，b为侧向。拖曳系统的数学模型以图1中的坐标系统为依据，全部的计算都在这2个坐标系下进行。惯性坐标系和拖曳缆的局部坐标系可以通过姿态角的关系进行相互转换，转换关系为：

$ [X,Y,Z]^{\text{T}}=\boldsymbol{\gamma}[b,n,t]^{\text{T}}。$

(1)

其中，${\boldsymbol{\gamma }}$为局部坐标系与惯性坐标系之间的转换矩阵，${\boldsymbol{\gamma }} = \left[ {\begin{array}{*{20}{c}} {{\text{cos}}\theta }&{\sin \theta \cos \phi }&{ - \sin \theta \sin \phi } \\ { - \sin \theta }&{\cos \theta \cos \phi }&{ - \cos \theta \sin \phi } \\ 0&{\sin \phi }&{\cos \phi } \end{array}} \right]$。

1.2 拖曳系统动力学模型建立与求解

在本研究中拖曳缆动力学模型的建立采用集中质量法，从拖曳缆的首端至末端离散为N-1段，得到N+1个节点，其中末端$ s = 0 $处和首端$s = L$处分别为拖曳缆的第$ i = 0 $和第$ i = N $个节点。使用牛顿第二定律在第i个节点上，可得到其基本的运动控制方程为：

$ {{\boldsymbol{M}}_i}{{\boldsymbol{\ddot x}}_i} = \Delta {{\boldsymbol{T}}_i} + {{\boldsymbol{B}}_i} + {{\boldsymbol{G}}_i} + {{\boldsymbol{F}}_{Di}}。$

(2)

式中：$ {{\boldsymbol{\ddot x}}_i} $为加速度；$ {{\boldsymbol{M}}_i} $为质量矩阵；$ \Delta {{\boldsymbol{T}}}_{i}、{{\boldsymbol{B}}}_{i}、{{\boldsymbol{G}}}_{i}、{{\boldsymbol{F}}}_{Di} $为节点i上的所有作用力，分别为拖曳缆张力、浮力、重力以及流体阻力，质量矩阵分为拖曳缆节点上自身的惯性质量与拖曳缆在流体中的附加质量$ {{\boldsymbol{M}}_{a,i}} $, 可分别表示为：

$ {{\boldsymbol{M}}_i} = ({\mu _{i - 1/2}}{l_{i - 1/2}} + {\mu _{i + 1/2}}{l_{i + 1/2}}){\boldsymbol{I}}/2 + {{\boldsymbol{M}}_{a,i}}，$

(3)

$ {{\boldsymbol{M}}_{a,i}} = ({{\boldsymbol{M}}_{a,i - 1/2}} + {{\boldsymbol{M}}_{a,i + 1/2}})/2，$

(4)

$ {{{\boldsymbol{M}}_a} = \rho l\sigma \left[ {\begin{array}{*{20}{c}} {1 - {{\sin }^2}\theta {{\cos }^2}\theta }&{ - \cos \theta \sin \theta {{\cos }^2}\phi }&{ - \sin \theta \cos \phi \sin \phi } \\ { - \cos \theta \sin \theta {{\cos }^2}\phi }&{1 - {{\cos }^2}\theta {{\cos }^2}\phi }&{ - \cos \theta \cos \phi \sin \phi } \\ { - \sin \theta \cos \phi \sin \phi }&{ - \cos \theta \cos \phi \sin \phi }&{{{\cos }^2}\phi } \end{array}} \right]。}$

(5)

式中：I为3×3的单位矩阵；μ、l、σ、ρ分别为拖曳缆单位长度的质量、节点间长度、横截面积、流体密度；节点$ \mathrm{i} $i与i+1之间的量：i+1/2为下标。

1）张力T，拖曳缆的应变在通常情况下为$\left| \varepsilon \right| \ll 1$，应力—应变关系做弹性化处理，使用胡克定律描述其本构关系，计算式为：

$ \Delta {{\boldsymbol{T}}_i} = {{\boldsymbol{T}}_{i + 1/2}} - {{\boldsymbol{T}}_{i - 1/2}}，$

(6)

$ {{\boldsymbol{T}}_{i + 1/2}} = E\sigma {\varepsilon _{i + 1/2}}{{\mathbf{\tau }}_{i + 1/2}} ，$

(7)

$ \varepsilon_{i+1/2} = \frac{\sqrt{(x_{i+1} - x_i)^2 + (y_{i+1} - y_i)^2 + (z_{i+1} - z_i)^2}}{l_{i+1/2}} - 1 。$

(8)

式中：τ为拖曳缆缆长方向上的单位切向量；E为拖曳缆的弹性模量。

2）浮力和重力B、G，拖曳缆节点i上所受到的浮力和重力可以表示为：

$ \begin{split} {{\boldsymbol{B}}_i} + {{\boldsymbol{G}}_i} =& - \rho {\boldsymbol{g}}({l_{i - 1/2}}{\sigma _{i - 1/2}} + {l_{i + 1/2}}{\sigma _{i + 1/2}})/2 + \\ &\rho {\boldsymbol{g}}({l_{i - 1/2}}{\mu _{i - 1/2}} + {l_{i + 1/2}}{\mu _{i + 1/2}})/2 。\end{split}$

(9)

式中：g为重力加速度。

3）流体阻力F_Di，将拖曳缆阻力由切向阻力和法向阻力组成，计算式为：

$ {{\boldsymbol{F}}_{Di}} = ({{\boldsymbol{F}}_{Di + 1/2}} + {{\boldsymbol{F}}_{Di - 1/2}})/2，$

(10)

$ {{\boldsymbol{F}}_{Di}} \approx - \rho ld\sqrt {1 + \varepsilon } ({C_n}\left| {{{{\boldsymbol{\dot x}}}_{n,b}}} \right|{{\boldsymbol{\dot x}}_{n,b}} + \pi {C_t}\left| {{{{\boldsymbol{\dot x}}}_t}} \right|{{\boldsymbol{\dot x}}_t}) 。$

(11)

式中：${\boldsymbol{\dot x}}$为节点的速度向量；d为缆直径；${C_n}$、${C_t}$别为其法向、切向阻力系数。

本研究中拖曳缆末端节点$ i = 0 $处连接拖曳体，并将拖曳体简化为质点模型，该模型主要受重力、浮力和水流阻力的作用，此处的边界条件为：

$ {{\boldsymbol{M}}_b}\frac{{{d^2}{{\mathbf{{\rm A}}}_s}}}{{{\rm{d}}{t^2}}} = {\boldsymbol{S}}。$

(12)

质量系数矩阵的非零元素为：

$ {{\boldsymbol{M}}_{b 1,1}} = m + {m_{ax}}\text{，} $

(13)

$ {{\boldsymbol{M}}_{b 2,2}} = m + {m_{ay}} \text{，} $

(14)

$ {{\boldsymbol{M}}_{b 3,3}} = m + {m_{az}}\text{，} $

(15)

$ {{\boldsymbol{A}}_S} = {\left[ {{X_b} {Y_b} {Z_b}} \right]^{\text{T}}} \text{，} $

(16)

$ {\boldsymbol{S}} = {[{X_F} {Y_F} {Z_F}]^{\text{T}}}\text{，} $

(17)

$ {X_F} = \rho {U^2}{L_x}^2/2\text{，} $

(18)

$ {Y_F} = \rho {U^2}{L_y}^2/2\text{，} $

(19)

$ {Z_F} = \rho {U^2}{L_z}^2/2 + mg - \rho gV \text{，} $

(20)

式中：m为拖曳体的质量；m_ax、m_ay、m_az均为拖曳体的附加质量；X_F、Y_F、Z_F分别为作用在拖曳体上的重浮力、惯性力、水流作用力、拖曳缆拉力等载荷；L_x、L_y、L_z为计算水流作用力所采用拖曳体的特征长度。

由于潜器与拖曳缆连接处节点$i = N$的位移和速度均与潜器的位移和速度相同，因此此处的边界条件为：

$ \left\{ \begin{gathered} {x_N} = {x_s}(t)，\\ {y_N} = {y_s}(t)，\\ {z_N} = {z_s}(t)。\\ \end{gathered} \right. $

(21)

$ \left\{ \begin{gathered} {{\dot x}_N} = {u_s}(t)，\\ {{\dot y}_N} = {v_s}(t)，\\ {{\dot z}_N} = {w_s}(t) 。\\ \end{gathered} \right. $

(22)

式中：N为节点$i = N$；x_s、y_s、z_s和u_s、v_s、w_s分别为潜器与拖曳缆连接处节点关于时间的坐标分量和速度分量。

将式(2)和式(6)联立可以拖曳缆的运动控制方程，即：

$ v = {\rm{d}}x/{\rm{d}}t，$

(23)

$ \left\{ \begin{gathered} \frac{{{\rm{d}}{{\dot x}_i}}}{{{\rm{d}}t}} = {M_i}^{ - 1} \cdot {F_i} \\ \frac{{{\rm{d}}{x_i}}}{{{\rm{d}}t}} = {{\dot x}_i} \\ \end{gathered} \right. (i = 0, \cdot \cdot \cdot ,N - 1)。$

(24)

通过将式(7)与相应的边界条件联立，并在给定初值条件下，应用四阶龙格-库塔方法进行数值积分，便可求解该方程。

2 水下拖曳系统参数与工况

潜器水下拖曳系统由提供驱动的潜器、作为拖曳系统主体的拖曳缆和执行具体任务的拖曳体3部分组成，每个组成部分的结构参数对系统的动态行为、稳定性和操控性能有着直接影响。为了统一研究的标准和建立对应的无量纲参数，本研究参考WANG等 ^[10]在研究中使用的拖曳缆和拖曳体的结构参数，如表1所示，水下拖曳系统中主要的结构参数。潜器在5000 s内从水深500 m的静水环境中，向上开展螺旋上浮运动，水平回转速度V_ζ和垂向潜浮速度V_t分别为4 kn和0.1 kn。本研究参考文献[10]和[17]建立关于潜器螺旋上浮运动中水下拖曳系统的无量纲参数，通过将各参数无量纲化，可以提高结论的适用性，进而揭示不同参数对系统响应的普遍影响，该系统的基准无量纲参数如表2所示。

如图2所示，潜器在螺旋上浮运动中，每隔100 s记录一次拖曳系统形态变化，从0 s开始，一直持续到5000 s，期间拖曳体的运动明显滞后于潜器的运动。图3为潜器螺旋上浮运动期间拖曳系统的运动响应，分别为拖曳缆与潜器连接处张力的变化历程和拖曳体的升沉变化历程。张力和拖曳体升沉的变化历程基本可分为瞬态阶段和稳态阶段，历程的整体变化趋势都呈现为先增大后减小再增大。在稳态阶段，张力趋于稳定并维持在一个固定值，而拖曳体的深度则表现为线性增长。

3 无量纲参数对运动响应的影响 3.1 潜器回转半径与拖曳缆总长度之比

本研究中潜器回转半径R与拖曳缆总长度L之比R/L的基准取值为0.40，反映了潜器的空间回转特性与拖曳系统整体尺寸之间的关系，通过调整潜器的回转半径R 分别为25、50、125和500 m，使R/L分别为0.05、0.10、0.25和1.00。如图4所示，潜器回转半径R与拖曳缆总长度L之比R/L对系统运动响应的影响。

对于张力变化历程，当R/L＜0.40时，瞬态阶段的张力出现明显的波动，且R/L越小，波动的频率越高；当R/L ≥ 0.40 时，瞬态阶段的张力波动现象消失，同时R/L越大，瞬态阶段持续的时间越短，稳态阶段的张力值越低。对于拖曳体的升沉变化历程，类似地，当R/L＜0.40时，瞬态阶段的拖曳体深度出现明显波动，且波动频率随着R/L的增大而降低，整体变化趋势呈现为先减小再增大；当R/L ≥ 0.40 时，瞬态阶段的波动现象消失。同时，随着R/L的增大，瞬态阶段持续时间缩短，而仿真结束时拖曳体的最终深度值逐渐减小。

3.2 拖曳缆总质量与拖曳体质量之比

本研究中拖曳缆总质量与拖曳体质量之比ω的基准取值为3.19，该参数描述了拖曳缆与拖曳体在质量分布上的相对权重对系统动态行为的影响，通过调整拖曳体质量 分别为751.0、25.0、18.8 kg，使ω分别为0.50、15.00、20.0。如图5所示，拖曳缆总质量与拖曳体质量之比ω对系统运动响应的影响。

对于张力变化历程，整体的变化趋势为先增大后减小再增大，ω的值越大则稳态阶段的张力值越小；当ω = 0.50时，整个历程中张力显著高于ω为其余值。对于拖曳体的升沉变化历程，仿真结束时拖曳体的最终深度值随ω增大而增大；当ω ≥ 15时，整体的变化趋势为先减小再增大；当ω < 15时，变化趋势为先增大后减小再增大。ω的大小对变化历程中从瞬态阶段过渡到稳态阶段的时间影响不显著。

3.3 拖曳缆单位长度质量与阻力之比

本研究中拖曳缆单位长度质量w与单位长度阻力r之比w/r的基准取值为0.066，该参数反映了拖曳缆质量与水动力阻力之间的平衡关系，通过调整拖曳缆法向阻力系数分别为2.644、1.312、0.313，使w/r分别为0.050、0.100、0.400。如图6所示，拖曳缆单位长度质量与单位长度阻力之比w/r对系统运动响应的影响。

对于张力的变化历程，其整体的变化趋势表现为先增大后减小再增大，随着 w/r的增大，瞬态阶段持续的时间逐渐缩短，而稳态阶段的张力值则随之增大。对于拖曳体的升沉变化历程，同样呈现出先增大后减小再增大的变化趋势，随着w/r的增大，瞬态阶段的时间缩短，但拖曳体的最终深度值逐渐减小。同时，w/r的大小对运动响应整体变化趋势的影响不显著。

3.4 水平回转速度和垂向潜浮速度之比

本研究中潜器水平回转速度V_ζ和垂向潜浮速度V_t之比V_ζ/V_t的基准取值为40，该参数反映了潜器在螺旋运动中的速度构型，通过调整潜器水平回转速度分别为2、6、8 kn，使V_ζ/V_t分别为20、60和80。如图7所示，潜器水平回转速度和垂向潜浮速度之比V_ζ/V_t对系统运动响应的影响。

对于张力的变化历程，当V_ζ/V_t为60和80时，其整体的变化趋势表现为先增大后减小；当V_ζ/V_t为40时，变化趋势表现为先增大后减小再增大；当V_ζ/V_t为20时，变化趋势则表现为多次波动。当V_ζ/V_t ≥ 40时，V_ζ/V_t的值越大则稳态阶段的张力值越大。对于拖曳体的升沉变化历程，当V_ζ/V_t为60和80时，其整体的变化趋势表现为瞬态阶段的陡增和稳态阶段较缓的线性增加；当V_ζ/V_t为20和40时，变化趋势表现为先增大后减小再增大；V_ζ/V_t的值越大则拖曳体的最终深度值越大。对于从瞬态阶段过渡到稳态阶段的时间来说，V_ζ/V_t的值越大则瞬态阶段的时间越短。

3.5 拖曳缆浮力与重力之比

本研究中拖曳缆浮力B_c与重力G_c之比B_c/G_c的基准取值为0.319，该参数反映了拖曳缆在水中受力的平衡状态，通过调整拖曳缆单位长度质量分别为2.400、0.480、0.240、0.048 kg/m，使B_c/G_c分别为0.100、0.500、1.000、5.000。如图8所示，拖曳缆浮力与重力之比B_c/G_c对系统运动响应的影响。

对于张力的变化历程，当B_c/G_c ≥ 0.319时，其整体的变化趋势表现为先增大后减小再增大；当B_c/G_c = 0.100时，变化趋势表现为先减小后增大再减小。当B_c/G_c = 1.000时，稳态阶段的张力处于最小值，而B_c/G_c = 0.100时，整个运动响应过程中张力值均明显高于其余参数的张力值。对于拖曳体的升沉变化历程，整体的变化趋势均表现为先增大后减小再增大，B_c/G_c的值越小则拖曳体的最终深度值越小；当B_c/G_c = 0.319时，瞬态阶段的持续时间显著长于其他B_c/G_c值的情况。

4 结　语

本文研究了潜器在螺旋上浮运动中的水下拖曳系统动力学特性，构建了适用于该系统的数学模型，并通过引入无量纲参数量化其动态响应规律。研究分析了无量纲参数——潜器回转半径与拖曳缆总长度之比R/L、拖曳缆总质量与拖曳体质量之比ω、拖曳缆单位长度质量与单位长度阻力之比w/r、潜器水平回转速度与垂向潜浮速度之比V_ζ/V_t以及拖曳缆浮力与重力之比B_c/G_c对运动响应的影响，其具体表现如下：

1）当R/L增大时，瞬态阶段的张力波动现象逐渐减弱，瞬态阶段的持续时间明显缩短。同时，稳态阶段的张力值逐渐降低。对于拖曳体的升沉变化，当R/L小于临界值时，拖曳体深度在瞬态阶段出现明显波动，而R/L 增大后，这种波动消失，最终拖曳体的深度值随着R/L的增大而减小。

2）随着V_ζ/V_t的增大，张力的变化历程表现出从多次波动逐渐转变为较为平滑的先增大后减小再增大的趋势，同时稳态阶段的张力值逐渐升高。对于拖曳体升沉的变化，当V_ζ/V_t较小时，其变化趋势表现为明显的多次波动，而随着V_ζ/V_t的增大，拖曳体深度变化趋势逐渐趋于平稳，最终深度值随 V_ζ/V_t增大而升高。

3）随着w/r的增大，瞬态阶段的持续时间逐渐缩短，张力在稳态阶段趋于更高的稳定值。而对于拖曳体的升沉变化，整体趋势表现为先增大后减小再增大，最终拖曳体深度值则随w/r的增大而逐渐减小。

4）当ω增大时，稳态阶段的张力值逐渐减小，而拖曳体深度值在仿真结束时呈现逐渐增加的趋势。对于B_c/G_c，当B_c/G_c = 1.000时，拖曳缆浮力与重力达到平衡状态，此时系统表现出显著的特性。对于张力变化历程，稳态阶段的张力值达到最低，瞬态阶段的波动较小且持续时间较短，拖曳系统整体运动响应较为平稳。

本研究通过建立潜器螺旋上浮运动下的拖曳系统动力学模型，系统揭示了无量纲参数对系统运动响应的调控机制，这些规律为潜器拖曳系统在复杂三维运动中的参数匹配策略提供了直接依据，并对水下拖曳系统的设计与优化提供了理论支撑。