0 引　言

海洋内波存在的条件之一是海水在多重因素的影响下产生密度分层，海洋内部跃层中动量源产生的内波垂向压力影响水面产生长时间的辐聚辐散作用，引起海洋表面粗糙度变化，这种变化难以躲避合成孔径雷达等能够在复杂情况中长时间、大范围的探测捕捉，从而暴露水下航行体的轨迹。因此，解释和探究海洋内波的存在和演化规律具有现实意义。

近年来，分层流研究逐渐受到国内外学者的重视，并逐步投入研究。对分层流内波尾迹研究最直观的方法是实验统计观察，分层流可查阅的研究也以实验为主。AH Schooley^[1]早在1963年使用水-甘油混合密度分层进行实验；Stockhausen等^[2]在1966年通过使用电导率探针测量尾流；Hopfinger等^[3]用荧光染料技术研究了3000以内雷诺数的不同傅汝德数中拖曳圆球尾流的内重力波和湍流尾迹；Voropayev等^[4]对线性分层中的射流进行实验，理论分析了偶极子涡结构的形成；Robey^[5]提出内波产生的机理将内波分为物体本身扰动影响产生的体积内波和湍流尾流产生的尾迹内波。Spedding等^{[6 − 7]}详细描述了湍流尾流在分层流场中的演化过程。随着21世纪计算机的飞速发展与应用，计算流体力学作为新的有效研究方法，与实验和理论相辅相成，推动了流体力学的发展。在传统分层流实验方法研究中，实验周期较长，人工模拟的水池密度分层条件与海洋自然形成的分层密度条件相比不具备强稳定性和回复能力，且制备分层过程对分层条件和技术要求苛刻。计算流体力学相比实验方法更加灵活，在验证实验基础上可以设置和更多复杂冗沉的工况而成为探究分层流场的重要手段。

对于密度分层流体中的CFD研究，构建密度分层主要包括多相流混合模型搭建的密度跃层及基于温度异重流模型搭建的线性分层。对于多相流密度跃层模型体积分数连续方程精度限制，收敛性弱，尤其是在同时存在自由面与晚尾迹的情况下存在模拟的不稳定性逐渐增加。而温度异重流的数值模拟将温度和密度进行耦合，引入热传导效应可忽略的假设，就能实现密度跟随温度的输运，理论上是可行性更强且稳定性更好，适用于晚尾迹的研究。国内外研究学者姚志崇^[8]等对涡对的演化过程进行研究，揭示分层流体对涡运动轨迹的影响；韩盼盼^[9]、牛明昌等^[10]、黄凤来等^[11]、马卫状等^[12]、曹留帅等^[13]采取温度异重流模型方法分别对圆柱等简单模型及潜航体进行分层流尾迹演化的研究、验证了温度异重流模型分层方法下RANS、LES等湍流模型对尾迹模拟的可行性；周根水等^[14]在商业软件中采用编译UDF输入的方法构建上层淡水、下层盐水、中间为密度线性分布的混合相；McMullan等^[15]在Openfoam中采用Boussinesq模型近似来计算流体温度变化产生的浮力，引入热膨胀系数建立不可压缩浮力方程，采用LES湍流模型进行计算，以上方法本质上是对热分层的求解。

目前文献中的工作主要集中于各几何体阻力源尾迹研究，缺少对水下动量源激发的大尺度涡旋尾迹的数值仿真，而在真实海洋的非声探潜研究中潜航器推进器动量源产生的尾迹效应不可忽略。因此，针对水下分层环境，本文采用温度异重流方法对射流及螺旋桨产生的涡对演化过程进行研究，分析该方法对动量源尾迹预报的准确性以及不同动量源产生尾迹演化的差异。

1 数值方法 1.1 控制方程

描述物理连续体基本特征时的控制方程由流体力学基本方程推导而出，包括质量守恒方程、动量守恒方程和能量守恒方程。

考虑水为不可压缩牛顿流体，流体密度$ \rho $作常数处理，质量守恒方程：

$ \nabla \cdot v=0。$

(1)

式中：v为流体的速度。

动量守恒方程：

$ \partial (\rho v)/\partial t+\nabla \cdot (\rho v\otimes v)=-\nabla \cdot (pI)+\nabla \cdot T+{f}_{{\mathrm{b}}}。$

(2)

式中：p为压力；T为粘性应力张量；f_b为作用于流体的体积力（如重力）。该模拟中浮力是驱动流动演化的关键因素，由Boussinesq近似限定密度的变化仅在式(2)中的体积力项考虑，其他项中的密度可视为常数。

