某无人水下航行器阻力特性数值模拟

梁 晶，康宝臣，冯丽娜

（1.海军装备部，山西 太原 030027；2.山西汾西重工有限责任公司，山西 太原 030027）

0 引言

无人水下航行器是在传统水下航行器的基础上，结合人工智能的特点发展起来的一种新型的航行器。其外形设计方法同传统的鱼雷设计方法相同，都是进行大量水池和风洞试验，需要大量的人力、财力，且周期长。随着流体力学数值计算技术（CFD）、计算机软硬件、网络和数据库的迅速发展，过去十几年来，不断探索和研究水下航行器水动力数值方法，使得分析水下航行器周围流动的CFD方法功能越来越强，在计算能力和实用化方面都发生了深刻的变化。

目前，水下航行器发展处于前列的国家，在大量系列基础研究和型号试验研究的基础上，逐步形成数据库和工程设计的快速估算方法，广泛用于水下航行器的研制；开展了水下航行器性能综合评估方法和仿真技术研究，为加快水下航行器的研制提供了新的手段。而各种阻力计算和流体动力参数计算精度均有不同程度的提高，采用的软件和流体模型也各不相同[1-9]。文献[1]采用Java语言和iSIGHT软件对航行器阻力性能及分析效率进行了详细计算和优化；文献[2]采用了FLUENT流体CFD软件对小型水下航行器进行了水动力仿真分析；文献[3]采用了FLUENT软件和SST剪切应力运输模型对低速水下循迹监测航行器进行了水动力学性能数值计算；文献[4]采用了FLUENT软件和ICEM网格划分软件对AUV进行了水动力性能数值仿真；文献[5]采用了 Hexpress软件对水下航行器进行网格划分，采用Finemarine软件对无人水下航行器的阻力进行了计算；文献[6]采用了 FLUENT流体CFD软件对水下航行体超空泡减阻进行了数值模拟，研究了水下航行体空泡形态和阻力系数的变化规律，分析了航行体受力和减阻特性；文献[7]利用计算流体力学（CFD）方法计算给定外形航行器在超空泡流型下的阻力计算，分析了航行器流体动力变化并给出影响其变化的相关因素；文献[8]根据标准k-ε双方程，采用混合模型对不同空化数下的回转体自然空泡特性及阻力系数在 FLUENT软件上进行了数值仿真；文献[9]通过数学描述方法对航行器型值先行计算，然后采用（1+K）因子公式计算法计算阻力，与图谱查找的方法得出回转体水下航行器最佳外形。无论采用何种方法和何种模型，均给水下航行器的水动力仿真计算提供了思路和方法。

本文通过对带鳍舵布局无人水下航行器绕流流场的数值模拟研究，探讨使用 FLUENT流体 CFD软件对该无人水下航行器阻力特性进行预测的可行性，以便于为水下航行器外形设计提供指导。

1 外形建模与网格划分

1.1 外形建模

外形建模采用ProE三维软件，抹除水下航行器外部细小特征，细小特征不会影响流体动力参数的计算。水下航行器外形和尾部放大外形如图1所示。

图1 带鳍舵布局的无人水下航行器模型Fig.1 Model of autonomous underwater vehicle with fins and rudders

1.2 网格划分

整个流畅的计算域网格与尾部网格如图 2所示。体网格总数约为 300万，网格质量标准中的skewness≤0.75。

图2 航行器面网格Fig.2 Vehicle body surface grid

2 边界条件与数值模拟

众所周知，直接求解N-S方程是非常困难的，运动方程（N-S方程）如下：

式中：U为速度矢量；t为时间；ρ为密度；p为压强；g为重力加速度；μ流体动力粘性系数；∇为散度；∇2为拉普拉斯算子。所以通常用 2种办法对湍流进行模拟，即对N-S方程进行雷诺平均和滤波处理。

这2种方法都会增加新的未知量，因此需要相应增加控制方程的数量，以便保证未知数的数量与方程数量相同，达到封闭方程组的目的。雷诺平均N-S方程是流场平均变量的控制方程，其相关的模拟理论被称为湍流模式理论。湍流模式理论假定湍流中的流场变量由1个时均量和1个脉动量组成，以此观点处理 N-S方程可以得出雷诺平均 N-S方程（简称RNS方程）如下。

连续方程：

动量方程：

