0 引　言

航行体水下运动时，随着运动速度的不断增大，在航行体壳体上具有最小压力 $ {p_{\min }} $ 点处，将具有最大的减压系数，且有

$ {p_{\min }} = {p_0} - {\xi _{\max }}\frac{{\rho v_0^2}}{2}\text{。} $

(1)

式中： $ {p_0} $ 为航行体上任一点的静水压力； $ {\xi _{\max }} $ 为最大减压系数； $ \rho $ 为海水的密度； $ {v_0} $ 为航行体的运行速度。

由此可见，随着航行体运行速度 $ {v_0} $ 的不断增大，壳体上最小压力点处的压力 $ {p_{\min }} $ 将不断下降。当 $ {p_{\min }} $ 下降到某一数值时，便会出现空泡^[1]。空泡的产生直接影响航行体的水动力特性及稳定性控制，当空泡进入高压区时，空泡迅速溃灭并产生极高的溃灭压力，对邻近溃灭域的固体壁产生巨大的破坏作用。因此，获取航行体在水中运动时空泡的分布情况对水动力、水下载荷研究具有重要意义，但在实际工程研制中，由于受试验设施、经费、时间等方面条件限制，很难通过试验方法获取。本文利用数值计算方法对某航行体水下不同运动速度的空泡状态进行研究，对该航行体水动力特性及水下载荷研究提供参考。

1 计算模型 1.1 控制方程

本文研究的水下航行体结构具有对称性，计算模型采用二维轴对称模型，以速度入口作为边界条件，假设海水为不可压流体，整个流场以连续性方程和N-S方程为控制方程^[2]：

$ \frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \left( {\rho u} \right)}}{{\partial x}} + \frac{{\partial \left( {\rho v} \right)}}{{\partial y}} = 0 \text{，}$

(2)

$ \frac{{\partial \left( {\rho u} \right)}}{{\partial t}} + u\frac{{\partial \left( {\rho u} \right)}}{{\partial x}} + v\frac{{\partial \left( {\rho u} \right)}}{{\partial y}} = \mu \left[ {\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {y^2}}}} \right] - \frac{{\partial p}}{{\partial x}} - \rho g\text{，} $

(3)

$ \frac{{\partial \left( {\rho v} \right)}}{{\partial t}} + u\frac{{\partial \left( {\rho v} \right)}}{{\partial x}} + v\frac{{\partial \left( {\rho v} \right)}}{{\partial y}} = \mu \left[ {\frac{{{\partial ^2}v}}{{\partial {x^2}}} + \frac{{{\partial ^2}v}}{{\partial {y^2}}}} \right] - \frac{{\partial p}}{{\partial y}}\text{。} $

(4)

式中：u和v分别为流体质点在x和y方向的速度分量；p为压力；g为重力加速度。

1.2 空化模型

采用混合物模型对计算模型中多相流进行计算，混合物模型将气、液混合物考虑成一个整体，并对一些本构进行假设。混合物质量守恒方程、动量守恒方程和能量守恒方程由考虑浓度(体积分数) 的扩散方程封闭^[3]，其质量方程为：

$ \frac{{\partial {\rho _m}}}{{\partial t}} + \nabla \cdot ({\rho _m}{C_m}) = 0\text{。} $

(5)

式中: 下标m代表混合物， $ {\rho _m} = \displaystyle\sum\limits_{n = 1}^2 {{\rho _n}{\alpha _n}} $ ； $ {C_m} = \displaystyle\sum\limits_{n = 1}^2 {\dfrac{{{\rho _n}{C_n}}}{{{\rho _m}}}} $ ； $ {P_m} = \displaystyle\sum\limits_{n = 1}^2 {{\rho _n}{\alpha _n}} $ 。动量方程为：

$ \begin{split} & \frac{\partial }{{\partial t}}({\rho _m}{C_m}) + {\rho _m}({C_m} \cdot \nabla ){C_m} = \hfill \\ & - \nabla ({p_m}) + \nabla \cdot (\overline{\overline \tau } + \overline{\overline {{\tau _t}}} ) + {M_m} + f \text{，} \end{split} $

(6)

式中: $ {M_m} $ 为界面动量源项。扩散方程为：

$ \frac{{\partial {\alpha _n}{\rho _n}}}{{\partial t}} + \nabla \cdot ({\alpha _n}{\rho _n}{C_m}) = {\varGamma _n} - \nabla \cdot ({\alpha _n}{\rho _n}{C_{12}}) \text{。} $

1.3 湍流模型

采用k-ε湍流模型。标准k-ε模型是基于湍流动能k和湍流能量耗散率ε的输运方程的半经验公式^[4]：

