潜体在密度分层流中航行时，其水动力特性较为复杂。众所周知，潜体在单一均匀流体中航行时，会在自由液面处兴起波浪，形成凯尔文波系；实际海水中由于温、盐差的存在会出现分层现象，潜体航行于多层流之间时，由于体积效应及湍流尾迹效应会在密度分层处产生内波^[1]，其周期较长，振幅较大，会对潜体的航行性能产生较大影响^[2]。故研究潜体在密度分层流这一特殊海况下航行时的受力特性具有十分重要的意义。

文献[3-6]对于运动物体在分层流体中激发内波的现象进行了研究。Gou等^[7]采用时域高阶边界元法研究了两层流体中的波衍射问题。Song等^[8]对内孤立波与海洋结构物的相互作用进行了数值模拟。Lin等^[9]、Bonneton等^[10]对密度分层流进行了试验研究。针对物体在密度分层流体中的受力情况，Motygin等^[11]基于势流理论推导了二维物体分别在上、下层流体中运动所产生的兴波阻力公式。Miloh等^[12]采用格林函数法得到了半潜扁球体在2层流体中阻力及自由表面扰动的解。Grue^[13]对“弗雷姆”号的航行过程进行研究，发现的受力变化规律与Ekman^[14]一致。勾莹等^[15]采用模型实验的方法研究了箱型结构在2层流中拖动时的阻力特性，并与单层流中拖航阻力实验结果进行了对比研究。

综上可知，针对潜体在分层流体中的研究主要以运动激发内波的研究居多，而对由于内波存在引起的潜体受力变化研究并不多见。与受力相关的研究也大多关注水面舰船，本文着重研究潜体航行于密度分层流中的受力特性。本研究基于粘流理论、压力-速度耦合项求解采用Simple算法，通过定义各流体项体积分数以及计算域的压力分布，建立了密度分层流数值水池模型，针对近水面潜体航行特性开展了模拟研究。

1 密度分层流数值模型基本理论 1.1 控制方程

本模型基于雷诺平均的N-S方程，不可压缩流体的控制方程为：

$ \frac{\partial U_{i}}{\partial x_{i}}=0 $

(1)

$ \begin{aligned} \frac{\partial U_{i}}{\partial t}+U_{j} \frac{\partial U_{i}}{\partial x_{j}}=&-\frac{1}{\rho} \frac{\partial p}{\partial x_{i}}+\frac{\partial}{\partial x_{j}}\left[\nu\left(\frac{\partial U_{i}}{\partial x_{j}}+\frac{\partial U_{j}}{\partial x_{i}}\right)\right]-\\ & \frac{\partial \overline{u_{i}^{\prime} u_{j}^{\prime}}}{\partial x_{j}}+f_{i} \end{aligned} $

(2)

式中：U_i是时均速度；u′_i为脉动速度；$ \overline {{{u'}_i}{{u'}_j}} $为雷诺应力项；f_i为质量力强度；p为压强；ρ为流体密度；ν为流体的运动黏度。

式(1)为连续性方程，式(2)为动量方程，两者均为时均计算处理后的形式。

1.2 湍流模型

需要引入湍流模型使方程(2)封闭。本研究采用Realizable k-ε湍流模型，该模型稳定性良好，压力梯度求解精度高，工程应用比较广泛，其湍动能及湍动能耗散率输运方程为：

$ \rho \frac{\mathrm{d} k}{\mathrm{~d} t}=\frac{\partial}{\partial x_{t}}\left[\left(\mu+\frac{\mu_{\mathrm{t}}}{\sigma_{k}}\right) \frac{\partial k}{\partial x_{i}}\right]+G_{k}+G_{b}-\rho \varepsilon-Y_{M} $

(3)

$ \begin{array}{r} \rho \frac{\mathrm{d} \varepsilon}{\mathrm{d} t}=\frac{\partial}{\partial x_{i}}\left[\left(\mu+\frac{\mu_{t}}{\sigma_{k}}\right) \frac{\partial \varepsilon}{\partial x_{i}}\right]+\rho C_{1} S \varepsilon- \\ \rho C_{2} \frac{\varepsilon^{2}}{k+\sqrt{v \varepsilon}}+C_{1 z} \frac{\varepsilon}{k} C_{2 z} G_{b}, \end{array} $

(4)

式中：μ为分子扩散所造成的动力粘性；μ_t为湍流粘性系数；C₁=max[0.43, η/(η+5)]；η=Sk/ε；S为平均应变率；C₂和C_1z为常数；C_2z为浮力对耗散率影响的函数；G_k为平均速度梯度产生的湍动能；G_b为由浮力产生的湍动能；Y_M表示湍流脉动膨胀对总耗散率的影响；σ_k和σ_z分别是湍动能k和耗散率ε的湍流普朗特数，其他参数取值为：C_1z=1.44，C₂=1.9，σ_k=1.0，σ_z=1.2。

