0 引　言

海洋是资源的宝库、交通的命脉和战略的要地，建设海洋强国也是我国的重大强国战略。随着国家对海洋开发力度的增加，无人水下机器人（UUV）也逐渐成为了一大热点。水下机器人是一种在水下环境中执行各种任务的自主（AUV）或远程操作（ROV）的机器人系统，被广泛用于海洋调查、资源开发、海底维护和军事应用等领域。由于用途不同，UUV的类别根据其尺度和排水量的大小也有超大型、大型、中型和小型之分。各类UUV的分类依据和用途如表1所示^[1]。

相比于中小型的UUV，大型和超大型的UUV具有自持能力强、多采用模块化设计以及作战能力先进的特点^[2]，还可作为潜艇的任务搭载模块，具有非常高的军事价值，也是未来各国水下战的角力点。

在水下自主航行器的研究方面，美国国防高级研究计划局（Defense Advanced Research Projects Agency，DARPA）提出了代号为“SUBOFF”的项目，SUBOFF模型总长4.356 m，最大直径0.508 m^[3]，全附体搭载十字型尾舵。该项目旨在通过使用计算流体动力学（Computational Fluid Dynamics，CFD）的技术来研究无人自主水下航行器的水动力性能。在泰勒水池中进行的SUBOFF试验，为潜艇的操纵性、流体动力性能提供了重要的实验数据，这些数据对于验证数值计算的准确性和潜艇设计的发展具有重要的参考价值。

针对传统的搭载十字型舵的UUV而言，搭载X型舵的UUV 具有更好的稳定性和机动性^[4]。尤其是对于大型或超大型的UUV（如俄罗斯的“大键琴”号、“波塞冬”号无人潜航器、美国“虎鲸”号无人自主水下潜航器等），为了有更好的稳定性和机动性，通常都会搭载X型舵。本文研究的对象正是如图1所示的搭载X型舵的大型UUV，其主要参数信息如表2所示。

计算流体力学（CFD）方法是获取UUV水动力系数最便捷和直观的方法，通过求解非定常雷诺平均模型（URANS）来获取船舶阻力，进而拟合出水动力系数，同时也可以直观地观察到流场变化。但是结合目前相关的文献[5]，针对X型尾舵以及大型UUV的相关研究较少，缺乏相关的水动力系数对运动控制仿真的支撑，因此本文将从这个角度对该模型进行研究。

由于物体所受的水动力与物体的形状及物体的运动方向有关^[6]，因此对本文所选的X舵大型UUV数值仿真获得的水动力系数用泰勒水池试验中SUBOFF号航行器的水动力系数作为参考。

1 基本理论与数值计算方法 1.1 控制方程

对于潜体水下运动，假设：

1） 流体不可压；

2） 流动时遵循质量守恒定律和动量守恒定律。

采用非定常雷诺平均模型（URANS），相应的控制方程如下：

