0 引　言

采用超空泡减阻技术的水下航行体，其减阻率可达95%以上，因此该项技术获得广泛应用研究^[1]。为了描述自然空化过程中汽-液之间的相变过程，目前发展出了多种空化模型，这其中包括基于Rayleigh-Plesset气泡动力学方程的Singhal模型^[2]、Zwart模型^[3]，基于经验的Kunz模型^[4]、Hosangadi模型^[5]等。不同的空化模型基本都包含可人工调节的相变系数，对于不同的空化模型、不同类型的空化问题，相变系数对数值计算结果影响很大。

研究人员针对半球头圆柱研究了Singhal空化模型中不同相变系数对计算结果的影响，指出不同空化数下相变系数的设置规范^{[6 – 7]}。Liu等^[8]在研究水泵叶轮的空化问题后认为，Zwart空化模型中的蒸发系数和冷凝系数的取值对叶轮上的空泡形态和尺度有重要影响。事实上，不同头型（平头、球头、锥头等）会导致周围流场产生不同的脉动特性和漩涡特性，致使初生空泡形态和空化产生的难易程度不尽相同^[9]。因此在利用空化模型模拟不同头型的空化现象时，相变系数也不是一成不变的。

本文采用基于CFD的数值计算方法，研究利用Zwart空化模型模拟不同头型的空化现象时相变系数对数值计算结果的影响，为进一步提高空化问题的数值计算精度提供指导性依据。

1 数值计算理论 1.1 N-S方程

针对汽-液两相流，混合介质的N-S方程可以表示为：

$\begin{array}{l}\displaystyle\frac{{\partial {\rho _m}}}{{\partial t}} + \frac{{\partial \left( {{\rho _m}{u_j}} \right)}}{{\partial {x_j}}} = 0\text{，}\\\displaystyle\frac{{\partial \left( {{\rho _m}{u_i}} \right)}}{{\partial t}} + \frac{{\partial \left( {{\rho _m}{u_i}{u_j}} \right)}}{{\partial {x_j}}} = \\ \displaystyle- \frac{{\partial p}}{{\partial {x_j}}} + \frac{\partial }{{\partial {x_j}}}\left( {\left( {\mu + {\mu _t}} \right)\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}} - \frac{2}{3}\frac{{\partial {u_i}}}{{\partial {x_j}}}{\delta _{ij}}} \right)} \right)\text{。}\end{array}$

(1)

式中：小标i和j分别为坐标方向； ${\rho _m}$ ，u和p分别为混合介质的密度，速度和压强； $\mu $ 和 ${\mu _t}$ 分别为混合介质的层流和湍流粘性系数。

1.2 湍流模型

计算中，湍流模型选择应用广泛的 $k - \varepsilon $ 两方程模型。标准 $k - \varepsilon $ 二方程模型是个半经验公式，主要是基于脉动动能k和能量耗散率ε，脉动动能k方程为：

$\frac{{\partial \rho k}}{{\partial t}} + \nabla \cdot (\rho \bar Uk) = \nabla \cdot [(\mu + \frac{{{\mu _{\rm{t}}}}}{{{\sigma _k}}})\nabla k] + {P_k} - \rho \varepsilon\text{，} $

(2)

耗散率ε方程为：

$\begin{split}&\displaystyle\frac{{\partial \rho \varepsilon }}{{\partial t}} + \nabla \cdot (\rho \bar U\varepsilon ) = \\& \displaystyle\nabla \cdot [(\mu + \frac{{{\mu _{\rm{t}}}}}{{{\sigma _\varepsilon }}})\nabla \varepsilon ] + \frac{\varepsilon }{k}({C_{\varepsilon 1}}{P_k} - {C_{\varepsilon 2}}\rho \varepsilon )\text{，}\end{split}$

(3)

式（2）与式（3）中：

${P_k} = {\mu _{\rm{t}}}\nabla \bar U \cdot (\nabla \bar U + \nabla {\bar U^{\rm{T}}}) - \frac{2}{3}\nabla \cdot \bar U(3{\mu _{\rm{t}}}\nabla \cdot \bar U + \rho k)\text{。}$

(4)

在 $k - \varepsilon $ 模型中，涡粘度 ${\mu _{\rm{t}}} = {C_\mu }\rho \frac{{{k^2}}}{\varepsilon }$ ，k表示湍动能，ε表示能量耗散率。式（2）和式（3）中的常数取值见表1。

有研究人员认为 $k - \varepsilon $ 空化区黏性耗散过强^[10]，因此普遍引入涡黏性系数的密度修正函数^[11]：

${\mu _t} = f\left( \rho \right){C_\mu }\frac{{{k^2}}}{\varepsilon }\text{，}$

(5)

式中：

$f\left( \rho \right){\rm{ = }}{\rho _v} + \frac{{{{\left( {{\rho _m} - {\rho _v}} \right)}^n}}}{{{{\left( {{\rho _l} - {\rho _v}} \right)}^{n - 1}}}},n \gg 1\text{。}$

式中：n一般取值范围为7~15，本文中n=10。

1.3 空化模型