1.3 平面进行波理论

航行体兴波存在首波系和尾波系，且两波系中的横波在艉部之后相遇而叠加，称为兴波干扰。艏横波的第1个波峰和艉横波第1个波峰之间的距离称为兴波长度，用mL表示，m为与船型和Fr有关的系数。由平面进行波理论可得，兴波阻力系数C_w为：

$C_{w} \propto \frac{R_{w}}{\frac{1}{2} \rho v^{2} S}=\left[C+D \cos \left(\frac{2 \pi m L}{\lambda}\right)\right] \cdot\left[\frac{v}{\sqrt{g L}}\right]^{4} $

(5)

式中：ρ为水密度；g为重力加速度；λ为波长；v为航速；S为湿表面面积；L为潜体长；C、D为常数。

由式(5)可知，对于给定船型，其兴波阻力系数C_w与长度Fr的4次方成比例。

2 密度分层流数值模型设置与验证 2.1 数值模拟相关设置

鉴于所研究的问题及结构均对称，故在潜体中纵剖面处采用对称面条件来提高计算效率。计算域总长为15 L，宽为3.5 L。基于欧拉多相流模型，通过UDF的方式分别给定空气、淡水、盐水三相的密度及初始体积分数。

计算域前端边界条件为速度入口，后端为压力出口，使用切割体网格，对自由表面、密度分界面以及潜体表面曲率变化较大等处进行网格加密以准确捕捉流场特性，在潜体重心处进行网格定位。调节网格增长率、棱柱层厚度等参数控制网格质量，经多次数值实验，最终总网格数为400万左右，其中棱柱层厚度为0.018 m，棱柱层数为6，增长率为1.2。计算域网格如图 1所示。

模型采用全附体标准SUBOFF标模，其总长L为3 m，最大直径0.174 m。计算域中流体密度沿垂向变化，分层情况如图 2所示。研究中SUBOFF潜体以不同速度U航行于不同位置处。潜体潜深为d₁，指潜体中心到自由液面的垂向距离；d₂为潜体位于淡水层中时其中心距内波面的距离；h₁、ρ₁分别为淡水层流体的深度及密度；h₂、ρ₂分别为盐水层流体的深度及密度。

2.2 网格及时间步长收敛性验证

采用航速Fr=0.5的全附体潜体进行收敛性研究，此时各参数为淡水层h₁=0.24L，密度取为997.561 kg/m³；盐水层高h₂=0.67L，密度为1 020 kg/m³，d₁=d₂=0.12 L。通过改变网格的基础尺寸来控制网格数量，分别取网格基础尺寸Mesh A为0.07 m，Mesh B为0.065 m，Mesh C为0.056 m，Mesh D为0.052 mm进行网格收敛性研究，取时间步长0.01 s。

表 1是网格数量改变对阻力计算结果的影响。其中误差是指采用某一网格量计算结果与采用最精细网格Mesh D所得计算结果之间的差别，以百分比形式表示；表中R_t代表总阻力，R_f代表摩擦阻力。由表 1可以看出，当网格基础尺寸由0.056 m降低到0.052 m级别时，对总阻力及摩擦阻力计算结果带来的影响并不大，说明此时网格收敛。考虑计算效率，在之后的计算中均采用Mesh C网格。

采用4种时间步长方案对上述计算模型进行了模拟，此时网格划分采用Mesh C方式进行，结果见表 2，表中误差为采用某一时间步长计算结果与采用最小步长Time D所得计算结果之间的差别。可见当时间步长从0.010 s降低到0.005 s，总阻力的误差为0.352%，摩擦阻力误差为0.022%，说明此时时间步长的改变对计算结果的影响不大，故后续的研究中将会采用Time C方案。

2.3 计算模型有效性验证

关于分层流中潜体受力特性的CFD研究较少，本文基于STARCCM+进行仿真，对数值模型进行验证。根据文献[15]中的实验装置参数进行建模。其中水深为0.6 m，上层淡水密度ρ₁=997 kg/m³；下层是密度ρ₂=1 024 kg/m³的盐水；上层淡水深h₁为0.3 m，下层盐水深h₂为0.3 m。实验模型为60 cm×45 cm×35 cm的方型浮箱，吃水为0.2 m，速度分别为0.14、0.16、0.18、0.2、0.24 m/s。

