0 引　言

为探索和开发海洋，世界各国研制出各种各样的海洋装备，如AUV、ROV、深海着陆器、深海运载器等。深海运载器是一种可携带任务载荷潜入水下进行着陆作业的潜水器，具备承载大、下潜深、执行任务多样、潜航时间长等特点，其未来将趋于大排量、大深度，装载更多任务载荷及能源储备，能在更长时间跨度下执行更多试验任务^[1]。为了能够携带更多大型任务载荷入水作业，考虑将深海运载器外形做成长方体外形。在下潜至海底过程中，为节约运载器本身携带的能源，考虑采用无动力方式下潜至近海底处。无动力下潜一般可分为线性下潜和空间螺旋下潜，采用无动力螺旋下潜可大大改善AUV在水平面内的漂移距离^[2]，但在着陆过程中需要调节AUV着陆姿态。深海运载器以初始水平姿态无动力下潜相比其他下潜方式更易于深海着陆坐底，避免螺旋下潜或垂直下潜后接近海底时的姿态变换问题。

在下潜过程中，深海运载器运动姿态影响其运动速度、受力情况和坐底状态。深海运载器需要具有良好的运动稳定性，运动姿态尽量维持小倾角范围变化状态，可以保持相对稳定速度和水平坐底姿态。有许多研究者对物体的自由下落姿态及力学特性做了深入的研究，比如平板、柱体模型等。由于下落过程中运动的不确定性和强非线性，导致该类问题求解十分复杂。Andersen等^[3]研究了低雷诺数的准二维流动中自由落体平板下落的动力学行为，并将实验测量值与数值计算值进行了定量比较，解释了其非定常气动力机制。Jin等^[4]采用新的运动网格对椭圆和矩形板进行数值计算，解释了实验与数值计算的差异。Wan等^[5]与Kumota等^[6]采用实验和数值模拟方法研究了平板自由下落中的流动涡结构和发展过程，并分析了其致力机理。周琪^[7]对自由下落平板典型运动的尾迹进行了分析，平板在非平稳下落过程中形成了卡门涡街，正负涡之间的相互脱落形成交叉力矩，使平板摆动下降。

与平板自由下落类似，长方体形深海运载器在无动力下潜时也会产生横滚现象。张翔^[8]对低雷诺数水下潜航器微缩模型下潜致力机理展开研究，发现载体尾部涡流脱落造成的涡激振动是着陆器下潜产生摆动的主要原因，提高雷诺数后，涡强度减小，稳定性提高。但对于高雷诺数水下潜航器，当雷诺数大于3×10⁶时，卡门涡街会自动出现，流体会对物体产生一个周期性交变作用力。

为探究不改变运载器外形的条件下净浮力、稳心高和初始横滚角等可变参数对运载器下潜姿态的影响，本文采用数值仿真方法对运载器深海中无动力下潜过程进行分析，鉴于深海环境特点，仅考虑运载器自身参数的影响，忽略洋流等扰动因素。

1 数值计算方法 1.1 控制方程和湍流模型

深海运载器外形为长方体，具有简单几何外形，其在水中运动时类似钝体绕流，具有复杂的流场结构。钝体绕流^[9]通常伴随着流动的分离、再附以尾迹区的非定常涡脱落等复杂特征。当流体流经运载器时，流体在运载器底部形成一个分离区域，这导致了流速减小和压力增加。然后在运载器顶部，流体会再次合流，形成漩涡区域，这些漩涡区域是由于分离和再合流引起的。

对于一般的RANS方法，尤其是基于线型涡黏模型的S-A、SST等方法，虽然在干扰区以外区域对流动求解具备很好的精度，但是在马蹄涡系影响范围内，尤其是分离区域很难得到较为理想的结果^[10]。

Spalart^[11,12] 指出传统的Boussinesq线性涡黏假设未能准确捕捉雷诺应力在不同方向（流向、法向和侧向）上的各向异性，无法精确预测方管内的流动特性。为克服这一问题，Spalart提出了一种二次本构关系，重新定义了雷诺应力：

