近自由面的超空泡流动的数值模拟研究

，，

(浙江理工大学机械与自动控制学院，杭州 310018)

0 引 言

超空泡现象在水中高速鱼雷、潜水艇发射战略战术导弹和反鱼雷鱼雷武器系统中有着重要的作用，水中高速运动物体在超空泡状态下运动时，可以减少阻力达90%。高速鱼雷作为对舰攻击的杀手锏武器，暴露出其在空气和水中连续机动隐身的性能，这也必然涉及超空泡流动。另外，空泡的溃灭又会给机器带来损伤，如产生水泵和水轮机的气蚀问题。因此，详细研究超空泡的流动机理，具有重要的实际意义，且一直是国内外的热门研究课题之一[1-2]。Chen等[3]通过求解RANS方程，模拟了超空泡在粘性浅水区域中的各种行为，包括靠近自由面处空泡和自由面的变形，并与实验数据进行了对比。超空泡产生时的空化数一般小于0.1，实现这个条件有三种途径[5]：a) 保证航行体的速度大于45 m/s；b) 降低空泡周围流场压力；c) 采用人工通气的方法增加空泡内压力。水下航行体在高速运行时周围液体汽化，产生包裹航行体的空化泡，即自然超空泡，此时运动体在空泡内运行。在深水位时，由于空泡形态受空腔内外压差及自由面的干扰较小，所以空泡形态才可以看成一个近似的椭球体，但是在接近水面时，自由面的影响将不会被忽略。自由面的波动将会引起空泡表面压力的分布，进而影响航行体的运动轨迹。

近年来，在国防科技技术和海洋领域以及自然科学领域方面，对超空泡的研究一直在不断进行。施红辉等[4-8]设计了最高速度可达100～120 m/s的模拟高速射弹实验装置，开展高速弹体水平入水实验，同时也进行了相关数值模拟，获得水平超空泡流动中自然超空泡的形态特性和一些水动力学特性。施红辉等[6]采用高速照相机记录了在280、34、18 mm三种不同水位下的弹体高速运动的过程，分析了水深对速度衰减的影响。施红辉等[7]通过大量实验对浅水区域中高速射弹进行研究，获得了近自由面处空泡形态的变化，并与半经验公式进行了对比。陈波等[8]通过采用Fluent软件模拟研究了浅水区时自由面对超空泡形状的影响，以及深水域时空化数对空泡形态的影响。马庆鹏等[9]对不同头型运动体高速入水的空泡形态进行了数值模拟研究，得到不同头型条件下高速入水运动参数及空泡形态发展规律和流场的速度分布规律。朱棒棒等[10]采用Fluent软件对不同头型物体高速入水过程进行了数值模拟，分析了入水空泡夹带空气的原因。黄海龙等[11]通过采用结构化网格对圆盘空化器进行了三维数值模拟，解决了带攻角圆盘空化器形成超空泡的问题，并分析了带攻角圆盘空化器生成超空泡的一些特征。Truscott等[12]对入水问题研究现状进行了综述，主要包括实验、理论和数值等方面的研究。

虽然目前国内外超空泡现象的研究已经有了许多研究成果，但是很少有文献考虑在浅水区自由面对超空泡流动及其诱导的波浪的影响。本文利用Fluent流体仿真软件，对物体在自由面附近的运动过程进行了数值模拟，将绘制的不同水深下的空泡直径的变化曲线与Logvinovich和Savchenko等[13-14]推导出的经验公式拟合的曲线作相关验证，进一步给出无量纲波浪高度随时间变化情况，通过对比不同工况下空泡形态及自由面的波动情况，分析自由面对超空泡流动的影响，给出先导波浪随时间的变化过程。模拟结果可为高速射弹超空泡的形态特性和流体动力学特性研究提供基础和参考。

1 数值方法

1.1 控制方程

本文为了解气、液、汽三相之间的相互作用，对近自由面超空泡流场进行了数值模拟。其中基本控制方程包括[8]：连续性方程和动量方程，选择VOF多相流模型和标准k-ε湍流模型；设置气、液、汽三相；采用Schnerr-Sauer空化模型。

连续性方程为：

(1)

动量方程为：

(2)

其中：μ表示混合物的动力粘度；ρ表示混合物的密度；SM表示自定义源；ui表示混合物的速度,i和j表示各分量的方向。

标准k-ε湍流模型基本形式为：

Gk+Gb-ρε-YM+Sk

(3)

(4)

其中：Gb表示浮力湍动能；Gk表示速度梯度湍动能；YM表示可压缩湍动能耗散率；C1ε;C2ε;C3ε为常数。

VOF模型的体积分数方程的一般形式为[8]：

(5)

Schnerr-Sauer空化模型的蒸气传输方程为：

(6)

1.2 数值模拟方法

