0 引　言

海洋分层流是海洋环境中普遍存在的物理现象，海洋分层流由温度、盐度差异引起的垂向密度梯度所主导，形成具有显著非均匀特性的连续介质，这种特性在深海环境中尤为突出^[1]。分层效应会使得潜艇在真实航行中产生额外的能量耗散，并会使潜艇周围的流场结构变得更加复杂，进而影响其阻力特性与尾流稳定性^[2]。随着深海探测与军事装备技术的快速发展，潜艇等潜艇在分层海洋环境中的高效航行与隐蔽性需求日益迫切，因此探究分层流场中潜艇的尾流演化数值模拟对提升潜艇推进效率和隐蔽性具有重大研究意义。

学者以Suboff模型为研究对象并开展相关实验测量和数值模拟研究，Liu等^[3]测量了不同速度梯度下Suboff模型的阻力和艇体上压力摩擦系数等分布；刘明坤等^[4]研究了Re=1.1×10⁶时潜艇尾流的演化特性；Jimemez等^[5]测量了Re=1.1×10⁶～6.7×10⁷范围内艇体尾流速度场的分布特征；董国华等^[6]研究了密度分层流体中潜艇阻力随航速、潜深增加而减弱的规律及其机理；Liu等^[7]研究了密度分层流中潜艇阻力与航速和航深之间的关系，结果表明航速和航深越大，密度分层所造成的阻力差异就越小；Huang等^[8]对比分析了潜艇在均匀流和分层流中的尾流及自由表面特征，重点揭示了分层流对潜艇水动力特性的影响；Kumar等^[9] 采用大涡模拟方法，分析轴对称旋转体绕流特性，揭示其复杂流动结构及演化规律。

关于潜艇尾流的研究，国内外做了许多实验和仿真模拟，但涉及分层流尾流的案例却很少，本文在上述实验和仿真研究的基础上，开展分层流尾流场的数值模拟研究。

1 数值模拟方法及研究对象 1.1 数值模拟方法 1.1.1 控制方程

本文基于RANS模拟方法，对于不可压缩流体其控制方程为：

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

(1)

$ {\displaystyle\frac{\partial\left({\rho u}_{i}\right)}{\partial {t}}+\frac{\partial \left({\rho u}_{i}{u}_{j}\right)}{\partial {x}_{j}}=-\frac{\partial {p}}{\partial{x}_{i}}+\frac{\partial}{\partial{x}_{j}}\left[\mu \left(\frac{\partial{u}_{i}}{\partial {x}_{j}}+\frac{\partial{u}_{j}}{\partial {x}_{i}}\right)-\rho \overline{{u}_{i}^{{'}}{u}_{j}^{{'}}}\right]-\frac{{f}_{i}}{\rho }。} $

(2)

式中：$ \rho $为流体的密度；$ \mu $为动力粘度；$ p $为压力；$ {u}_{i} $、$ {u}_{j} $分别为速度在$ x $、$ y $方向上的时均值分量；$ {f}_{i} $为体积力；$ {u}_{i}{u}_{j} $为雷诺应力项；$ \rho \overline{{u}_{i}^{{'}}{u}_{j}^{{'}}} $为雷诺平均应力。雷诺应力张量可以由Boussinesq模型给出：

$ \overline{{u}_{i}^{{'}}{u}_{j}^{{'}}}=-{v}_{t}\left(\frac{\partial{u}_{i}}{\partial {x}_{j}}+\frac{\partial{u}_{j}}{\partial {x}_{i}}\right)+\frac{2}{3}{\delta }_{ij}k ，$

(3)

$ {v}_{t}=\frac{{C}_{u}{k}^{2}}{\varepsilon } 。$

(4)

式中：$ {v}_{t} $为涡黏系数；$ k $为湍动能；$ {C}_{u} $为经验常数，本文$ {C}_{u} $取0.09。

1.1.2 湍流模型

本文水动力定常计算选用SST(Shear-Stress-Transport)k-omega湍流模型，该模型稳定性良好，压力梯度求解精度也比较高。其方程为：

$ \frac{\partial}{\partial t}(\rho k) + \frac{\partial}{\partial{x}_{i}}(\rho k{u}_{i}) = \frac{\partial}{\partial{x}_{j}}\left({\Gamma }_{k}\frac{{\partial}_{k}}{\partial{x}_{i}}\right) + {G}_{k} - {Y}_{k} + {S}_{k} + {G}_{b}，$