数值研究结果与试验结果^[15]的对比见图 3，其中数值模拟计算结果为箱体在双层流中拖航时所受的总阻力。由图 3可以看出，数值模拟结果与文献中的试验数据吻合较好，说明本文所采用的CFD计算模型对潜体受力的模拟具有较好的精度。

3 模拟结果展示与参数影响分析 3.1 航速变化 3.1.1 航速变化对潜体所受阻力的影响

潜体位于淡水层中，潜体与交界面的距离会对其水动力性能带来较大影响^[16]。为系统地研究潜体所处位置对其阻力曲线的影响，本文模拟采用3个工况来展示与各交界面距离不同时潜体的受力特征，见表 3，各算例的Fr均取为0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9、1.0。

设计工况A、B对比来观察淡水层厚度对阻力曲线的变化影响；工况B、C的对比可以观察潜体在与内波面距离不变情况下，与自由液面距离(即潜深)变化对阻力曲线变化的影响如图 4所示。

由图 4(a)可知，在3种工况下，随着航速增加，潜体的总阻力均呈现不断增加的趋势。对于工况A，即淡水层的厚度较小、潜体距离2个交界面均较近时，其总阻力较其他工况最大，且其增幅会在高航速区(0.7≤Fr≤0.9)变得平缓，当速度继续加大，阻力会继续加速上升。结合图 4(b)、(c)分析原因为在密度分层流中，当潜体在与交界面距离较近处航行时会兴起较大波浪，进而引发较大的兴波阻力，兴波阻力为剩余阻力的主要成分，当速度较低时，兴波阻力与摩擦力均随速度增加而加大，故总阻力加速上升，而当航速超过Fr=0.5时，兴波阻力系数不断下降，这会在一定程度上导致总阻力的增速降缓。而随着航速的继续加大，摩擦阻力的持续增加又导致了总阻力的持续上升。

由图 4(b)可知，随着航速增加，潜体所受摩擦阻力会不断上升；由于摩擦阻力主要受航速及流体密度的影响，故3个工况下摩擦阻力随速度的变化几乎一致，阻力曲线重合度很高。由图 4(c)可见，3个工况下的兴波阻力系数表现出相同的趋势，即随着潜体航速的增加，阻力曲线先上升后下降，在Fr=0.5附近达到峰值，这是由于在式(6)中有一余弦项cos(2πmL/λ)的值是在-1.0~+1.0变动，这是由潜体艏、艉的兴波干扰情况所决定的，故兴波阻力系数C_w与Fr曲线存在波动，当Fr=0.5时潜体艏、艉兴波发生了明显的不利干扰。

为直观的展现图 4(c)中现象，图 5为工况B中的潜体在以Fr=0.3，0.5，0.7航行时在自由液面处的兴波。可见，当Fr=0.5时，自由液面处的兴波具有较高的峰、谷值，说明此时在潜体艏、艉兴波发生了不利干扰，这与图 4(c)中的分析一致。

3.1.2 航速变化对潜体表面压力的影响

针对工况A、B、C，分别选择Fr=0.3，0.5，0.7速度展示潜体表面压力分布情况。图 6中潜体沿X轴负方向前进；图 6为3种工况下潜体整体压力随航速的变化情况，3个工况中潜体表面的一些拐点处压力值变化相近且集中，不易区分，故通过图 7展示潜体以工况B航行时，其艏、艉及附体周围存在的一些拐点处(图 6(b)中方框)的压力随航速的变化情况以直观分析。

综合图 6、图 7可见，不同工况下，潜体表面压力分布展现出相同的形状。潜体表面压力主要可分为静压、动压2部分，动压主要由速度变化引起，而静压则是指潜体表面单位面积上方流体所具有的静压能，主要是由潜深决定。同一航速下，工况A的潜体的表面压力会低于工况B，而工况C时潜体表面压力最大；分析其原因主要在于，工况A相比B、C，其潜体距离自由液面比较近，是潜深的加大增加了潜体表面的静压能从而导致了压力的增加。

此外，随着航速的增加，潜体表面压力逐渐加大，其中潜体艏、艉部以及驾驶舱、尾翼等表面曲率变化较大处压力随航速的变化较为剧烈，而潜体主体部分的压力基本维持稳定；分析原因为速度的改变主要影响潜体表面一些拐点处的流体动压，进而使得这些部位的压力变化较大。

3.2 潜体与上下交界面距离的变化

较单一流体，航行于密度分层流中潜体，涉及到自由液面(上交界面)与内波面(下交界面)；本节研究了潜体与上交界面距离d₁(即潜深)及与下交界面的距离d₂变化对其受力特征的影响。

3.2.1 潜体与上下交界面距离对阻力的影响