对于二维流场来说，空泡形态可以用弹体轴线截面处空泡轮廓来描述,陈波等[8]采用的二维数值模拟方法，其模拟结果验证了二维流场的合理性，因此，本文也采用二维数值模拟方法对超空泡流场进行计算。本文给出了数值计算的物体模型，如图1所示。数值模拟中所使用的物体头型统一采用平头型，弹身为圆柱体[8]。计算区域及其网格划分，如图2所示。计算区域所取的是1900 mm×950 mm的矩形区域，采用“H”型进行网格加密处理。如图2所示，计算区域的左端设置成空气和水速度入口，右端设置成空气和水压力出口，上端设置为symmetry，下端和弹体周边均为壁面。压力耦合和速度耦合均采用SIMPLE算法，而梯度的求解则采用Least squares cell based(最小二乘法)，压力的求解则采用PRESTO!算法，至于动量离散则采用Second order upwind(二阶迎风格式)。

图1 弹体模型尺寸

图2 计算区域及网格划分

2 网格无关性验证

图3 数值模拟结果中波浪高度H和水深h以 及超空泡的长度Lc和直径Dc的确定

为了验证网格对数值模拟的影响，针对水深26 mm的工况，采用对自由面及弹体周围进行网格细化，得到相应的四个工况case 1、case 2、case 3、case 4，网格数分别是525832、550792、600712、700552。通过Fluent软件及后处理得到这四种工况的水相图，如图4(a)-(d)所示。

图4 不同网格数的水相图(v=60 m/s; t=6.0 ms)

将数值模拟结果所得的水相图导入到AutoCAD软件中，经过多次测量获得波浪高度的数据，将这些数据导入Origin软件中，输出无量纲的波浪高度的曲线图，如图5所示。通过对不同网格数的数值模拟，从图5中可以得出，在相同的水深等条件下，不同网格数对波浪高度的影响；对比四种工况，case 1和case 2有相同的走势；case 3和case 4有相同的走势，从而得证网格无关性。但是考虑到实际计算周期等其他条件的限制，又要为了使计算精度较高，同时又能更好地捕捉自由液面波动，对子弹周围和自由面附近做了网格加密处理，选取网格数为600712万，时间步长为2×10-6s。

图5 不同网格数下无量纲波浪 高度随时间的变化曲线

3 数值模拟结果与分析

3.1 水深h的影响

文中计算了四个水深(即h=26、106、206、406 mm)下的超空泡形状及自由面随时间的变化过程。图6是四个水深下的水相图。如图6所示,深水处的超空泡的外形轮廓是一个上下基本对称的椭球形(见图6(c)-(d))。然而，近自由面的工况下，自然超空泡为上表面形状出现上凸状态，下表面略显平坦(见图6(a))。

图7(a)-(d)给出了在四个水深下的，空泡上下壁面的流线；图8给出了了在四个水深下的空泡外形轮廓。比较图7和图8可知，是空泡形状引起了自由面的抬升，Wang等[16]也提出了类似的观点。当然，本文计算也指出，当水深大于206 mm后，空泡轮廓对自由面的抬起作用大大减小。

图7 v=60 m/s，t=6 ms时，不同水深下自由面及空泡周围的典型流线

图8 不同水深下，弹体上下两条流线痕迹形状(空泡外形轮廓)的对比曲线(v=60 m/s; t=6 ms)

从图6(a)可知，数值模拟算出了超空泡在自由面上拉出了一个倾斜的先导波浪，波浪的前端位置超过了物体或超空泡的前端位置，这是由于快速运动的物体表面形成的超空泡对浅水层挤压造成的。施红辉等[7]已经通过实验证实这一结论。另外，从图6(a)中的计算结果可知，空泡右上部分的水层已经很薄，预示在实际情况中空泡很可能破裂，而空气进入超空泡中并扩张其尺寸。这一点在实验中也得到了证实[4]。在图3中给出了自由面上的起波位置的定义，即波前距离S。图9(a)-(c)给出了S的计算结果，它们大致呈现这样的趋势：随着水深的增加，自由面上的起波位置距离物体头部越远；随着速度的增加，自由面上的起波位置距离物体头部越近。

图9 不同水深h和速度v下先导波浪随时间的变化曲线

自由面的变化过程用自由面波动引起的最大高度，即波浪高度H(见图3)来描述。图10给出了波浪高度随时间的变化，由图中可知：在深水位时，自由面的波动不是太明显，而在浅水区时，自由面的波动很大。针对h=26、106、206、406 mm四个水深，分别模拟了初速度为100 m/s和200 m/s的情况，进一步分析速度对自由面的影响，给出自由面引起的波浪高度的曲线图，如图11(a)-(d)所示。在图11中可以得出：在相同水深的条件下，随着速度的增加，自由面的波动较大。所以，通过模拟可以得出：自由面的变化不仅要受到水深的影响，而且还要受到物体速度的影响。由于水的声速为1500 m/s，当物体运动速度为200 m/s时，马赫数只有0.13，属于不可压流动。又因为本文采用的理想流体，所以这三种介质的压缩性不用考虑。

图10 在v=60 m/s时不同水深h下无量 纲波浪高度随时间变化曲线

图11 不同速度v下无量纲波浪高度随时间变化曲线

3.2 与半经验公式的比较

(7)

σv为空化数[13]，即

(8)