能量守恒方程：

$ \partial (\rho E)/\partial t+\nabla \cdot (\rho Ev)={f}_{{\mathrm{b}}}\cdot v+\nabla \cdot (v\cdot \sigma )-\nabla \cdot q+{S}_{E}。$

(3)

式中：E为单位流体质量的总能量；q为热通量；$ \sigma$为应力张量。能量方程中忽略粘性耗散产热，仅保留温度的输运方程。

以上控制方程及假设及流场限定描述了该流体质点的受力、运动、能量转换关系。是计算流体力学得以实现的基础方程。

1.2 密度分层模型

本文采用温度异重流模型方法，设置温度与垂向高度的线性耦合以及温度与密度线性耦合。

分层强度用浮频率N描述

$ {N}={\left[-\left(g/{\rho }_{0}\right)\left({\mathrm{d}}\rho /{\mathrm{d}}z\right)\right]}^{1/2}。$

(4)

式中：g为重力加速度；ρ₀为参考密度。

以N=2.1 rad/s为例，如下：

$ T(z)=-45z+293.15，$

(5)

$ \rho (T)=-10T+3929.5，$

(6)

$ \rho (T)=\rho [T(z)]=450z+998。$

(7)

式中：z为垂向高度；$ z\in [0,0.2] $；T为温度。

以上密度、温度耦合设置均通过商业软件STAR-CCM+实现。

1.3 湍流模型

叶片在高负载状态下，近场存在的湍流脉动瞬态细节无法被RANS模型捕捉，对涡尾迹演化的扩散、强度预测不足导致涡速度分布、涡核强度等关键物理量失真。理想的LES计算湍流模型在满足网格分辨率的情况下对螺旋桨这类复杂模型的计算所需计算资源庞大。

为精确捕捉螺旋桨表面复杂流动分离及分层流环境下的尾迹演化特性，本研究采用延迟分离涡模拟（Delayed Detached-Eddy Simulation，DDES）模型。DDES在近场利用RANS模型高效处理，降低网格需求，而在尾迹涡旋形成演化部分采用LES模型精确解析大尺度涡结构，混合RANS部分采用EB k-episilon模型。

$ \partial k/\partial t+\nabla \cdot (kv)=\nabla \cdot ({\mu }_{{\mathrm{eff}},k}\nabla k)-(\varepsilon -{\varepsilon }_{0})+{S}_{k}+{P}_{{{}_{k}}}。$

(8)

式中：k为湍动能；μ_eff,k 为参考运动粘度；$ \varepsilon $为耗散率；S_k为指定源项；P_k为计算项。

$\begin{split} &\partial \tilde{\varepsilon }/\partial t+\nabla \cdot (\tilde{\varepsilon }v)=\nabla \cdot ({\mu }_{{\mathrm{eff}},k}\nabla \tilde{\varepsilon })+\\&\qquad\qquad{C}_{\varepsilon 1}{P}_{\varepsilon }/\rho {T}_{e}-C_{\varepsilon 2}^{*}\rho (\tilde{\varepsilon }/{T}_{e}-{\varepsilon }_{0}/{T}_{0})+{S}_{\tilde{\varepsilon }}。\end{split}$

(9)

式中：DDES输运方程中湍流耗散率$ \tilde{\varepsilon } $替换$ \varepsilon $，$ {C}_{\varepsilon 1} $、$ C_{\varepsilon 2}^{*} $为模型系数。

$ \tilde{\varepsilon }={k}^{3/2}/\tilde{d}。$

(10)

其中：$ \tilde{d} $的计算式为：

$ \tilde{d}={k}^{3/2}/\varepsilon -{f}_{d}\max \left[0,{k}^{3/2}/\varepsilon -\psi {C}_{{\rm{DES}}}\Delta \right]。$

(11)

式中：$ {C}_{{\rm{DES}}} $为DES模型中的长度尺度系数系数；$ \psi $为校正因子；$ \Delta $为相邻网格中心最大距离；$ {f}_{d} $、$ \psi $的计算式为：

$ \psi ={[{{C}_{\varepsilon 1}}/{{C}_{\varepsilon 2}}{{C}_{\mu }}\varphi ]}^{3/4}，$

(12)

$ {f}_{d}=1-\tanh [{(8{{r}_{d}})}^{3}]。$

(13)

相比DES模型，DDES模型加入了延迟因子，使得混合模型能够更准确判定RANS/LES转换的交界区域。

2 计算模型 2.1 几何模型

本文对射流、螺旋桨2种不同类型的动量源进行数值模拟。