${\tau }_{ij，QCR} = {\tau }_{ij} - {C}_{r1} \left[ {O} _{ik}{\tau } _{jk} + {O} _{jk}{\tau } _{ik} \right] - {C}_{r2}{\mu }_{t} \sqrt{2{S} _{mn}{S} _{mn}}{\delta }_{ij} 。$

(1)

式中：${\tau _{ij}}$为Boussinesq线性涡黏假设的雷诺应力；${\mu _t}$为涡黏系数；${O_{ik}}$为无量纲旋转张量。定义如下：

$O_{ik}=\displaystyle\frac{\displaystyle\frac{\mathrm{\partial}U_i}{\partial x_k}-\displaystyle\frac{\partial U_k}{\partial x_i}}{\sqrt{\left(\displaystyle\frac{\partial U_m}{\partial x_n}\displaystyle\frac{\partial U_m}{\partial x_n}\right)}}。$

(2)

${S_{mn}}$为应变率张量，定义如下：

${S_{mn}} = \frac{1}{2}\left( {\frac{{\partial {U_m}}}{{\partial {x_n}}} + \frac{{\partial {U_n}}}{{\partial {x_m}}} - \frac{2}{3}\frac{{\partial {U_k}}}{{\partial {x_k}}}{\delta _{mn}}} \right)。$

(3)

式中：${C_{r1}}$和${C_{r2}}$为模型常数，分别为0.3和2.5。这种重新定义的雷诺应力${\tau _{ij,QCR}}$不仅是应变率的函数，还同时考虑了其他方向应变的影响，因此引入了雷诺应力的各向异性。

本文使用STAR-CCM+软件对深海运载器下潜进行数值模拟。根据钝体下潜特点，本文采用SST k-omega模型，本构选项为二阶QCR方法，使用Gamma-ReTheta转换模型，选择二阶迎风对流格式，自定义壁面函数，设置自由流边界为边界层厚度。合理使用二次本构关系的湍流模型能够得到更加合理的压力分布和流动分离情况^[13]。

1.2 缩比模型试验基础与数值仿真精度验证 1.2.1 缩比模型试验及问题现象

为探究运载器下潜过程中的物理现象，课题组制作了1/8缩比模型，为长方体钝体外形，尺寸L₁×B₁×H₁为1.72 m×0.3 m×0.35 m，如图1(a)所示。使用该模型进行水池试验，试验过程如图1(b)所示。运载器模型在下潜过程中存在横滚现象，而纵倾角逐渐趋于0，如图1(c)和图1(d)所示。

1.2.2 数值仿真精度验证

与试验条件保持一致，使用1/8缩比模型进行仿真，通过仿真数据与试验数据对比，验证仿真方法准确性。

设置仿真条件：运载器净浮力为3.2 kg，下潜深度为11 m，稳心高为0.03 m，初始横滚角为0°。仿真结果与试验情况对比如图2所示。可以看出，三者变化趋势基本一致，k-epsilon模型仿真结果前期较为接近试验值，但随距离增加逐渐偏离试验值；线性本构k-omega模型仿真结果中期较为接近试验值，但起始横滚角度较大，随距离增加也逐渐偏离试验值；QCR SST k-omega模型仿真结果中期横滚角偏大，总体接近试验值，相对其他模型准确性更高。

试验最终横滚幅值为41.93°；QCR SST k-omega模型仿真结果的最终横滚幅值最接近试验值，为36.11°，横滚幅值误差为13.8%；线性本构k-omega模型横滚幅值为33.3°，误差为20.5%；k-epsilon模型横滚幅值为31.82°，误差为24%。QCR SST k-omega模型误差低于15%，使用该模型进行后续运载器下潜运动规律分析。

1.3 坐标系和计算域

本文研究的深海运载器模型及坐标系，如图3所示，运载器为长方体钝体外形，尺寸L×B×H为13.7 m×2.4 m×2.8 m；排水体积V为89.997 m³；排水质量T为92328.82 kg。流体计算域及边界条件如图4所示，其中，计算域为直径10L、高240 m的圆柱体。在计算域中间设置网格加密区和运动区，下侧和周围侧面为速度进口，进口速度设定为[0，0，0] m/s，使得运载器在自身净浮力作用下无动力下潜。上侧为压力出口边界，出口压力设定为0 Pa，运载器表面边界设置为壁面条件。