${\left\{\begin{gathered} \nabla\cdot{{\boldsymbol{u}}}=0 \\ \frac{\partial{{\boldsymbol{u}}}}{\partial t} + \left({{\boldsymbol{u}}} - {{\boldsymbol{u}}}_b\right)\cdot\nabla{{\boldsymbol{u}}} = - \frac{1}{\rho}\nabla p - \nabla\cdot\left(\nu_{eff}\left(\nabla{{\boldsymbol{u}}} + \left(\nabla{{\boldsymbol{u}}}\right)^{\mathrm{T}}\right)\right) + {{\boldsymbol{S}}}。\end{gathered}\right.}$

(1)

式中：$\boldsymbol{u}$为流速；$p$为压力；${\nu }_{eff}$为有效黏度；$\boldsymbol{S}$为源项；$\rho$为流体密度；${\boldsymbol{u}}_{b}$为网格运动速度。有效黏度${\nu }_{eff}=\nu +{\nu }_{t}$，其中$\nu$为运动黏度为常数；${\nu }_{t}$为湍流涡黏度；由湍流模型求解获得。

1.2 湍流模型

为了使URANS方程封闭，必须通过引入新的湍流模型方程来将湍流的脉动值和时均值等联系起来^[3]。目前2种常用的湍流模型分别是涡粘模型和雷诺应力模型。其中涡粘模型常用于近壁湍流模拟，而雷诺应力模型则用于描述湍流中的应力分布。本文采用SST $k-\omega$湍流模型，这是一种基于湍动能$k$和比耗散率$\omega$的模型，用于高雷诺数流动，尤其在壁面附近具有较高的预测精度。剪切应力输运$k-\omega$考虑了基本的湍动剪切力对湍动粘度进行了修正，同时也修正了对$\omega$输运方程的交叉扩散项和混合函数，保证在近壁区域和远场都有很好的预测结果^[3]。其中动能$k$和单位耗散率$\omega$的传输方程为：

$\left\{\begin{split}\frac{\partial }{\partial t}\left(\rho k\right) + \nabla \cdot \left(\rho k\bar{\boldsymbol{u}}\right) = & \nabla \cdot \left[\left(\mu +{\sigma }_{k}{\mu }_{t}\right)\nabla k\right] + {P}_{k} - \\ & \rho {\beta }^{*}{f}_{{\beta }^{*}}\left(\omega k-{\omega }_{0}{k}_{0}\right) + {\boldsymbol{S}}_{k}, \\ \frac{\partial }{\partial t}\left(\rho \omega \right) + \nabla \cdot \left(\rho \omega \bar{\boldsymbol{u}}\right) = & \nabla \cdot \left[\left(\mu + {\sigma }_{\omega }{\mu }_{t}\right)\nabla \omega \right] + {P}_{\omega } - \\ & \rho \beta {f}_{\beta }\left({\omega }^{2}-{{\omega }_{0}}^{2}\right) + {\boldsymbol{S}}_{\omega }。\end{split}\right.$

(2)

式中：$\bar{\boldsymbol{u}}$为平均速度；${\mu }_{t}$为动力粘度；${\sigma }_{k}、{\sigma }_{\omega }$为$k$和$\omega$的模型系数；${P}_{k}、{P}_{\omega }$为结果项；${f}_{{\beta }^{*}}\mathrm{、}{f}_{\beta }$分别为自由剪切修正因子和涡流延伸修正因子；${\boldsymbol{S}}_{k}、{\boldsymbol{S}}_{\omega }$为自定义源项；${k}_{0}、{\omega }_{0}$为环境湍流值。

1.3 边界条件

边界条件的设置如下：速度入口、压力出口和对称平面。速度入口边界用于流体进入的条件，其中速度分布和流体特性是已知的。在流入边界预设条件，以计算入口的体积、动量和能量通量。压力出口边界用于设置流出时的工作压力，即流体进入环境的静压。对称平面边界则表示模拟中的对称虚拟平面。在几何和流场具有对称性的位置，可利用对称平面边界条件来缩小计算域范围。对称平面边界还可作为“滑移壁面”来处理，不存在剪切应力。

在进行拘束模型回转试验的模拟时，采用计算域全局网格运动方法，即整个计算域内的网格随物体一起运动，之间不存在相对运动，而是将运动模型与流动求解器连接在一起，并在每个时间步更新控制网格的空间位置和边界条件函数。

1.4 求解方法

采用半隐式方法（简单SIMPLE）对压力速度耦合方程组求解，通过求解压力校正方程满足对速度场的质量守恒约束。压力校正方程是基于连续性方程和动量方程构建的，通过调整压力来求解满足连续性方程的预测速度场。此方法也称为预估校正法，可通过压力校正方程获得作为变量的压力。对流项采用二阶迎风格式，时间项采取一阶离散格式。

2 大型UUV模型CFD数值仿真

本次计算工作主要基于STAR-CCM+平台展开，该软件在进行CFD计算时相对误差较小且在工程上可接受的范围内，认为可以运用在水下自主航行器后续的水动力研究中^[6]。数值仿真内容主要包括斜航试验仿真、回转试验仿真以及平面运动机构（PMM）平台试验仿真，分别用于潜体线速度水动力导数、角速度水动力导数以及加速度水动力导数的计算。计算工况设置如表3所示。

2.1 计算域设置和网格划分

在内部计算域大小对计算结果影响可以忽略的条件下^[7]，将计算域设置为矩形区，船首和四周设置为12 m，船尾设置为30 m。在计算域内采用计算较稳定的多面体网格进行划分，网格量约120万。在壁面网格处理上，划分时保证潜体大部分区域的无因次第一层网格距离y+约为30～60。

划分后的计算网格如图2所示。

2.2 结　果 2.2.1 垂直面斜航试验仿真

各工况下计算获得的流场速度矢量（部分）如图3所示。

各工况下计算获得的侧向力以及转动力矩汇总如表4所示。

2.2.2 垂直面垂荡试验仿真

各工况下计算获得的流场速度矢量如图4所示。

各工况下计算获得的侧向力以及转动力矩汇总如表5所示。

2.2.3 垂直面回转试验仿真

各工况下计算获得的流场速度矢量图如图5所示。

各工况下计算获得的侧向力以及转动力矩汇总如表6所示。

2.2.4 垂直面纵摇试验仿真

各工况下计算获得的流场速度矢量如图6所示。

各工况下计算获得的侧向力以及转动力矩汇总如表7所示。

2.2.5 垂直面固定攻角纵摇试验仿真

各工况下计算获得的流场速度矢量图如图7所示。

各工况下计算获得的侧向力以及转动力矩汇总如表8所示。

3 数据处理和分析 3.1 数据拟合

将以上获得的数据用Matlab中的cftool工具进行拟合，并进行无因次化，获得的垂直面水动力导数如表9所示。

对比该UUV水动力导数与SUBOFF试验获得的水动力导数可见，大部分水动力导数在符号以及量级上较为一致，但由于X舵与艇型的细微差异，在水动力性能上具有明显差异。

如表10和表11所示，线性水动力导数可由斜航或回转获得，同时也可由PMM运动获得，由于2种计算方式采用的动力分解模型不同，计算出的水动力导数也会存在一定差异。经由PMM试验获得的线性导数由于工况设置较少，得到的数值精度不如惯性类导数，因此在稳定性判定中将采用斜航与回转试验的系数进行分析。

3.2 垂直面动稳定性分析

在水下航行器设计的过程中，为了使其在垂直面具有直线稳定性，并不要求潜体是静稳定的，而是使衡量静稳定性的相对倾覆力臂${l}_{\alpha }^{\text{'}}$为适当正值，因此本文仅讨论UUV的动稳定性。

在扰动$\Delta \alpha 、\Delta \theta$的作用下，扰动运动的特征方程如下：

$\begin{array}{c}{A}_{3}\left(\Delta \dddot{\alpha }\right)+{A}_{2}\left(\Delta \ddot{\alpha }\right)+{A}_{1}\left(\Delta \dot{\alpha }\right)+{A}_{0}\left(\Delta \alpha \right)=0。\end{array}$

(3)

其中，

$\left\{ \begin{array}{l}{A}_{3}=\left({I}_{y}^{{'}}-{M}_{\dot{q}}^{{'}}\right)\left({m}^{{'}}-{Z}_{\dot{w}}^{{'}}\right)-{Z}_{\dot{q}}^{{'}}{M}_{\dot{w}}^{{'}}，\\ {A}_{2} = - {M}_{q}^{{'}}\left({m}^{{'}} - {Z}_{\dot{w}}^{{'}}\right) - \left({I}_{y}^{{'}} - {M}_{\dot{q}}^{{'}}\right){Z}_{w}^{{'}}-{M}_{w}^{{'}}{Z}_{\dot{q}}^{{'}}-({m}^{{'}}+{Z}_{q}^{{'}}){M}_{\dot{w}}^{{'}}，\\ {A}_{1}={M}_{q}^{{'}}{Z}_{w}^{{'}}-{M}_{\theta }^{{'}}\left({m}^{{'}}-{Z}_{\dot{w}}^{{'}}\right)-{M}_{w}^{{'}}\left({m}^{{'}}+{Z}_{q}^{{'}}\right)，\\ {A}_{0}={M}_{\theta }^{{'}}{Z}_{w}^{{'}}。\\\end{array}\right.$

将式(3)中的$\Delta \alpha$置换成$\Delta \theta$，即是关于纵倾角$\theta$的扰动运动方程。由古尔维茨判据可知，关于攻角和纵倾角$\alpha 、\theta$动稳定的充要条件是特征方程的系数满足：

${A}_{i} > 0\left(i=\mathrm{0,1},\mathrm{2,3}\right), \begin{array}{c}{A}_{1}{A}_{2}-{A}_{0}{A}_{3} > 0。\end{array}$

(4)

将系数代入式(4)，则可改写成如下不等式：

$\begin{array}{c}{K}_{vd}=\left(\dfrac{{l}_{q}^{{'}}}{{l}_{\alpha }^{{'}}}+k\dfrac{{l}_{FH}^{{'}}}{{l}_{\alpha }^{{'}}}\right) > 1。\end{array}$

(5)

式中：${K}_{vd}$为垂直面动稳定系数；${l}_{\alpha }^{{'}}$为无因次水动力中心臂（相对倾覆力臂）；${l}_{q}^{{'}}$为无因次相对阻尼力臂；${l}_{FH}^{{'}}$为无因次扶正力矩的相对力臂；k为常系数。其定义为：

${l}_{\alpha }^{{'}}=-\frac{{M}_{w}^{{'}}}{{Z}_{w}^{{'}}} ，{l}_{q}^{{'}}=-\frac{{M}_{q}^{{'}}}{{m}^{{'}}+{Z}_{q}^{{'}}}，{l}_{FH}^{{'}}=\frac{{M}_{\theta }^{{'}}}{{Z}_{w}^{{'}}} ，$

$k=\frac{-{Z}_{w}^{{'}}\left({I}_{y}^{{'}}-{M}_{\dot{q}}^{{'}}\right)\left({m}^{{'}}-{Z}_{\dot{w}}^{{'}}\right)}{\left({m}^{{'}}+{Z}_{q}^{{'}}\right)\left[{M}_{q}^{{'}}\left({m}^{{'}}-{Z}_{\dot{w}}^{{'}}\right)+\left({I}_{y}^{{'}}-{M}_{\dot{q}}^{{'}}\right){Z}_{w}^{{'}}\right]}+\frac{{m}^{{'}}-{Z}_{\dot{w}}^{{'}}}{{m}^{{'}}+{Z}_{q}^{{'}}} 。$

分别将本潜体的水动力系数和SUBOFF的水动力系数代入式(5)，得出本潜体的${K}_{vd}=1.591$，SUBOFF的${K}_{vd}=1.238$，满足潜体动稳定性衡准要求，且搭载X舵的本潜体具有更好的动稳定性。

4 结　语

本文利用STAR-CCM+对搭载X舵的大型UUV进行了仿真分析，分别通过垂直面斜航、垂荡、回转、纵摇以及固定漂角纵摇试验确定了垂直面相关的水动力系数，为X舵大型UUV在垂直面进行深度控制等各类操纵仿真奠定了基础。由于在雷诺数在高于1.4×10⁷时可以忽略尺度效应^[8]，得到的垂直面水动力系数在X舵UUV中均具有参考价值。

分析表9的数据并结合文献[9]可以得出，经由PMM试验获得的黏性类水动力系数与SUBOFF的数值差异较大，而经由斜航试验获得的数据差异较小，虽PMM试验可获得较多的试验数据但是准确度有所不足。在斜航试验中虽然需要多次改变UUV的攻角进行多次试验，但是计算结果的准确度较高。

在和SUBOFF进行对比时，斜航试验中获取的速度系数${Z}_{w}^{{'}}$与${M}_{w}^{{'}}$相差较小，均在15%以内；在回转试验中获取的角速度系数${Z}_{q}^{{'}}$和${M}_{q}^{{'}}$相差在30%以内，且数值的绝对值均小于SUBOFF，由水动力系数的意义可知当艇体的俯仰角速度改变时X舵大型UUV所受的垂向力和力矩较小；在PMM试验中获取的加速度系数${Z}_{\dot{w}}^{{'}}$、${M}_{\dot{w}}^{{'}}$、${Z}_{\dot{q}}^{{'}}$和${M}_{\dot{q}}^{{'}}$相差较大，且数值的绝对值要高于SUBOFF，即在垂向加速度和角加速度改变时垂向力和力矩的变化较大，UUV具有更高的敏感性和快速响应的性能。

在运动稳定性分析中，搭载X舵的大型UUV满足动稳定性衡准条件具有更高的垂直面动稳定系数，因此其直线稳定性更加良好。