在引入Boussinesq假设，即认为湍流雷诺应力与应变成正比之后，湍流计算就归结为对雷诺应力与应变之间的比例系数（即湍流粘性系数）的计算。为了真实模拟航行器表面的微观流场，同时提高仿真计算的精度，本文采用雷诺平均N-S方程加湍流模型的方法，这种计算方法能够提供比较精确的流场数值解。

数值模拟采用ANSYS FLUENT 18.5软件求解雷诺平均N-S方程，条件设置如下[10-11]：

1）解器：三维耦合（Coupled）隐式稳态求解器，Coupled Implicit Solver算法是其它所有商用CFD软件都不具备的。该算法也对 N-S方程组进行联立求解，由于采用隐式格式，因而计算精度与收敛性要优于 Coupled Explicit方法，但是占用较多的内存。该算法另一个突出的优点是可以求解全速度范围，即求解范围从低速流动到高速流动。

2）湍流模型：Realizable k-ε模型。

3）材料：水，密度：998.2 kg/m3；粘性：0.001 003 kg/（m·s）。

4）边界条件：入口速度条件（υ=8.231 111 m/s，方向平行于航行器中心轴线），出口为自由流。

5）方程离散方法：采用压力与速度的耦合（Coupled）方法，参数的离散采用二阶精度的迎风格式。

6）收敛标准：FLUENT默认。

阻力系数计算公式如下：

3 数值模拟结果与分析

3.1 数值模拟结果

航行器表面及空间的速度分布如图3所示。

图3 航行器表面及空间的速度分布Fig.3 Velocity distribution on the surface and space of the vehicle body

水下航行器头部表面及空间的压力分布如图4所示。

图4 航行器头部表面及空间的压力分布Fig.4 Pressure distribution on the surface and space of the vehicle bow

尾端面的空间速度矢量分布如图5所示。

图5 尾端面的空间速度矢量分布Fig.5 Space velocity vector distribution of stern end plane

航行器表面静态压力分布曲线如图6所示。

图6 航行器表面静态压力分布曲线Fig.6 Curve of static pressure distribution on the vehicle body surface

图3-6给出的速度与压力分布表明：当流体流经航行器表面时，在头部附近受到阻滞，速度迅速降低，动压减少。因此，在头部附近会形成一个局部静压高于来流静压的正压驻点区（见图 4）。来流沿着头部平面对称的排开加速，压力迅速降低，在图6中，压力线陡然下落。来流绕过头部时，继续加速，压力继续降低，在平直段转折点处压力达到最小值、速度最大。此后，流速开始减小，压力开始恢复增加，在航行器平直段，压力保持常值，该值与前方来流的静压值相当。在航行器尾部，平直段与收缩段过渡区域，即尾部开始收缩时，有一个较大低压区，说明来流流过尾部曲线收缩段时存在加速。随着航行器尾部继续收缩，来流减速，压力开始恢复升高，加速区之后为轴对称扩张流动。在航行器末端，形成包含涡环在内的死水分离区如图5所示。

以上分析结果都很好地符合了回转体表面压力和速度分布的真实规律[12]。

3.2 数值模拟结果与实验结果比较

航行器阻力与阻力系数的数值模拟结果与实验结果比较如表1所示。

表1 阻力系数的数值模拟结果与实验结果比较Table1 Comparison of numerical simulation results of resistance coefficient with experimental results

从表1中可见，基于FLUENT流体CFD软件的仿真结果与风洞试验结果的误差在3%左右。初步分析，仿真计算的误差可能与网格生成的类型及质量，求解器、湍流模型、湍流参数的选择，边界条件的设置，方程离散方法的选取等因素有关。因此，通过提高网格质量和更合理的配置仿真计算参数可进一步提高仿真的精度。以上数值模拟结果可满足航行器前期设计阶段流阻力系数估算的需求。

4 结束语

采用FLUENT流体CFD软件对带鳍舵布局无人水下航行器阻力特性的预测精度为3%左右，数值模拟结果形象、直观，可利用该方法指导风洞、水洞、水池等实验，甚至可以部分代替无人水下航行器阻力特性的初步估算，对水下航行器的外形设计具有一定的参考价值。为进一步提高 FLUENT等商用CFD软件在国内水下航行器领域的应用水平，应深入探讨降低其预测流体动力参数误差的方法研究。