0 引　言

水下潜器是一种具有水下观察和作业能力的海洋工程装置，被广泛应用于观测海洋环境、开发海洋资源、检查和维护水下设施等任务。水下潜器通常由潜器主体以及作为其运动控制机构的导管螺旋桨组成，主体内部可以根据不同的水下作业任务搭载不同类型的参数监测传感器和水下作业装备，水面工作人员可以通过脐带缆对控制机构发出指令来操纵潜器的轨迹与姿态^[1-3]。在海洋研究与开发需求不断增大的背景下，水下潜器作为一种在海洋监测与资源开发等领域不可或缺的装置，受到越来越多的关注。

水下潜器在作业时，由于作业任务不同、水域复杂环境等因素的影响，其运动形式存在多种变化^[4]。利用数值手段准确地模拟出水下潜器在复杂运动下的水动力特性，可以为设计出控制性能优良的水下潜器提供前提。

现阶段国内外学者对水下潜器在各种运动形式下的水动力特性研究中，主要以一维直线和二维平面运动为主。Jason等^[5]运用CFD方法设计一种针对C-SCOUT型自主式水下潜器的动力学模型，并在加入了通体推进器的条件下，利用相对运动原理，确定了360°攻角范围内作用在水下潜器上的外力。金晓东^[6]将缆绳、潜器主体、导管螺旋桨耦合，建立一个三维水动力数学模型，使用CFD软件Fluent模拟了整个拖曳系统的运动过程，但运动过程限制在二维平面上。Bellingham等^[7]通过采用势流分析和三维粘性CFD建模设计了一个低阻力的Tethys型自主式水下潜器，探讨了螺旋桨与主体的相互作用问题，给出了功耗和效率的初步试验结果。胡维峰^[8]运用动网格技术，数值模拟了水下潜器在自航中遭遇海流干扰时通过偏转舵保持其航向的过程，验证了转舵操纵数值模拟的可行性与可靠性，运动过程同样也仅限于二维平面。文献[9]中利用动网格与滑移网格结合的方法对水下潜器二维回转运动状态的整个过程进行模拟，这种方法虽然能够模拟水下潜器实际的运动过程，但网格质量会随每一次网格更新而下降，特别是在高速旋转螺旋桨的附近流场等大变形区域，如果计算时间步长选取不当，容易产生负体积网格导致计算失败。上述文献的研究大多采用传统的动网格方法，且运动形式仅限于单个或2个维度之中的简单运动，未对三维立体空间中复杂运动形式下水下潜器的水动力特性进行研究。

重叠网格法是近年来逐渐兴起的解决船舶领域CFD模拟问题的一种方法。重叠网格法相比传统动网格法，其优势在于能实现与之相同的功能而不会出现负体积网格，在模拟多体相对运动方面具有很好的普适性。鉴于重叠网格法的以上优点，采用多重重叠网格技术在由“潜器主体+导管螺旋桨”组成的水下潜器系统进行三维回转运动的流场内求解其RANS方程，观察水下潜器系统在三维回转运动中的动力特征，分析导管螺旋桨在三维回转运动过程中的推力特性。计算中，采用用户自定义函数(UDF)的方式定义水下潜器主体回转运动方式，通过定义叠加运动对各导管螺旋桨的旋转运动进行指定，程序中运用了uclib.h头文件，进而实现所需要的数值模拟。

1 计算方法 1.1 控制方程与湍流模型

设定流体是不可压的，则水下潜器在运动过程中的连续性方程和动量方程分别为^[10]：

$ \frac{\partial {U}_{i}}{\partial {x}_{i}}=0 ，$

(1)