射流参考Flór J B和Van Heijst^[16]的设置情况。实验在密闭水箱中进行，射流管直径为2 mm，浮频率N=2.1 rad/s，实验时水箱中水总深度在0.16～0.2 m，射流管位于水箱高度的中间位置。射流数值模拟参考实验标准设计，计算域如图1所示。

螺旋桨选择国际标准桨模INSEAN E1619七叶大侧斜螺旋桨（见图2）。该模型是典型的水下航行器推进器模型，具有确定设计参数及广泛且研究基础详实的基本参数，该模型对水下航行器隐身性能研究具有关键价值。

2.2 网格设置

在射流计算网格（见图3）中，射流孔周围网格最小尺寸为3×10⁻⁵ m，射流孔表面过3400个网格单元，逐渐过渡至演化区域，演化区域网格单元尺寸为3.75×10⁻³ m，分层流z方向上网格设置了加密区，该方向上最小尺寸为1.875×10⁻³ m，整体网格量在320万左右。

敞水E1619模型计算采用DDES湍流模型网格量约450万。在螺旋桨叶梢、叶根部均设置有加密。

螺旋桨的计算域被切割为静止域和旋转域2个部分。旋转域为桨盘面前0.4D、桨盘面后0.4D、直径为1.2D完全包裹螺旋桨表面的圆柱体，其中D为螺旋桨直径。计算采用滑移网格，旋转域和静止域内外边界网格大小一致，便于数据传递（见图4）。通过运动模块设置旋转域螺旋桨表面旋转速度实现螺旋桨的运动。

2.3 计算设置

模拟中直接设置射流孔、简化射流管道的建模，可节省大量计算资源。射流孔剖面流入液体密度大小与该高度下液体密度保持一致，入流速度参考一般层流的管道速度剖面公式。在已知流率的情况下，根据Poiseuille流动的解：

$ {u}_{\max }=\frac{{R}^{2}}{4\mu }\frac{\partial p}{\partial z}=\frac{2Q}{\text{π} {R}^{2}}=2\overline{u}。$

(14)

式中：R为管道半径；Q为流率；$ \mu $为水的粘度。

射流管道速度剖面公式为：

$ u(r)=\frac{2Q}{\text{π} {R}^{2}}\left(1-\frac{{r}^{2}}{{R}^{2}}\right)。$

(15)

式中：r为径向距圆心的距离（$ 0\leqslant r\leqslant R $）。

射流与螺旋桨均为方形计算域，射流孔设置为速度入口，右端为出口，另外四面均为对称平面。

相关文献显示在远离入流口（x/D>50）的情况下湍流脉动能量衰减至射流能量的10⁻³量级，远场的演化与入口湍流无关，且本次模拟聚焦于低雷诺数下（Re<10⁴）射流晚尾迹的演化及尺度结构特征，而非近场的湍流影响^{[17 − 18]}。故计算求解以及射流入口速度剖面的计算均采用理想层流入口，即不采用湍流模型，用准DNS方法，螺旋桨计算中采用DDES湍流模型，均采用隐式非稳态计算，时间离散格式为二阶，离散方程求解采用SIMPLE算法。

由于计算从动量源输入到晚尾迹研究时间和运动变化跨度大，迭代时间步长随计算进行调整射流动量源输入期间时间步长最低为0.0001 s，随着输入停止、流速下降、尾流演化过程逐步增加至0.04 s；螺旋桨动量源运动以180时间步每转设置起始时间步长，逐渐过渡到0.04 s。

3 结果与讨论 3.1 验证实验

本研究对照还原Flór J B等^[16]工况进行射流数值模拟，分别对：流量Q=6 mL/s，浮频率N=2.1 rad/s工况；流量Q=12 mL/s，浮频率N=2.2 rad/s工况展开数值模拟，讨论射流晚尾迹演化结构过程及尺度特征。

在分层流射流晚尾迹演化中，流体受到垂向密度梯度产生抑制作用，在水平尺度上演化出相反的一对尺度相当的偶极子涡旋，这与采用混合紫外线敏感染料曝光照片下的射流实验对比吻合，如图5所示。

另外，研究还对比了涡旋结构的尺度演化及速度演化。在流量Q=12 mL/s，浮频率N=2.2 rad/s工况中提取了涡旋结构的尺度演化与涡旋结构速度演化对比。