Zwart空化模型是在Kubota^[12]空化模型的基础上做了一定的修正，其描述液-汽质量传递率的表达式如下：

${\dot m^{\rm{ - }}} = {F_{{\rm{vap}}}}\frac{{3{\alpha _{nuc}}\left( {1 - {\alpha _v}} \right){\rho _v}}}{{{R_{nuc}}}}\sqrt {\frac{2}{3}\frac{{\left| {{p_v} - p} \right|}}{{{\rho _l}}}}\text{，} $

(6)

${\dot m^{\rm{ + }}} = {F_{{\rm{cond}}}}\frac{{3{\alpha _v}{\rho _v}}}{{{R_{nuc}}}}\sqrt {\frac{2}{3}\frac{{\left| {{p_v} - p} \right|}}{{{\rho _l}}}}\text{，} $

(7)

式中： ${\dot m^{\rm{ - }}}$ 为蒸发过程； ${\dot m^{\rm{ + }}}$ 为冷凝过程； ${R_{nuc}}{\rm{ = }}1$ μm，为液体内简化的汽核半径； $ {r_{nuc}}{\rm{ = }}$ $5{\rm E}{\rm{ - }}4$ ，为单位体积液体内汽核的体积分数；下标v表示水蒸气；l表示液体水。

F_vap和F_cond为有关相变的经验系数（相变系数），分别代表蒸发系数和冷凝系数，默认的经验值为50, 0.01^{[3, 13]}。本文主要研究对于不同头型空化问题，这2个经验系数的设置规范。

另外，研究表明湍动能对空化也有重要影响，目前普遍采用以下方式对汽化压强修正^[2]：

${p_{turb}} = 0.39\rho k\text{，}$

(8)

汽化压强：

${p_v} = \left( {{p_{sat}} + \frac{{{p_{turb}}}}{2}} \right)\text{。}$

(9)

式中： ${p_{sat}}$ ，k分别为当地饱和蒸汽压强和流场的当地湍动能。

2 计算模型与参数

采用三维结构网格划分不同头型的圆柱模型所在流场，以半球头圆柱为例，其网格模型如图1所示。其中，边界层首层厚度保证y+=50左右。计算中采用速度入口，流速U=10 m/s；压力出口，相对压力0 Pa。

流场中无穷远处的参考压力根据流速和不同空化数确定。空化数计算公式为：

$\displaystyle\sigma = \frac{{{p_\infty } - {p_{sat}}}}{{\frac{1}{2}{\rho _l}{U^2}}}\text{。}$

(10)

3 计算结果

表面压力系数计算方法为：

$\displaystyle{p_{coe}} = \frac{{p - {p_\infty }}}{{\frac{1}{2}{\rho _l}{U^2}}}\text{。}$

(11)

Hunter等^[14]采用试验方法研究了不同头型在不同空化数条件下模型的表面压力系数分布情况，这类头型可以分为两类：流线型和非流线型。流线型头型的圆柱表面光顺，流场不易产生分离涡，如半球头圆柱、尖顶圆角圆柱等。非流线型头型会使流场产生强烈的脉动和流动分离，如平头圆柱、不同锥角的锥头圆柱等。本文研究了不同相变系数对这两类头型空化计算结果的影响。

3.1 流线型头型的空化计算

采用不同的相变系数计算空化数 $\sigma {\rm{ = }}0.4$ 时，半球头圆柱表面压力系数分布和蒸汽体积分数云图及空泡形态，如图2所示，其中表面压力系数试验值来自于文献[14]，以蒸汽体积分数等于0.5的等值面作为空泡表面。

当发生空化时，空泡内压力为饱和蒸汽压p_sat，因此空化区内的压力系数与空化数互为相反数，也是半球头圆柱体表面的最低压力系数。图2（a）中未显示出空泡表面，说明流场内部气相体积分数都在0.5以下，蒸发速率较慢。从表面压力系数分布看，最低压力低于饱和蒸汽压，也说明蒸发速率不够。主要原因是半球头圆柱表面光顺，周围流场相对稳定，脉动和涡特性不明显，不易产生空化现象。因此需要提高蒸发速率，使头部周围在低压环境下快速形成空化区，目标是使压力提高到饱和蒸汽压。保持F_cond不变，增大F_vap，计算了如图2（b）~图2（f）所示的算例。

从图2（b）~图2（f）中可以看出，当蒸发系数F_vap达到5 000时，表面压力系数最小值达到空化区内饱和蒸汽压产生的压力系数值，半球头圆柱头部圆周出现稳定空泡。然而圆柱表面压力达到最低值后呈缓慢趋势增加至零压，未出现正压力的尖峰突变。这种现象可以理解为冷凝速率较低，回射流不及时，空泡呈较慢的速度溃灭消失。因此需提高冷凝系数，验证冷凝系数对空化现象的影响，计算图2（g）~图2（k）所示的算例。