1.4 网格无关性验证

考虑到计算资源及效率，在验证网格无关性时，使用定常模型进行模拟，验证模拟结果对网格的依赖性大小，从而确定是否达到了网格无关性。确定模拟结果在不同网格上变化较小且具有网格无关性后，采用动网格技术，在已验证的网格上进行模拟，可获得更精确的动态行为预测。

网格划分采用切割体网格生成器生成六面体网格，边界层处加入8层棱柱网格以有效模拟近壁面处的流动。若网格基本尺寸为BS，加密区域的网格尺寸为0.5 BS，运载器表面平面部分和曲面部分网格尺寸分别为0.1 BS、0.08 BS。为选取合适的网格基本尺寸，使用8种不同网格基本尺寸计算来流流速为2 kn，来流方向竖直向上工况下的阻力，计算结果如表1所示。

可以看出，在网格基本尺寸为1.00 m时，相对误差较小，网格数量合适，可以获得相对较高的计算精度。最终选用网格基本尺寸为1.00 m进行运载器下潜模拟，网格划分场景如图5所示。

2 深海运载器下潜仿真分析 2.1 深海运载器下潜姿态变化及力矩分析 2.1.1 深海运载器下潜姿态变化

设置运载器下潜的初始横滚角${\phi _0} = {0^{\text{o}}}$；稳心高h=0.05H；净浮力W=0.03T。其中，H为运载器高度，T为运载器排水质量。以此为研究对象进行仿真模拟，得出横滚角与纵倾角变化如图6所示。运载器无动力水平下潜时由于流体作用，会产生周期性横滚现象，而纵倾角逐渐趋于0。

深海运载器下潜是一种钝体绕流问题，此时运载器运动，流体静止。运载器下潜运动时，流体流经运载器表面向上运动，运载器两侧会周期性地脱落出旋转方向相反、有规则的双列线涡。开始时，这两列线涡分别保持自身的运动前进，接着它们互相干扰，互相吸引，干扰越来越大，形成非线性的涡街。

运载器下潜仿真涡量场如图7所示，取运载器上方A、B两处流体微团为研究对象，根据伯努利原理，可得：

${P_{\rm A}} + \frac{1}{2}\rho v_{\rm A}^2 + \rho g{h_{\rm A}} = {P_{\rm B}} + \frac{1}{2}\rho v_{\rm B}^2 + \rho g{h_{\rm B}} 。$

(4)

B处存在的涡使得B处流体速度${v_{\rm B}}$比A处速度${v_{\rm A}}$大，由伯努利原理可知A处压力${P_{\rm A}}$＞${P_{\rm B}}$，运载器上方压力不平衡。当A处发生涡旋脱落时也是同样原理。运载器上方发生周期性的涡旋脱落，引起运载器上方压力发生周期性变化，造成运载器下潜过程中产生周期性横滚运动。

2.1.2 深海运载器横滚力矩分析

运载器下潜时，在运载器顶部的流动区域，由于湍流的不稳定性，流体会产生旋转的涡流。这些旋转的涡流以交替的方式分离出来，形成了卡门涡街，使得运载器周期性摆动下潜。摆动下潜中流体旋涡与运载器相互作用引起的力矩，称为涡诱导力矩。

运载器在下潜过程中产生横滚主要受到以下4种力矩：复原力矩M_R、阻尼力矩M_f、惯性力矩M_L及涡诱导力矩M_V，如图8所示，由力矩平衡条件，可得：

${M_{{R}}} + {M_{{f}}} + {M_{{L}}} + {M_{{V}}} = 0 。$

(5)

复原力矩是运载器产生横滚角时，重力与浮力不在同一垂直线上产生，使运载器恢复到原来位置的力矩，在水下运载器产生横滚时不会发生重浮心位置变化，不会产生非线性项，其值为：