工况D、E中相同的设置参数有ρ₁=997.561 kg/m³，ρ₂=1 020 kg/m³，h₂=0.67 L，Fr=0.5。不同的参数设置见表 4，工况D、E分别用于研究潜体与自由液面距离d₁、与内波面距离d₂的变化对其水动力性能的影响。

图 8为潜体与上、下界面的距离(d₁、d₂)改变时其阻力变化情况，其中各图实线为d₁改变时的阻力曲线，虚线为d₂改变时的阻力曲线。

由图 8中d₁变化可见，随着潜深的增加，潜体所受总阻力逐渐下降；而在位置变化时并没有改变流体的密度以及潜体的航速，故潜体所受摩擦阻力变化极小；此外，随着潜深的增加，潜体的兴波阻力系数逐渐下降，结合图 8(a)中总阻力的下降，说明当潜深改变时，兴波阻力为总阻力的主要成分。

由于兴波阻力是潜体对其兴起的波浪提供能量而产生的，故随着潜深的增加，潜体的兴波也将逐渐减弱。这也说明对于额定功率下航行的潜体，下潜越深，其航速越高。可以预计，当潜深足够大时，其阻力将变得特别小。

由图 8中d₂变化可见，随着d₂的加大，总阻力出现下降的趋势，但相比d₁改变带来的阻力变化，d₂改变导致的阻力下降幅度很低。与d₁改变时的情况相同，潜体所受的摩擦阻力并未随潜体与内波面距离d₂的改变发生较大的变化。由图 8(c)可见，潜体的兴波阻力系数随着d₂的增加也出现降低的趋势，兴波阻力系数曲线与总阻力曲线一致的变化说明d₂变化时兴波阻力同样为总阻力的主要成分。与图 8(a)中一致，d₂改变带来的阻力系数变化幅度低于d₁。

同时发现，在距离d₂达到0.15L之前，潜体所受的总阻力及兴波阻力系数下降速率远大于d₂超过0.15L之后，当d₂超过0.15L后，再继续增加d₂，总阻力及兴波阻力系数几乎没有变化。这说明，当潜体与内波面之间的距离较小时，距离的变化对其阻力的影响比较明显，而当潜体航行位置离内波交界面较远时，海洋环境中的密度分层将几乎不会对潜体阻力产生影响。

综上可知，潜体与2个交界面的距离改变对其阻力特性影响较大，尤其当潜体靠近交界面时，其兴波阻力会变大，进而导致总阻力的增加。比较来看，潜体与自由液面的距离d₁(即潜深)变化对阻力所带来的影响幅度要大于潜体与内波面距离d₂变化所带来的影响。

3.2.2 潜体与上下交界面距离潜体表面压力的影响

图 9为针对工况D、E研究潜体表面压力分布情况。由图 9可见，不同位置工况下潜体表面压力分布形状基本相同，观察每一个工况下潜体的表面压力分布情况，潜体艏、艉及潜体附体处压力数值会比较高，而其余部分比较低。

由图 9(a)可知，当d₂固定时，随着潜深d₁增加，潜体表面压力整体逐渐增加，这主要是由于潜深增加，作用在潜体表面的静水压强增加导致的。与图 6、7中航速增加所带来的影响不同，潜深的增加并没有使潜体艏、艉等拐点处的压力数值产生较大的改变，而是使潜体压力整体上增加。说明航速主要影响动压，进而对潜体表面拐点处的压力影响较大；而潜深主要影响静压，进而影响潜体整体压力数值；相比来说，动压带来的压力数值变化比较明显。

由图 9(b)可见，当d₁不变时，d₂的改变对潜体表面压力的影响很小。这是由于d₁不变使得潜体所受静压不变，且航速固定，动压也不会产生较大的变化。

4 结论

1) 潜体的受力特性随航速的影响较大。随着航速的增加，潜体总阻力及摩擦阻力均不断上升，由于潜体艏、艉兴波产生了不利干扰使得在Fr=0.5附近时潜体的兴波阻力系数最大，此时潜体会兴起较大的波浪。

2) 潜体与上下交界面之间距离的改变同样会对潜体水动力性能产生较大的影响。当潜体靠近交界面时，其遭受的阻力逐渐增加，兴波阻力为总阻力的主要成分。潜体与自由液面的距离即潜深的变化对潜体受力的影响要大于其与内波面距离变化时带来的影响。

3) 潜体航速及位置的不同会对潜体表面压力带来影响。其中速度变化主要通过改变动压进而对潜体艏、艉及附体等表面曲率变化较大处的压强产生影响；而位置的变化主要通过改变静压，进而影响潜体整体压强。

4) 本数值模型具有较好的精度，可为密度分层流中潜体受力研究提供一种有效的分析手段。