$ \frac{{\partial {U_i}}}{{\partial t}} + {U_j}\frac{{\partial {U_i}}}{{\partial {x_j}}} = {f_i} - \frac{1}{\rho }\frac{{\partial p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left( {\nu \frac{{\partial {U_i}}}{{\partial {x_j}}} - \overline {{u_i}{u_j}} } \right)，$

(2)

式中： $ {U}_{i} $ 和 ${U}_{j}\left(i,j=1,2,3\right) $ 为时均速度分量， $ p $ 为压强， $ {f_i} $ 为质量力分量， $\ \rho $ 为流体密度， $ \nu $ 为流体动力粘度， $ - \overline {{u_i}{u_j}} $ 为雷诺应力。由于雷诺应力为未知量，此时方程中的未知量数大于方程数。使用SST $ k - \omega $ 湍流模型来封闭方程组，其方程如下^[11]：

$ \frac{{\partial (\rho k)}}{{\partial t}} + \frac{{\partial (\rho k{u_i})}}{{\partial {U_i}}} = \frac{\partial }{{\partial {U_j}}}[{\varGamma _k}\frac{{\partial k}}{{\partial {U_j}}}] + {G_k} - {Y_k} + {S_k}，$

(3)

$ \frac{{\partial (\rho \omega )}}{{\partial t}} + \frac{{\partial (\rho \omega {u_i})}}{{\partial {U_i}}} = \frac{\partial }{{\partial {U_j}}}[{\varGamma _\omega }\frac{{\partial \omega }}{{\partial {U_j}}}] + {G_\omega } - {Y_\omega } + {S_\omega }，$

(4)

式中： $ k $ 为湍动能； $ \omega $ 为湍流耗散率； $ {\varGamma _k} $ 和 $ {\varGamma _\omega } $ 分别为 $ k $ 和 $ \omega $ 的有效扩散项； $ {G_k} $ 为由于平均速度梯度导致的湍动能 $ k $ 的产生项， $ {G_\omega } $ 为 $ \omega $ 的产生项； $ {Y_k} $ 和 $ {Y_\omega } $ 分别为 $ k $ 和 $ \omega $ 的耗散项； $ {S_k} $ 和 $ {S_\omega } $ 为自定义源项。有效扩散项定义为：

$ {\varGamma _k} = \mu + \frac{{{\mu _t}}}{{{\sigma _k}}}，$

(5)

$ {\varGamma _\omega } = \mu + \frac{{{\mu _t}}}{{{\sigma _\omega }}}，$

(6)

$ {\mu _t} = \frac{{\rho k}}{\omega }\frac{1}{{\max \left[ {1,\dfrac{{S{F_2}}}{{{a_1}\omega }}} \right]}}。$

(7)

式中： $\mu $ 为流体运动粘度； ${\mu _t}$ 为湍流粘度比； ${\sigma _k}$ 和 ${\sigma _\omega }$ 分别为 $ k $ 和 $ \omega $ 的湍流普朗特数； $ S $ 为应变率； ${F_2} = \tanh $ $ \left( {\varPhi _2^2} \right)$ ； ${\varPhi _2} = \max \left[ {2\dfrac{{\sqrt k }}{{0.09\omega y}},\dfrac{{500\mu }}{{\rho {y^2}\omega }}} \right]$ ； $ y $ 为到下一表面的距离； $ {a_1} $ = 0.31。

1.2 几何模型

本文所研究的水下潜器系统由 “水下潜器主体+3个标准Ka4-70/19A导管螺旋桨”所组成。水下潜器主体为立式对称倒三角型壳体结构，潜器主体关于中纵剖面对称，结构两侧有2个鱼雷状浮体，浮体轴线相互平行，浮体之间通过水平翼板连接，2个浮体尾部各设置1部标准Ka4-70/19A导管螺旋桨。主体下部为一纵剖面形状为Myring型的舱室，舱室后设置1部标准Ka4-70/19A导管螺旋桨，该桨的桨毂中心线位于中纵剖面上。在潜器系统中所设置的3部导管螺旋桨的几何参数一致。规定从潜器尾部向前看位于潜器左侧的桨为左桨，位于右侧的桨为右桨，位于Myring型舱室尾部的桨为下桨。潜器系统几何模型如图1所示，图中同时给出了计算中所使用的直角坐标系统，坐标原点位于立柱前缘距离舱室上表面126.55 mm处，立柱前缘位于潜器主体中纵剖面上。表1和表2给出了水下潜器主体和导管螺旋桨的主要几何参数。

1.3 计算域及边界条件

为了分析图1所示水下潜器系统在三维回转运动中的水动力特征，构造了图2计算域I～计算域V这5个相互耦合的计算域。

计算域I为一个长度方向与来流速度平行的长方体计算域，为避免边际效应对水下潜器的运动过程产生影响，计算域I的尺度需要足够大。基于这样的考虑，计算中长方体的尺度取为 $\left(\right. {\text{34}}{\text{.0}}L \times 28.5D \times 22.5H \left.\right)$ ，这一计算域作为水下潜器运动的背景网格。

计算域II为水下潜器及导管螺旋桨附近的矩形流场计算域，长度方向与水下潜器纵中剖面方向平行。这一计算域随着水下潜器作回转运动而旋转，用于实现水下潜器三维回转运动。计算域II不包括3个导管螺旋桨周围的流场计算域III、计算域IV和计算域V。综合考虑数值计算精度及计算量，本文计算域II的尺度为 $\left(\right. {\text{2}}{\text{.0}}L \times 1.5D \times 1.5H \left.\right)$ 。

计算域III、计算域IV以及计算域V均为包含导管螺旋桨在内尺度相同的3个圆柱形计算域。计算域III为左桨计算域，计算域IV为右桨计算域，计算域V为下桨计算域。这3个计算域用于实现导管螺旋的三维回转运动以及螺旋桨的旋转。每一个圆柱形计算域的轴线与螺旋桨桨毂中心线重合。计算域的直径取150%导管入口直径，长度取150%导管长度。这样的取值是为了保证重叠网格在划分时区域边界距离壁面至少有4层网格单元。

计算域I与计算域II之间通过重叠网格边界条件进行耦合，其中计算域I作为计算域II的背景网格；计算域II与计算域III、计算域IV以及计算域V之间通过重叠网格边界条件进行耦合，而计算域II又作为计算域III、计算域IV以及计算域V的背景网格。通过这样的重叠网格嵌套，实现多重重叠网格的计算。计算域III、计算域IV以及计算域V之间互相独立，计算域I速度入口为流速V=0.2 kn的无空泡均匀来流。各边界条件定义见表3。

1.4 网格划分

使用STAR-CCM+切割体网格生成器对每个计算域单独划分网格，每个网格单元均为六面体。由于外围计算域与水下潜器的尺度差异较大，需要对计算域进行逐级加密，加密区域（见图2）集中在水下潜器预定的运动范围内。这一加密操作不仅能最大限度减少网格数量、降低计算量，还能辅助检验运动模型的正确性。为最大程度保证重叠网格的插值精度，在每个重叠网格边界两侧使用相同密度数量级的网格。在螺旋桨的导边、随边、叶梢、叶根以及导管边缘处曲率变化较大，需要额外设置网格加密。本文对不同位置的3个导管螺旋桨使用完全相同的网格尺寸控制。图3为潜器主体附近网格分布图，图4为导管螺旋桨附近网格细节。图3与图4中均可见重叠网格边界。

1.5 水下潜器及导管螺旋桨的运动方式

利用构造的水下潜器主体及导管螺旋桨几何模型和计算域对水下潜器系统以回转角速度 $ \omega $ 沿一回转半径的圆周轨迹运动，同时水下潜器主体上附属的3个导管中螺旋桨（见图1）以一定的转速旋转时所表现出来的水动力特性进行分析，计算水下潜器系统在沿水平面作圆周轨迹运动（计算工况Ⅰ）和沿与水平面有一夹角的斜面（计算工况Ⅱ）这2种计算工况下作圆周轨迹运动时，3个导管中螺旋桨所发出来的推力特性以及水下潜器主体所受的水动力载荷。其中，计算工况Ⅰ可以看做是计算工况Ⅱ在 $ \varphi = 0 $ 的一种特例；计算工况Ⅰ所描述的水下潜器系统回转运动本质上属于一种二维问题，而计算工况Ⅱ则是一种三维问题。

计算中，将水下潜器及导管螺旋桨这一整体的初始运动中心位置选定为坐标原点，水下潜器以及导管螺旋桨进行三维回转运动的示意图如图5所示，图中同时给出了本文所使用的直角坐标系。

可知，水下潜器及导管螺旋桨绕一固定轴 $l$ 进行角速度恒定为 $ \omega $ 的逆时针圆周运动，回转中心为P点，A点为水下潜器及导管螺旋桨的初始运动中心。定义水下潜器及导管螺旋桨从初始位置开始绕回转中心P运动过的角度为回转角 $\theta $ （图5中的 $ \angle AP{A_1} $ ），B点为回转角 $\theta = \text{ π} $ 时的运动中心。定义A点为运动最高点，B点为运动最低点，回转半径 $R = PA = PB$ 。定义升沉角 $\varphi $ 为 $AB$ 与水平面的夹角，由于所使用的坐标系中xOy平面与水平面平行，因此 $\varphi $ 等于 $AB$ 与y轴的夹角。采用定义UDF来控制水下潜器及导管螺旋桨的运动方式，回转中心P、回转轴l以及B点的位置不需要人为定义，只需定义 $ x，y，z $ 方向运动速度、回转角速度 $ \omega $ 以及升沉角 $\varphi $ 。

UDF是一段C语言程序，通过STAR-CCM+的UserFunctions.lib进行编译，将需要实现的函数导入STAR-CCM+中。本文所使用的UDF中关于水下潜器及导管螺旋桨运动的主要代码如下：

$ {\rm Velocity}\left[ i \right]\left[ 0 \right]{\text{ }} = {\text{ }} - \frac{{r\omega }}{2}\sin \omega t，$

(8)

$ {\rm Velocity}\left[ i \right]\left[ 1 \right]{\text{ }} = {\text{ }} - \frac{{r\omega }}{2}\cos \omega t，$

(9)

$ {\rm Velocity}\left[ i \right]\left[ 2 \right]{\text{ }} = {\text{ }}\frac{{2r\omega \cdot temp}}{\text{π} }\sin \varphi。$

(10)

其中：i为时间步；Velocity[i][0]，Velocity[i][1]，Velocity[i][2]分别为x，y，z方向速度；t为运动时间；temp的取值规则为：当 $0 \leqslant \theta < \text{π} $ 取1，当 $\text{π} \leqslant \theta < 2\text{π} $ 取–1。

使用STAR-CCM+内置工具中的叠加运动功能对螺旋桨的旋转进行定义。在导管螺旋桨随潜器主体一起进行回转运动已经由UDF定义的前提下，只需计算出螺旋桨的旋转轴原点（取螺旋桨的盘面中心作为旋转轴原点）并指定螺旋桨旋转速率。

2 数值模拟结果与分析

对在计算工况Ⅰ和计算工况Ⅱ这2种水下潜器系统回转运动中水动力特性进行计算分析。计算中，水下潜器系统回转半径R、回转角速度 $ \omega $ 以及导管中螺旋桨的转速N分别取为： $R = 3.0\;{\rm{m}}$ ， $\omega = {\text{π} }/{{15}}\; {{\rm rad}}/{\rm s}$ ， ${N} = 900.0\;{\rm{r/min}}$ 。对于计算工况Ⅰ， $ \varphi = 0 $ ；对于计算工况Ⅱ， $\varphi = \text{π} /15$ 。

在2种计算工况的回转运动数值分析中，所使用计算域划分与计算域边界条件均相同。在网格划分上，环形加密区域的位置随升沉角改变，但使用的网格密度相同。2种运动情况的网格数量分别为423.39万和419.49万。从运动最高点出发回到最高点为一个完整的圆周运动，用时为一个周期。图6给出了 $\varphi = \text{π} /15$ 时水下潜器在初始时刻及完成一个运动周期并继续运动时的计算域网格变化示意图。

从图6计算网格的变化情况可以看出，水下潜器及导管螺旋桨完成了一个周期（回转角 $\theta = 2\text{π} $ ）的三维回转运动过程。由于使用了多重重叠网格，在运动过程中，网格质量一直保持在与初始时刻相同的较高水平，在每个重叠网格边界两侧的网格拥有相同密度的数量级。数值模拟结果证实了提出的运动模型能够实现水下潜器及导管螺旋桨的三维回转运动。

图7～图9给出了水下潜器3个导管螺旋桨在 $\varphi = 0$ 及 $\varphi = \text{π} /15$ 运动情况下，完成2个回转周期的过程中发出的推力随时间的变化曲线。可以看到，在2种运动情况下，各导管螺旋桨发出的推力最终都趋于稳定的周期性变化，且变化周期与水下潜器进行圆周运动的周期相同。这是由于在回转过程中，导管螺旋桨轴向与来流V的方向存在一个夹角，由图5(a)可知这一夹角在数值上等于回转角 $\theta $ ，而各导管螺旋桨轴向来流分量 $V\cos \theta $ 的变化将导致轴向推力的变化，因此在回转角速度 $\omega $ 恒定这一条件下，导管螺旋桨轴向推力呈现与 $\theta $ 同周期的变化。

分别观察图7～图9可以发现，在同一时刻下， $\varphi = $ $ \text{π} /15$ 情况下各桨发出的推力值要小于 $\varphi = 0$ 情况下各桨发出的推力值。事实上，由式(10)可知，当 $\varphi = 0$ 时各导管螺旋桨在Z方向上的速度分量为0，而 $\varphi = \text{π} /15$ 时各桨在Z方向上的速度分量不为0，因此当 $\varphi = \text{π} /15$ 时，这就导致各桨的进速在同一时刻要大于 $\varphi = 0$ 时的进速，进而使桨叶切面的水动力螺距角增大，攻角减小，产生的推力值较 $\varphi = 0$ 时小。

图10给出了水下潜器主体在 $\varphi = 0$ 及 $\varphi = \text{π} /15$ 平面运动情况下Z方向上的流体阻力随时间变化曲线。可以发现，水下潜器主体所受到的流体阻力在Z方向上的分量在一个回转周期内呈先减小后增大的变化趋势，但对比图10与图7～图9可知，在水下潜器三维回转运动过程中，水下潜器所受流体阻力的变化会出现剧烈的波动，这与导管螺旋桨推力的稳定变化规律有较大差异。主要原因在于水下潜器在进行三维回转运动过程中存在一个速度的升沉运动分量（即式（10）），这一升沉运动分量使得水下潜器主体下表面对水流的流动存在阻碍效应，阻碍效应加剧了主体周围流场的复杂程度，导致流体阻力的波动剧烈，与从水中提起一块平板时的感受到水流的不规律扰动相符，这一点在图11中有所体现。图11给出了 $\varphi = \text{π} /15$ 、回转角 $\theta = \text{π} /6$ 时水下潜器周围流线，从图中可以看出水平翼板下表面以及舱室下表面对水流明显的阻碍作用，有一部分水流在流经水平翼板下表面以及舱室下表面时，由于受到非流线型表面的阻挡，导致水流流向发生大幅度的改变。这种阻碍效应还会加剧后方导管螺旋桨所在区域流场的复杂程度，导致螺旋桨发出的推力产生明显的波动，这与图9中 $\varphi = \text{π} /15$ 条件下的下导管螺旋桨推力曲线相符。

上述现象还可以从图10中分析得到。在三维回转运动( $\varphi = \text{π} /15$ )条件下，水下潜器主体所受Z方向阻力值要大于二维回转运动( $\varphi = 0$ )条件下的阻力值，且阻力曲线的波动相较二维回转运动( $\varphi = 0$ )条件下更剧烈。在二维回转运动( $\varphi = 0$ )条件下，水下潜器不具有垂向的运动速度分量，水平翼板下表面以及舱室下表面始终与来流方向平行，对水流的阻碍作用不明显，因此主体Z方向所受阻力在一个运动周期内大致呈现先减小后增大的规律性变化趋势。但在三维回转运动( $\varphi = \text{π} /15$ )条件下，由于水下潜器的垂向的运动速度分量的存在，水下潜器在回转运动过程中，非流线型表面不再与来流方向保持平行（见图11），对来流存在阻碍作用，使得主体周围流场的复杂程度增加，所受阻力存在明显的波动。这一点还可以从图12与图13中观察比较得到。图12与图13分别给出了回转角 $\theta = $ $ \text{π} /6$ 时 $\varphi = \text{π} /15$ 及 $\varphi = 0$ 运动情况下水下潜器主体表面压力在同一尺度下的分布云图。对比2种运动情况下水下潜器主体上下表面压力可知， $\varphi = \text{π} /15$ 时的压力差显然要大于 $\varphi = 0$ 时，因此水下潜器主体所受Z方向阻力前者要大于后者，这与图10的计算结果一致。

3 结　语

本文利用多重重叠网格方法对水下潜器系统进行了二维与三维回转运动的数值模拟，对在无潜器影响和存在潜器影响下不同位置的导管螺旋桨的推力特性进行了计算，重点对水下潜器系统沿斜面作三维圆周轨迹运动时所表现出来的导管螺旋桨推进力以及潜器主体所受到的水动力特点进行计算观察。数值计算结果表明，提出的三维回转运动的数学模型能够正确地模拟水下潜器及导管螺旋桨在水下进行三维回转运动的完整过程。

通过分析数值模拟结果，有以下结论：

1）利用多重重叠网格法，水下潜器及导管螺旋桨可以实现指定的三维回转运动水动力特性计算分析，计算域网格质量随运动时间的变化一直保持在初始运动时刻的较高水平，可以验证提出的三维回转运动数学模型的准确性。

2）在水下潜器及导管螺旋桨进行单向性的三维回转运动过程时，各导管螺旋桨所产生的推力均呈周期性变化。在每个运动周期内，各导管螺旋桨的推力特性随回转角度的改变发生变化，原因在于各导管螺旋桨轴向来流分量 $V\cos \theta $ 的变化将导致螺旋桨所感受到的进速不同，使得螺旋桨发出不同的推力值。

3）相比于二维回转运动( $\varphi = 0$ )，水下潜器在进行三维回转运动过程中，由于存在一个垂向的运动速度分量，导致水下潜器所表现出来的导管螺旋桨推力特性及阻力特性与二维回转运动下有较大差异。三维回转运动条件下各桨发出的推力值较二维回转运动条件下小，潜器主体所受Z方向的流体阻力较二维回转运动条件下更大且波动更加剧烈。这一结论有助于提高对水下潜器在进行复杂三维运动过程中水动力特性的认识，并为性能优良的水下潜器运动控制系统设计提供参考。