${M_{\text{R}}} = - Dh\phi 。$

(6)

式中：D为运载器排水量，t；h为运载器稳心高，m；$\phi$为横滚角，rad。

阻尼力矩是运载器与流体之间相对转动产生的，目前精确地确定阻尼力矩十分困难，采用大角度横摇时的经验公式^[14]，得：

${M_{\text{f}}} = - W{\dot \phi ^2}。$

(7)

式中：W为横滚阻尼力矩系数，$\text{tf}\cdot\text{m}\cdot{\text{s}}^{\text{2}}$；$\dot \phi$为横滚角速度，rad/s。

阻尼力矩引入的非线性项可用等效线性化方法转化为线性项，降低求解难度。

惯性力矩是主要由运载器本身的惯性力矩和附加惯性力矩组成，其大小一般与角加速度成线性关系，为：

${M_{{L}}} = - \left( {{I_{xx}} + {J_{xx}}} \right)\ddot \phi = - I{'_{xx}}\ddot \phi 。$

(8)

式中：$I{'_{xx}}$为运载器本身惯性矩和附加惯性矩之和，kg·m²；$\ddot \phi$为横滚角加速度，rad/s²。

涡诱导力矩的大小和方向取决于涡流的排列、大小、速度以及物体的形状和方向，不同的排列和流动特性会导致不同方向和大小的力矩。由于诱导力矩的非线性特征和涡诱导力矩难以表达，本文以数值方法近似求解运载器无动力下潜运动姿态及下潜速度变化。

2.2 控制参数对无动力下潜运动影响

在外形保持不变的情况下，深海运载器可变参数有净浮力W、稳心高h与初始横滚角$\phi$。研究运载器在深海无动力下潜，洋流扰动及浮力变化对运载器影响较小，因此忽略洋流及浮力变化的扰动。设置运载器下潜条件为无洋流扰动并且无浮力变化，采用控制变量法对3个参数进行探究，分析它们对运载器水平下潜运动状态的作用规律。

2.2.1 稳态最大横滚角与净浮力

保持运载器下潜的初始横滚角${\phi _0} = {0}$°与稳心高h=0.05H不变，净浮力取0.01T～0.05T内固定值，进行仿真模拟。仿真主要参数变化如表2所示。不同净浮力条件下运载器横滚运动如图9所示，忽略运载器稳定下潜部分相位角影响，便于对比不同净浮力条件下运载器横滚角幅值与横滚周期的变化。其中横滚周期为运载器左右摆动各一次后回到初始位置所用的时间。

从图9可知，随着净浮力的增大，运载器横滚角幅值逐渐增大，横滚周期逐渐变长。增大净浮力是通过增大重力实现，相当于增大了惯性力矩，同时也增大了稳态平均下潜速度。同时，流体流过运载器表面的速度也会增加，涡旋交替脱落速度增大，涡诱导力矩增大，横滚角随之增大。角速度的增加量相比横滚角的增加量小，因此随着净浮力的增大，横滚周期也相应增大。

如图10所示，稳心高不同，横滚角幅值随净浮力变化的速率不同。稳心高较大时，横滚角幅值随净浮力增加相对缓慢；稳心高较小时，横滚角幅值的变化对净浮力的敏感程度逐渐降低。运载器的稳态平均下潜速度与净浮力近似呈线性正相关关系。

2.2.2 稳态最大横滚角与稳心高

保持运载器下潜的初始横滚角${\phi _0} = {0}$°与净浮力W=0.03T不变，稳心高取0.05H～0.25H内固定值，不同稳心高引起的参数变化如表3所示。不同稳心高条件下运载器横滚运动如图11所示，忽略运载器稳定下潜部分相位角影响，便于对比不同稳心高条件下运载器横滚角幅值与横滚周期的变化。

从表3和图11可知，随着稳心高的增大，运载器下潜的横滚角幅值逐渐减小，横滚周期逐渐变短。由于运载器横滚的减弱，横滚运动所消耗的能量减小，使得下潜速度增加。增加稳心高可以增大运载器的复原力矩，使运载器能够更容易恢复到平衡状态，因此运载器在涡诱导力矩作用下发生横滚偏转时，有更大的复原力矩平衡涡诱导力矩，使得运载器偏转幅度减小。

如图12所示，净浮力不同，横滚角幅值随稳心高变化的速率不同。在净浮力较小时，稳心高增大，横滚角幅值减小，且减小速率逐渐减慢；在净浮力较大时，横滚角幅值与稳心高的关系近似于线性负相关关系，这表明在净浮力较小时，更易于通过调节稳心高来改善运载器的横滚姿态。在净浮力较大时，运载器的稳态平均下潜速度与稳心高呈线性正相关关系。在净浮力为0.01T下，由于运载器没有周期性大角度横滚运动，因此稳心高变化对稳态平均下潜速度大小的影响较小。

2.2.3 稳态最大横滚角与初始横滚角

保持运载器下潜的稳心高h=0.1H与净浮力W=0.03T不变，初始横滚角取1°～ 5°内固定值，仿真结果如图13所示。

改变运载器下潜的初始横滚角，相当于改变了钝体绕流初始迎流面的几何外形，并没有改变运载器自身参数，对初始状态有影响，对运载器水平下潜的最终运动状态没有影响。从图13可知，初始角度对运载器水平自由下潜运动状态的影响仅反映在自由下落初期2个周期内，随着时间的增加，运载器逐渐趋于稳定横滚下落状态，不同初始横滚角的稳态横滚角幅值及稳态平均下潜速度基本不变。

2.3 下潜运动参数设计原则

深海运载器携带任务载荷进行着陆坐底作业，要求下潜过程中保持一定的速度且姿态稳定，因此在设计深海运载器时，要综合考虑净浮力、稳心高、横滚角幅值及下潜速度。将仿真结果通过插值获得近似曲面，如图14和图15所示。

为满足运载器观测及着陆需求，要求深海运载器着陆时姿态角绝对值小于5°。选取图15中横滚角幅值小于5°的点，组成可选区域，如图16所示，该区域中所有稳心高与净浮力值均满足横滚角幅值小于5°。将该区域简化为二维平面，其斜边近似为直线，如图17所示，斜边方程拟合为：

$h = 7.828W + 0.01021。$

(9)

该可选区域表达式为：

$\left\{ {\begin{array}{*{20}{l}} {{\text{0}}{\text{.05}H} \leqslant h \leqslant {\text{0}}{\text{.25}H}}，\\ {0.01{\text{}T} \leqslant W \leqslant 0.032{\text{}T}} ，\\ {h \geqslant 7.828W + 0.01021} 。\end{array}} \right.$

(10)

在该区域内，下潜速度与稳心高、净浮力关系如图18所示。当稳心高为0.25H，净浮力为0.03T时，存在最大下潜速度，为2.14 m/s，并且下潜速度随净浮力的减小而减小，随稳心高的增大而减小。

3 结　语

本文基于二次本构关系，对深海运载器无动力下潜进行数值仿真，探究净浮力、稳心高和初始横滚角等可变参数对运载器下潜姿态的影响，并得出以下结论与展望：

1）二次本构关系QCR SST k-omega模型仿真结果与试验值横滚幅值误差13.8%，可用于深海运载器无动力下潜数值模拟。

2）深海运载器以初始水平姿态无动力下潜会产生横滚现象，在外形确定的情况下净浮力与稳心高会影响横滚角幅值与稳态平均下潜速度，横滚角幅值随净浮力增大而增大，且不同稳心高条件下增大速率不同；横滚角幅值随稳心高增大而减小，且不同净浮力条件下减小速率不同。初始横滚角仅在下潜初期2个周期内有影响，不影响其稳态横滚角幅值。

3）根据仿真结果拟合出横滚角幅值小于5°的净浮力与稳心高可选区域，可指导运载器与着陆器设计。当运载器净浮力与稳心高难以控制在可选区域范围内时，可考虑使用减摇机构改善运载器姿态变化。