(5)

$ {\displaystyle\frac{\partial}{\partial t}(\rho{\omega }) + \frac{\partial}{\partial{x}_{i}}(\rho { \omega }{u}_{i}) = \frac{\partial}{\partial{x}_{j}}\left({\Gamma }_{\omega }\frac{{\partial}_{k}}{\partial{x}_{i}}\right) + {G}_{\omega } - {Y}_{\omega } + {S}_{\omega }+{G}_{\omega b}} 。$

(6)

式中：$ {G}_{k} $为由平均速度梯度产生的湍流动能生成量；$ {G}_{\omega } $为$ \omega $的生成量；$ {\Gamma }_{k} $和$ {\Gamma }_{\omega } $分别为k和ω的有效扩散率；$ {Y}_{k} $和$ {Y}_{\omega } $为湍流引起的k和ω的耗散；$ {S}_{k} $和$ {S}_{\omega } $为自定义源项；$ {G}_{b} $和$ {G}_{\omega b} $考虑了浮力项。湍流粘度$ {\mu }_{t} $定义为：

$ {\mu }_{t}=\frac{{\rho }_{k}}{{ \omega }}\displaystyle\frac{1}{\max\left[\displaystyle\frac{1}{{\alpha }^{*}},\frac{{SF}_{2}}{{\alpha }_{1}\omega }\right]}。$

(7)

式中：$ {\alpha }_{1} $的值取0.31；$ {\alpha }^{\mathrm{*}} $为模型系数；$ {F}_{2} $和$ S $分别为混合函数和应变幅度。

1.2 研究对象

本文以Suboff航行器标准缩比模型为研究对象，模型结构如图1所示，模型主要参数见表1。

计算域如图2所示，以Suboff模型艇首位置为坐标原点建立无界长方体数值水池进行计算，图2中，艏部与入口截面之间的距离为3D而尾舵与出口界面之间的距离为17.2D，艇身与周向界面之间的距离均为6D。计算域设置参考^[9]。在数值水池中入流界面设置为速度入口，出流界面设置为压力出口，艇身左右两侧设置为对称平面，该设置能更好地避免流场边界的壁面反射效应

计算域采用棱柱层网格进行离散，为了能更加精准地捕捉到尾流场内的变化，依次对艇体指挥台以及尾舵进行面加密，并对指挥台后缘尾舵后缘以及数值水池中指挥台和尾舵后位置进行体加密。如图3所示。

影响海水密度的因素有很多，在浅海海域内海水密度变化主要由温度和盐度共同影响，而海洋中气候多变，不同气候下海水密度差别很大，为使得模拟工况更接近现实，提取近一年来的西北太平洋某处深度-密度数据，对取得数据进行处理得到海水密度均值，拟合该数据得到海水密度-深度函数式(8)和函数曲线如图4所示。

$ \rho \left({z}\right)=1\;022.144\;85+az-b{{z}}^{2}+c{{z}}^{3}-d{{z}}^{4}。$

(8)

式中：z为深度；a= 0.01774；b= 1.16328×10⁻⁵；c= 4.46117×10⁻⁸；d= 6.11289×10⁻¹³。

2 仿真计算结果验证 2.1 阻力对比验证

为验证数值模拟物理模型选取的可靠性，以及仿真计算的真实性，在本节中对不同航速下不增加海水密度函数的工况进行仿真计算，计算完成后将不同航速下Suboff潜艇阻力仿真计算数值与海军水面作战中心（NSWC）在1998年实测数据^[3]做对比（见图5）。航行速度参考实测数据，分别为3.0507、5.144、6.0962、7.1611、8.2311、9.152 m/s。

从表2可以看出，仿真计算与试验所测得的阻力均随着航速的增加而增加，在较高航速下仿真计算与实验结果误差较大，但其误差也在4%之内，由此证明仿真计算能够很好地模拟预测流场内航行器受力的状况。

2.2 尾流场时均速度对比验证

在对尾流场时均速度进行验证时，选取航速3.0507 m/s、密度998 kg/m³的工况与文献[5]对比，如图6所示，根据文献[5]对以艇尾为原点的航行体下游X/D=3、6、9、12、15等相同测量截面上的数据进行提取并加以比较得到如图7的结果。其中，r为流场内某点到航行体尾舵下游的径向距离，u为当地流体速度，U为航速。