从图2（g）~图2（k）可以看出，当冷凝系数F_cond达到0.4后，表面压力系数分布已经很接近试验值，出现正压力峰值，圆柱头部有一个主空泡，主空泡后出现次生空泡脱落现象，与局部空化的试验现象相符。当冷凝系数F_cond达到0.7～1时，表面压力系数与试验值基本吻合。因此，对于半球头圆柱，当空化数 $\sigma {\rm{ = }}0.4$ 时，在数值计算过程中F_vap和F_cond宜分别设定为5 000，0.7～1。为验证这样的参数设置方案是否满足其他空化数条件下的计算精度要求，计算如图3所示的算例。

在一定空化数下，计算得到的表面压力系数出现一段较平稳的最低值，这个压力系数值与空化数互为相反数，代表这段区域压强为饱和蒸汽压，即出现了空泡，这段压力系数最低值的长度范围在一定程度上表明了空泡的长度尺度。由图3可以看出，将蒸发系数设定为5 000，冷凝系数设定为0.7，1时，不同空化数下表面压力系数的分布与试验值吻合的很好，间接地验证了数值计算得到的空泡尺度与试验结果一致，只有当空化数为0.2时，在空泡溃灭区压力分布与试验值相比有些异常，但总体来看压力变化趋势与试验值吻合，这样的对比结果验证了计算半球头圆柱空化问题时相变系数设置的合理性。

为了验证以上结论对其他流线型头型的空化问题是否适用，基于以上设置规范计算了尖顶圆角圆柱在不同空化数小的压力系数分布和空泡形态如图4和图5所示。

图4中不同空化数下，数值计算得到的表面压力系数与试验值变化趋势吻合，证明计算方法的准确性和相变系数设置的合理性。图5为不同空化数下空泡的形态变化。空化数为0.4时空泡初生，当空化数达到0.24时，空泡特征为一个主空泡后夹杂即将脱落的次生空泡。

以上计算结果表明，对于流线头型的空化问题，这样的相变系数设定原则合理。

3.2 非流线型圆柱的空化计算

以平头圆柱作为计算对象，计算其在空化数为0.4时不同相变系数下的压力分布情况和纵剖面蒸汽体积分数分布云图，如图6所示。

从图6可以看出，相变系数的改变对于平头圆柱表面压力系数分布影响不大。蒸发系数和冷凝系数的取值影响了流域内蒸汽体积分数的分布，并且流域内的蒸汽体积分数在平头圆柱周围的分布不均匀，具有强烈的不规则性。由于在获得表面压力系数分布曲线时是取纵剖面与圆柱表面交线的其中一侧，因此空化的不规则性在一定程度上会影响表面压力系数的分布情况。但是综合来看，不同相变系数下计算的压力分布结果与试验值变化趋势基本一致。

取F_vap=5000 ，F_cond=0.7（流线头型空化计算建议值）和F_vap=50，F_cond=0.01（空化模型默认值）2组相变系数计算平头圆柱在不同空化数下的表面压力分布如图7所示。

图中显示，不同相变系数下计算得到的表面压力系数变化趋势符合试验现象，但是与试验值存在一定的偏差，在空泡溃灭区压力上升段，计算值相对试验值有提前的现象。尤其当空化数较大时，试验中平头圆柱头部周围未达到空化压力却已经出现空化，这主要是由于平头圆柱头部周围产生严重的流动分离导致过早的出现空化现象，当空化数较小时（ $\sigma \leqslant 0.3$ ）最小压力才接近饱和蒸汽压^[13]。而满足较小的空化数条件时，表面压力系数分布对相变系数并不是很敏感，保证蒸发系数和冷凝系数的数量级关系即可模拟空化现象。

以上计算说明，对于非流线型圆柱的局部空化计算，只有当空化数小于0.3时，数值计算结果才能更接近试验结果。令F_vap=5 000，F_cond=1计算90°锥头圆柱在不同空化数条件下的空化现象，如图8所示。当空化数 $\sigma {\rm{ = }}0.3$ 时，计算结果的空化区长度与试验结果整体变化趋势一致，不同的是在空化溃灭区压力变化存在一定差异。当空化数 $\sigma {\rm{ = }}0.6$ 时，同样存在流动分离导致空化过早出现的问题，而数值计算中的空化模型无法捕捉这一现象，因此在空化初生区域数值计算得到的空化区压力系数比试验结果小。

4 结　语

采用基于CFD的数值计算方法结合Zwart空化模型求解不同头型圆柱在不同空化数条件下的表面压力系数，并与试验数据对比，可以得到以下结论：

1）流线型头型周围流场稳定，不易产生空化，增大相变系数才能保证计算结果接近试验结果，建议值为F_vap=5 000，F_cond=0.7～1.

2）非流线型头型周围流场流动分离和脉动强烈，在计算空化问题时对相变系数的变化不敏感，当空化数较小（ $\sigma \leqslant 0.3$ ）时，保证蒸发系数和冷凝系数的数量级关系即可较准确模拟空化现象。

3）对于非流线型头型，当空化数较大时（ $\sigma > $ $0.3 \sim 0.4$ ），流动分离会导致空化现象过早发生，从而使初生空泡区压力高于饱和蒸汽压，空化模型在捕捉这一特征方面略显不足。

后续工作将基于现有空化模型结构，考虑流场内流动分离强度因素，对汽化压强再次修正，使其能够适应非流线头型在较大空化数条件下空化问题的求解。