图6展示了涡度剖面分布在沿着穿过涡旋中心的横截面演化情况，结果显示在t=90 s时涡旋半径较实验偏小、涡度更为集中。这种情况可能是由于层流模型忽略湍流输运效应，导致涡量扩散不足而涡核集中度偏高。而在远离入流口（x/D>50）的情况下，湍流效应衰弱至可忽略，在t=140 s和t=220 s时数值模拟与实验误差显著缩小。该方法能够有效预测射流涡旋晚尾迹的演化过程。

图7为涡移动速度的演化，其中U_exp为实验中照片拍摄下的测量速度；U₀为Lamb-Chaplygin模型的预测速度；U₀^*为测量射流中心（或附近）最长例子路径相关速度再根据Lamb-Chaplygin模型预测的结果；U_fit为通过优化点与与规定曲线以最小二乘法拟合标准差最小化为三阶的结果。在t=50～210 s内，数值模拟足够预测涡旋移动速度衰减趋势。

3.2 动量源计算对比

本次模拟设计了2组动量进行对比，射流动量与螺旋桨动量计算式为：

$ {J}_{{\rm{jet}}}=\frac{\rho \cdot {Q}^{2}\cdot t}{{S}_{{\rm{jet}}}}，$

(16)

$ {J}_{{\rm{propeller}}}=\rho \cdot {n}^{2}\cdot {D}^{4}\cdot {K}_{T}\cdot t。$

(17)

式中：$ {S}_{{\rm{jet}}} $为射流孔面积；$ {K}_{T} $为螺旋桨推力系数；$ n $为螺旋桨转速；$ D $为螺旋桨直径。

模型位置均在垂向中心，因模型大小导致密度差异影响忽略不计。

数值计算设置了2组相同动量、推进时间的射流动量源与螺旋桨动量源进行对比（见表1），水平涡量演化如图8所示。

在相同动量的动量源尾涡对比中，射流与螺旋桨均在水平面内产生了尺度相当的大型涡旋结构，涡旋结构持续时间长，特征清晰。

射流与螺旋桨产生涡的尺度相当，但涡量具有明显差异。螺旋桨产生的尾流较射流水平扩散更多，晚尾迹中涡量较射流相比较小；射流产生的动量沿运动方向距离更远，涡度贴合程度更高且完全形成涡结构的速度较慢。这种现象产生的可能原因是射流的速度分布几乎在轴向上，而形成涡旋需水平两侧方向上的动量转化，螺旋桨相比射流具有更多径向动量，能更早具备涡旋产生的条件。

图9所示为2种不同动量源产生涡旋结构的涡量分布，不同动量源涡量大小具有明显差异。但射流与螺旋桨产生涡的尺度相当，且稳定性和对称性良好，涡旋结构半径大小和涡量峰值分布基本重合。

为进一步探究涡量分布，将动量源的涡量分布作归一化处理，并与Lamb-Chaplygin理论模型和实验进行比较。涡量分布归一化曲线中Lamb-Chaplygin理论模型曲线更加饱满，偶极子涡中单个涡分布自对称，涡度分布曲线呈关于涡旋中心对称的标准抛物线型，与计算和实验结果有一定差异。这种现象产生是由于实际动量源实验中产生对称涡对，涡旋单侧受另一涡旋影响，涡内外侧条件不完全对称，所以在实验中偶极子涡的单个涡分布不满足自对称，且涡对的涡量不完全相等。而通过归一化曲线对比，涡量分布数值模拟结果更加贴合实验结果。

4 结　语

本文结合温度异重流方法，采用层流及DDES湍流模型对分层流体中的不同动量源尾迹演化展开了数值模拟和讨论。该方法采用温度-密度耦合，通过设置垂向密度温度梯度来模拟海洋中的密度垂直分层。该方法中计算射流演化结果与实验相比吻合较好。

密跃层中的垂向恢复力引起动量源尾流垂向运动受限，从能量角度出发，在分层流中尾流垂向上的势能转换为水平方向上的动能，从而导致动量在水平分布上更广泛。动量源尾流长时间演化为大型相反对称的涡结构。同动量情况下，射流尾流产生的涡旋强度较大，螺旋桨则动量分布更为广泛，且相比射流能够更早形成大型涡结构。射流与螺旋桨动量源产生的涡结构分布归一化尺度与实验结果相比吻合较好。研究说明动量源类型及尺度与涡对形状、形成时间、涡量大小等结果密切相关，而不同动量源产生涡旋结构的涡量归一化分布相同，可以从更多角度如控制动量源特征尺度、改变分层强度条件等方面进行更多讨论。

所得结果对数值模拟方法研究分层流体尾迹、其余动量源尾迹研究对比、水下航行体尾迹探测等提供参考。