将尾流速度场无量纲化后将仿真数据与实验数据相对比，实验与仿真数据尾迹宽度以及扩展速率几乎保持一致，但在某些位置又有所区别，这些位置集中表现在r/D=0与r/D=−0.5处。图7仿真计算曲线中x/D=12、15在r/D=0处尾流场速度时均曲线呈现出倒M形分布这是由潜艇尾舵的对称性引起的，并且仿真计算曲线在r/D=−0.5处还呈现出突变现象，出现这种现象的原因是仿真计算中考虑了指挥台的影响（见图3）。指挥台对尾流场造成了干扰引起了图中所示的突变现象。实验数据中没有上述所描述的现象是因为文献[5]的实验仅考虑了艇体与指挥台尾流相互干涉的影响，而未考虑尾舵等部件对尾流的影响如图8所示，故在图7中可明显看出仿真结果相较于实验数据在r/D=0和r/D=−0.5处尾流衰减更为明显。

另一方面，仿真计算所得到的结果要略小于实验结果，在尾流近场x/D=3和x/D=6时这种情况最为明显。文献[11]认为这种现象是湍流模型在整个流场中高估了剪切应力值并且湍粘模型普遍对逆压梯度产生的特定应力各向异性不敏感，此外，壁面附近的强流动非均匀性除却模型经过调试的平衡流动以外，任何涡粘模型都无法很好地复现这种不均匀流动现象。

通过上述潜艇阻力以及尾流速度场演化相对比，仿真计算虽与实测数据存在一定偏差，但误差均在可接受范围内，其随流向的速度亏损与尾流宽度相较于实验也比较精准，故仿真计算能够很好的模拟尾流速度场的演化特性。

3 仿真计算结果验证 3.1 分层流模型验证对比

为验证潜艇在海洋不同密度下的尾流演化，在此选取50、100 m水深为计算工况，并导入拟合海水深度-密度曲线(式8)得到在海水深度50 m下密度为1023.003 kg/m³见图9(a)，100 m水深下海水密度为1023.807 kg/m³见图9(b)。将以上2种工况在3.0507 m/s的航速下进行计算。

将分层流模型与初始计算模型（航速3.0507 m/s，密度998 kg/m³）进行阻力对比如表3所示，当航深增大时由于海水密度增大所以潜艇所受阻力也逐渐增大，50 m航深下潜艇所受阻力与实验值相差约为1.67 N，100 m航深下潜艇所受阻力与50 m工况下差值约为1.88 N，可以看出潜艇所受阻力与海水密度之间存在线性关系，但由于初始工况与实验值之间存在−0.58%的误差，将该误差考虑在50 m航深和100 m航深工况下得到预估实际值分别为105.20 N和107.10 N如表4所示，可以预测该值为潜艇在航深为50 m和100 m时实际所受阻力值。

图10(a)为流展向(x-z)截面下游流场的速度云图，从速度云图上观测3种工况差别不是很明显，因此以下游尾舵为原点，建立尾舵至出流界面之间的一条直线，用以观测不同密度下尾流速度场之间的变化，如图10(b)所示。模拟建立在相对运动的基础下，即航行器静止，流体运动，图中可以直观反映出在密度较高的流体介质中速度的亏损略小于密度较低的流体介质中。

3.2 分层流模型尾流速度场无量纲化对比

图11(a)和图11(b)展示了航深为50 m和100 m下X/D=3、6、9、12、15时数值模拟和实验尾流速度场无量纲化数值对比，2张图中航深为50 m和100 m尾流速度时均值与图7存在相同的变化趋势，但又有细微差别，如在r/D=0时100 m航深下X/D=3、6、9、12、15尾流速度场无量纲化之后的数值均略大于50 m航深工况对应数值，相对的证明在密度更高的水域中速度亏损应略小于密度较低的水域。这与图10(b)结论一致。

3.3 尾流速度亏损验证对比

为了验证不同密度下尾流场速度亏损以及艇体附加部件（如指挥台、尾舵等）对尾流的干涉，本文用艇体不同周向角和尾舵下游流展向径向距离所代表的特征截面进行分析如图6(b)所示，其中0°特征截面尾流速度亏损主要受尾舵影响，45°特征截面尾流亏损主要受艇体影响，90°特征截面尾流速度亏损主要受指挥台和尾舵影响。