$ \begin{split} \frac{{\partial (\rho k)}}{{\partial t}} + \frac{{\partial (\rho k{u_i})}}{{\partial {x_i}}} =& \frac{\partial }{{\partial {x_j}}}\left[\left(\mu + \frac{{{\mu _t}}}{{{\sigma _k}}}\right)\frac{{\partial k}}{{\partial {x_j}}}\right] + {G_k} +\\ & {G_b} - \rho \varepsilon - {Y_M} + {S_k} \text{，} \end{split} $

(7)

$ \begin{split} \frac{{\partial (\rho \varepsilon )}}{{\partial t}} &+ \frac{{\partial (\rho \varepsilon {u_i})}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_j}}}\left[\left(\mu + \frac{{{\mu _t}}}{{{\sigma _\varepsilon }}}\right)\frac{{\partial \varepsilon }}{{\partial {x_j}}}\right] + \\ & {C_{1\varepsilon }}\frac{\varepsilon }{k}({G_k} + {C_{3\varepsilon }}{G_b}) - {C_{2\varepsilon }}\rho \frac{{{\varepsilon ^2}}}{k} + {S_\varepsilon }\text{。} \end{split} $

(8)

其中： $ {G_k} $ 为由于平均速度梯度引起的湍动能 $ k $ 的产生项； $ {G_b} $ 为由于浮力引起的湍动能 $ k $ 的产生相； $ {Y_M} $ 代表可压湍流中脉动扩张的贡献； $ {S_k} $ 和 $ {S_\varepsilon } $ 为用户定义的源相。

另外，在该模型中，湍动粘性系数 $ {\mu _t} $ 可以表示为：

$ {\mu _t} = \rho {c_\mu }\frac{{{k^2}}}{\varepsilon }\text{，} $

(9)

$ k $ 定义为：

$ k = \frac{1}{2}({u'_x}{u'_x} + {u'_y}{u'_y})\text{，} $

(10)

式中： $ {u'_i} $ 为i（i = x，y）方向的脉动速度。

$ \varepsilon $ 定义为：

$ \varepsilon = {c_D}\frac{{{k^{\frac{3}{2}}}}}{l}\text{，} $

(11)

式中： $ {c_D} = 1 $ ， $ l $ 为湍流尺度。

标准 $ k - \varepsilon $ 模型中的其他系数为：

$ {c_\mu } $ =0.09， $ {\sigma _k} $ =1.00， $ {\sigma _\varepsilon } $ =1.30， $ {C_{1\varepsilon }} $ =1.44， $ {C_{2\varepsilon }} $ =1.92。

2 数值计算与结果分析 2.1 模型建立及边界条件

计算模型以某航行体实际尺寸为依据，并根据计算需要顺时针旋转90°，对称轴为EF, 模型网格图如图1所示。GHIJKLM区域为航行体，曲面GHI为航行体头部，JKLM为航行体尾部，整个航行体处于海水中，AB为速度入口，CD为压力出口，重力加速度方向为 $ \overrightarrow {FE} $ 方向。数值计算采用Mixture多相流模型和标准k-ε湍流模型，压力-速度采用Simple方法进行迭代求解,为提高计算精度，所有方程的对流项都采用2阶迎风格式进行离散。

2.2 计算结果及分析

数值计算中，航行体头部及尾部采用Quad-Map类型网格，其他区域采用结构化四边形网格对整个流场进行计算区域的离散化，并对航行体附近区域进行局部加密。本文计算了航行体5种速度下的空泡分布情况，如图2～图6所示。

由图2～图6可知，随着航行体速度不断增加，空泡越来越长，对航行体的水动力特性影响也变大，空泡溃灭时对航行体的冲击载荷也将变大。当速度为26 m/s时，其肩部产生的空泡已经发展至航行体中部；当速度为29 m/s时，其肩部产生的空泡几乎将整个航行体包围，航行体尾部也出现了较大空泡；当速度为33m/s时，其肩部空泡比29 m/s时更长，已经与尾部空泡连成一体。航行体各状态下的阻力系数见表1。

由表1可知，速度为20～26 m/s时，其阻力系数变化不大，当速度达到29 m/s时，受尾部空泡影响，航行体头部与尾部压差增大，其阻力系数明显增大。

3 结　语

本文采用Mixture方法数值模拟了某航行体水下不同运动速度下的空泡状态，清晰模拟出各种状态下空泡的分布情况，并计算了各工况下的阻力系数。通过计算表明，当航行体尾部出现空泡时，其阻力系数有明显变化。但是受条件限制，没有相应的空泡试验数据对计算模型进行校核，其真实分布情况有待于进一步研究。但由于空泡数值计算近年来发展较为成熟，其误差范围可控制在15%之内，因此本文可对该航行体水动力特性及水下载荷研究提供一定参考。