其中：P∞为周围压力；Pc为空泡内部压力；v∞为运动体的速度；ρ为流体的密度。

(9)

其中：cx0=0.82；k=0.9～1.0。该公式适用于空化数σ的范围为0～0.25。

而根据乌克兰流体研究所的Savchenko[14]提出了适用于空化数范围0.012～0.057的经验公式:

(10)

最后，根据Logvinovich空泡截面独立膨胀原理推导的空泡外形轮廓公式:

(11)

可以根据公式(8)-(11)绘制V=60 m/s,Dn=12 mm时的超空泡外形轮廓图。根据数值计算的结果，选取t=6 ms,h=26、106、206、406 mm时的数值模拟的空泡图片，导入AutoCAD软件中进行测量，得到相关空泡直径的数据，再进行无量纲化处理，在Origin软件中绘制超空泡的轮廓曲线如图12所示。

图12 不同水深h下的模拟空泡外形轮廓与Logvinovich 公式及Savchenko公式对比曲线(v=60 m/s; t=6 ms)

由图12可以得出，在浅水位h=26 mm时，计算出最大空泡直径要明显大于Logvinovich和Savchenko的半经验公式的值。这是由于当水平运动超空泡经过浅水区时，自由面将向上抬起，空泡形状上下不对称；而且上部薄水层受空泡拉动，不再满足Logvinovich的空泡独立膨胀原理。然而从数值量级上，计算结果与半经验公式的是一致的。

4 结 论

本文选择平头射弹模型，计算了水深对超空泡形状的影响。弹体在h=26、106、206、406 mm四种不同水深下高速运行时，随着水深的增加，空泡对称性越来越好。超空泡运动诱发自由面上产生波浪，其最大高度随着水深的增加而减小；随着物体速度的增加而增加。数值模拟还计算出了在自由面上的波浪，是一个先导波浪，即波浪的前端位置超过了物体或超空泡的前端位置，这与施红辉等[7]的实验结果相吻合。

本文首次分析了先导波的起始距离以及水深对它的影响。在浅水区时，自由面的起波位置靠近空泡前端，速度越大越靠近；在深水时，自由面的起波位置远离空泡前端，趋于无穷远处，且空泡形状更接近一个椭球体，可为高速射弹超空泡的形态特性和流体动力学特性研究提供参考。

参考文献：

[1] Saranjam B. Experimental and numerical investigation of an unsteady supercavitating moving body[J]. Ocean Engineering,2013,59(2):9-14.

[2] 孙士明,陈伟政,颜开.通气超空泡泄气机理研究[J].船舶力学,2014,18(5):492-498.

[3] Chen X, Lu C J. Properties of natural cavitation flows around a 2-D wedge in shallow water[J]. Journal of Hydrodynamics Ser B,2011,23(6):730-736.

[4] 施红辉,周杨洁,贾会霞,等.水深和弹体长径比对超空泡弹体阻力系数及空泡形状影响的实验研究[J].兵工学报,2016,37(11):2029-2036.

[5] 施红辉,周素云,张晓萍,等.水下超空泡流体机械的机理和技术研究综述[C]//第十五届全国激波与激波管学术会议论文集.杭州,2012:554-564.

[6] 施红辉,胡青青,胡俊辉,等.水平运动超空泡在近自由面区域的水动力学特性研究[C]//第九届全国实验流体力学学术会议论文集.杭州,2013:357-365.

[7] 施红辉,罗喜胜.可压缩性和高速多相流动[M].合肥:中国科学技术大学出版社,2014:228-231.

[8] 陈波,胡青青,施红辉,等.离自由面不同水深下水平运动超空泡的数值模拟研究[J].浙江理工大学学报,2015,33(3):375-381.

[9] 马庆鹏,魏英杰,王聪,等.不同头型运动体高速入水空泡数值模拟[J].哈尔滨工业大学学报,2014,46(11):24-29.

[10] 朱棒棒,施红辉,侯健,等.高速射弹入水时空气携带量的数值模拟[J].浙江理工大学学报,2017,37(3):402-408.

[11] 黄海龙,黄文虎,周峰.圆盘空化器超空泡形态三维数值模拟研究[J].兵工学报,2008,29(1):78-84.

[12] Truscott T T, Epps B P, Belden J. Water entry of projectiles[J]. Annual Review of Fluid Mechanics,2014,46(1):355-378.

[13] Logvinovich G V. Hydrodynamics of flow with free boundaries[M].施红辉,译.自由边界流动的水动力学[M].上海:上海交通大学出版社,2012,97-111.

[14] Savchenko Y N. Experimental investigation of supercavitating motion of bodies[C]// RTO-AVT and VKI. Brussels,2001:43-66.

[15] 易文俊,王中原,熊天红,等.水下射弹典型空化器的超空泡形态特性分析[J].弹道学报,2008,20(2):103-106.

[16] Wang Y, Wu X, Huang C, et al. Unsteady characteristics of cloud cavitating flow near the free surface around an axisymmetric projectile[J]. International Journal of Multiphase Flow,2016,85:48-56.