文献[10]指出远尾流应处于自相拟状态（指流场的空间分布如速度、温度在特定坐标系下呈现出尺度不变性，具体内容可参见文献[12]第5.2.3节）并变得与雷诺数无关。在远尾流区，流动处于运动平衡状态，其平均流量及湍流量可仅由2个局部尺度参数完全确定：一个是局部速度尺度u₀（即最大亏损速度），另一个是半尾长度l₀（即从中心线到速度亏损值降至u₀/2处的距离），见图7，在本文中半尾宽度l₀约为0.125 m。文献[10]指出平面射流达到自相拟状态时，亏损速度应满足Gaussian分布特征：

$ f\left({\eta }\right)=\text{exp}(-\text{ln}2\times {{\eta }}^{2})。$

(9)

式中，$ {\eta }=r/{l}_{0} $。

图12～图14显示了不同密度下各特征截面上以入流速度进行无量纲化后的尾流速度亏损分布，从图中可以看出3种特征界面上尾流强度都随着下游距离的增加而迅速衰减，且3种工况下速度亏损值与式(9)所拟合的函数曲线结果较为相近，证明3种工况下远尾流均达到自相拟状态，并且在这种状态时尾流亏损基本满足Gaussian分布。

图12是0°特征界面上的尾流亏损特征，在该特征截面上，尾流主要受尾舵的影响，在图(13)中可看出在$ {r}/{l}_{0} $=−1～−2和1～2的位置速度亏损曲线波动较为明显，这处正是艇体尾舵后的速度亏损曲线，说明了尾舵给尾流造成了额外的速度亏损，并且由于艇体的对称性使得尾流亏损也相应具备了对称性，从尾流速度分布的自相拟特征分析，流体介质密度对尾流远场的自相拟状态也产生了影响，在航深为50 m和100 m工况下当X/D=12时候基本已处于自相拟状态，此时密度为998 kg/m³的工况下为尾流场还较不稳定。另一方面分层流还会影响尾流核心区的范围，图(13)中密度为998 kg/m³的工况在X/D≥12时尾流核心区处于自相拟状态的范围在$ {r}/{l}_{0} $=−1～−1区间。而航深为50 m和100 m的工况下在X/D≥12时尾流核心区的处于自相拟状态的范围$ {r}/{l}_{0} $=−1.5～1.5区间。

图13为45°特征截面上的尾流速度亏损分布，该特征截面上尾流速度亏损收到所有舰体、指挥台以及所有尾舵的相互干涉，该工况下X/D=12和X/D=15时候航深为50 m和100 m工况其速度亏损值要小于密度为998 kg/m³，另一方面从自相拟的角度来看，这3种工况都只有在尾流核心区达到了自相拟的状态，而远离尾流核心区的位置，由于流体受到艇体不同部件的干扰较多，所以尾流边缘区域只是近似的达到了自相拟状态。

图14为90°特征界面上的尾流速度亏损分布，该截面上尾流主要受指挥台和尾舵的相互干涉。主要表现在$ {r}/{l}_{0} $=−4～−2的区域速度亏损值较之函数曲线更大，该工况下高密度的流体介质较之密度相对较低的流体介质速度亏损也会更少。这3种工况也呈现出尾流核心区的位置达到自相拟的状态，在尾流边缘区由于受到的影响因素较多故未能达到自相拟状态。

4 结　语

本文基于RANS方程对suboff模型在密度为998 kg/m³、50 m航深（对应的海水密度为1023.003 kg/m³）、100 m航深（对应的海水密度为1023.807 kg/m³）条件下的定常流场进行了数值模拟研究。主要验证了密度分层流对潜艇阻力的影响，并在此基础上通过和已有实验结果对比深入研究分层流中潜艇尾流场的数值演化特性，得到以下结论：

1）分层流中随着流体介质密度的增加潜艇所受阻力会呈现出不同程度的增高，而在同航速下随着流体介质密度的增加也会使得尾流场的速度亏损值减小。

2）潜艇指挥台和尾舵对尾流之间存在复杂的干涉作用，使得航行器尾流更难进入自相拟状态，在分层流中，潜艇在密度较高的水域中尾流会相对更容易进入自相拟状态。

3）尾舵和指挥台之间的相互干涉作用对潜艇尾流的最大亏损速度造成的影响很小，但从尾流速度场无量纲化图来看二者对尾流近场也存在较为明显的影响，而在工程实例中泵喷对尾流近场的影响也很大。基于上述原因若能在设计过程中考虑艇体和指挥台、尾舵以及泵喷之间的相互干涉作用可以对提升潜艇推进效率和隐蔽性具有重大意义。