射弹跨声速入水初期阶段多相流场特性数值研究

陈 晨， 魏英杰， 王 聪， 毕殿方

(哈尔滨工业大学 航天学院，哈尔滨 150001)

入水过程是指运动体以一定的初始速度从空气中穿越自由液面并进入水中的过程。对于高速入水过程的研究一直以来都极具挑战，主要是因为在该过程中伴随着强烈的相变、湍动、流体可压缩等复杂的流动问题。另外当运动体的速度较高导致附近流体的压力低至饱和蒸汽压时，此时就会发生空化现象[1]。高速入水问题的一个重要应用就是超空泡射弹，如高效防御武器机载灭雷系统RAMIC(Rapid Airborne Mine Clearance System)，入水后弹体周围形成了超空泡，这使得射弹入水后仍能保持高速运动，提高了射弹的射程及毁伤能力，在打击鱼雷与水雷方面扮演着重要的角色。因此高速入水问题逐渐获得越来越多的关注。

对于入水问题的研究主要有两个方面：一是入水初期的运动体拍击问题，另一个是超空泡流的流体动力问题。对于入水拍击问题，较多的学者开展过诸如平底结构自由落体[2]、楔形体自由落体[3]、圆球倾斜入水[4]以及弹性圆柱体倾斜入水[5]等低速入水的研究。然而当速度较高时就需要考虑运动体与自由面之间的空气垫效应以及水的可压缩性的影响，尤其是对平头及钝头的运动体尤为重要。Karman[6]最先提出考虑了水的可压缩性的预测最大拍击压力的公式，p=ρcV，其中:ρ为水的密度;c为水中声速;V为入水速度。然而在之后的各种试验研究中发现，运动体实际受到的压力远小于应用Karman给出的预测公式得到的数值。这是由于该公式并没有考虑滞留在运动体与自由面之间的空气产生的空气垫效应。Verhagen[7]对滞留在平板与自由液面之间的空气的可压缩性进行了研究并给出了压力预测模型。Korobkin[8-9]应用声学近似法对钝体进入可压缩液体问题进行了研究，并给出了冲击波后的速度场及压力场分布。对于超空泡流的流体动力研究，Neaves等[10-11]应用时间导数预处理算法结合低扩散迎风方法对射弹以420 m/s垂直入水可压缩多相流进行模拟，并得到空泡形态、阻力系数及空泡外的压力分布。黄闯等[12-13]基于VOF(Volume of Fluid)模型对水下超空泡射弹进行了数值模拟，并针对水的可压缩性对阻力系数及空泡形态的影响进行了研究，结果表明在超声速范围内随着运动速度的提高，超空泡的尺寸快速减小。Shi等[14-15]应用小口径步枪以及光学观测技术对速度约为342 m/s的子弹开展了垂直入水试验，且观察到一系列如自由液面的运动、喷溅、水下空泡流动及射流等非定常流动现象。另外试验还发现，入水后的弹道偏移与入水深度有关，且子弹的冲击能量随着入水的深入而减少[16]。Truscott[17]针对入水弹道稳定性问题开展了不同形状子弹的小角度倾斜入水试验，试验表明具有大长细比、钝头的模型在水中运动时具有较好的弹道稳定性。

低速入水且忽略流体可压缩性的空泡流体动力学问题自十九世纪末[18]开始研究至今获得了大量的研究成果。但是高速入水试验由于受到试验设备以及测量条件的限制很少得到开展，少量的高速入水试验也是在军方的资助下完成的，因此公开的资料很少。另外，一般在对入水问题进行理论分析与数值仿真时都是忽略能量方程，并假定流体是不可压缩的。然而当初始入水速度达到声速甚至超声速时，空气、水及水蒸汽的密度均不再是定值，而是变成了随压力改变的变量，此时以上的假设也就不成立了。因此，由于问题本身的复杂性使得关于考虑流体可压缩性的高速入水问题的研究很少。

随着计算流体力学与计算机技术的发展，使应用数值模拟方法对高速入水多相可压缩空泡流问题进行研究得以实现。本文针对缩比的射弹模型跨声速垂直入水过程开展了二维轴对称数值模拟研究，考虑气、汽、液三相流体的可压缩性，对入水初期的流场发展进行分析，并研究流体可压缩性以及入水高度对入水初期流场结构的影响。

1 数值计算方法

数值模拟采用VOF多相流模型，控制方程组包括连续性方程、动量方程、湍流模型及空化模型。另外由于考虑了气、汽、液三相流体的可压缩性，为使方程组闭合所以还需添加能量方程及流体的状态方程。在建立数值计算方法后，应用文献中的试验数据及理论分析结果对本文的数值方法进行验证。

1.1 连续性方程与动量方程

本文基于VOF模型对雷诺时均方程组RANS(Reynolds-averaged Navier-Stokes)进行求解，其中连续性方程和动量方程为

(1)

(2)

式中：混合密度ρm=ρvαv+ρgαg+ρl(1-αv-αg)，下标v, g及l分别为水蒸汽、空气及水；u为速度；μm=μvαv+μgαg+μl(1-αv-αg)为混合相的动力黏度。

1.2 湍流模型

(3)

(4)

式中：Gb为由浮力产生的湍动能；YM为可压缩流中的脉动膨胀；Sk和Sε为源项；C1ε,C2ε,σk,σε及Cμ均为模型常数，它们的取值分别为C1ε=1.44，C2ε=1.92,σk=1.0,σε=1.3,Cμ=0.09；Gk为由平均速度梯度引起的湍动能，μt为湍动黏度，两者的表达式为

(5)

(6)

1.3 空化模型

当运动体以高速运动时会引起空化现象，因此需要引入空化模型用于描述水相与水蒸汽相之间的质量转换，本文选用Zwart-Gerber-Belamrik空化模型[19]用于计算，其表达式为

(7)

式中：气泡半径RB=1×10-6m；汽化核心体积分析αnuc=5×10-4；汽化系数Fvap=50；凝结系数Fcond=0.001。

1.4 能量方程

由于考虑了流体的可压缩性，因此还需对能量方程进行求解，能量方程的表达式为

(8)

式中：keff为有效传热率；认为能量E为平均质量的变量，其表达式为

(9)

式中：Ei为各相的比热容且共用同一温度值；Sh为由空化产生的源项，其表达式为

(10)

式中：Lev为汽化潜热。

1.5 状态方程

此外还需要添加流体的状态方程对能量方程进行封闭。对于水相，Saurel等[20]在标准状态方程的基础上考虑了温度的影响，同时认为参考压力是与饱和状态有关的变量，修改后的水的状态方程为

(11)

饱和压力与饱和密度是温度的函数，且函数中的Oldenbourg系数见表1。

(12)

1+b1θ1/3+b2θ2/3+b3θ5/3+b4θ16/3+b5θ43/3+b6θ110/3

(13)

式中：θ=1-T/Tc；水相的临界条件Pc=2.264×107，Tc=647.14 K，ρc=332 kg/m3。

表1 Oldenbourg系数

对于空气相及水蒸汽相，将其视为理想气体，理想气体的状态方程为

p=ρgRT

(14)

式中：R为气体常数；ρg为气体密度。

2 数值方法验证

对于跨声速入水问题的研究在公开文献中很少，为了验证本文数值方法的正确性，分别应用文献[21]中得到的试验结果以及Savchenko等研究中提供的计算公式与数值计算结果进行对比验证。

2.1 阻力系数验证

文献中试验装置示意图如图1所示，主要由火炮、水池、测速仪及防护板组成，当试验模型运动至试验区域内测速仪即触发，两组测速仪分别位于起始位置及距起始位置20 m处，防护板置于50 m处。试验模型见图2，试验结果见表2。

图1 试验系统示意图Fig.1 Experimental system

图2 试验模型Fig.2 Photograph of experimental model

试验编号发射速度/(m·s-1)20 m处速度/(m·s-1)平均速度/(m·s-1)阻力系数1900.45738.35819.40.012 8552508.63417.13462.880.012 864

图3为数值模拟结果与试验结果的阻力系数对比

图3 阻力系数对比Fig.3 Comparison of drag coefficients

图，数值模拟得到的阻力系数在入水前较小，在触水时刻阻力系数骤然增大，之后产生稳定空泡包裹试验模型从而使阻力系数趋于平稳；试验所得结果即为试验模型入水后稳定航行阶段的阻力系数。由图可知，考虑到试验中各项影响误差，数值模拟结果与试验结果吻合较好。

2.2 空泡形态验证

本文用于验证空泡形态以及后续研究入水多相流场的计算模型，如图4所示，该模型参考文献[22]进行了缩比与简化。计算域为二维轴对称模型，如图5所示，其中D为计算模型圆柱段直径。坐标原点取在模型头部中心点，+X为重力方向，Y轴为水平方向，模型头部距离自由液面10D，计算域半径为200D。图6为局部网格示意图，为保证空泡边界的计算精度，在模型周围应用非结构网格进行了加密处理，加密区之外应用结构网格沿着轴向与径向对网格进行渐疏处理。初始入水速度为400 m/s。

图4 计算用模型Fig.4 Photograph of computation model

图5 计算域示意图Fig.5 Computational domain

图6 局部网格示意图Fig.6 Local mesh

由于试验得到的数据有限，因此本文采用近似计算方法与基于细长体理论的方法再对数值模拟得到的空泡形态进行对比验证。

(1)近似计算方法

空泡前端具有较大的曲率并且仅由空化器形状决定而与空化数无关，因此空泡头部形态“1/3”规律[23]对于平头空化器是适用的，即

(15)

式中：R(x)为空泡半径；Rn为空化器半径。

Savchenko等[24]通过大量的速度范围在50～1 400 m/s内的试验结果，给出了适用于圆盘空化器的预测空泡形态的半经验公式

(16)

(2)细长体理论

另一个用于预测空泡形态的方法是基于细长体理论的针对圆盘空化器[25]的方法，其计算公式为

(17)

(18)

Cd=0.82C*(1+σ),Ma<1

(19)

Cd=0.82C*+σ,Ma>1

(20)

图7为射弹入水深度为L～4L时的空泡形态对比图，其中L为模型长度。由图7可知数值结果与理论结果吻合程度较高。综上，我们可以认为本文的数值模拟方法适用于射弹跨声速入水多相流场特性的研究。

图7 入水空泡形状对比Fig.7 Comparison of cavity shapes

3 结果分析

针对图4中的计算模型，开展初始入水速度为400 m/s的射弹垂直入水过程数值模拟，模拟过程包括射弹从空气中运动直至接触水面，然后历经撞击阶段、流动形成阶段与开空泡阶段。之后依据计算结果，分别对入水初期的流场发展、流体可压缩性以及入水高度对流场的影响进行分析。

3.1 入水初期流场发展分析

图8为来流马赫数为2时的圆柱绕流的头激波试验照片以及激波前后参数[26]，其中激波前的参数下标为1，激波后的参数下标为2。图9为计算得到的射弹在空气中运动时的压力云图，从图中可以看到在射弹的头部形成了与图8相似的弓形激波，在头部正前方为一小段正激波，但是在弓形激波其他位置激波斜角β随激波的向后延展而不断减小。另外射弹头部前方出现高压区，压力峰值在头部正前方，其他位置压力值沿着激波延展方向递减。这是由于对于跨声速运动，由正激波前后参数比的计算公式可知，激波后方即靠近射弹头部一侧的压力p2，密度ρ2以及温度T2均有所增加。低压区出现在射弹头部附近、锥段与柱段连接的肩部后方以及尾部边缘，这些位置均为流动分离处，由于此处流动的过膨胀使这些区域成为了低压区

图8 激波外形与前后参数Fig.8 Shock wave and parameters

图9 入水前压力云图Fig.9 Contour of pressure

。

为了研究在入水过程中流场结构的变化规律，本节给出了入水初期阶段的流线图，如图10所示，其中t=0为射弹运动至原未扰动自由液面高度的时刻。射弹头部具有两道激波，较前端的激波为弓形激波，后端的激波为再压缩激波。随着射弹进入水中，头部附近的液体获得射弹传递的动能而沿周向流动，此时空泡边缘的流线是以头部为原点四散开来。同时由于传播介质的改变，弓形激波的流线从空泡内指向水中时发生偏转，且随着射弹入水的深入，空泡内的弓形激波向后的弯曲程度增加。而再压缩激波则随着射弹入水的深入而不断后退，且波形趋于水平。另外由于空泡壁上的反射激波以及透射激波的共同作用下，空泡口喷溅位置附近产生了旋涡结构，旋涡结构将促进喷溅处不同流体的混合。

图10 入水初期阶段的流线图Fig.10 Streamlines of the fluid domain

图11给出了入水初期阶段不同时刻下空泡壁面的速度矢量图。当射弹入水后，头部周围的液体在得到动能后开始向四周流动，从而形成了与大气相连接的空泡，此时空泡壁面的速度方向由空泡指向水中，因此空泡将继续扩张，且扩张速度沿重力方向逐渐增大，速度最大值均出现在头部流动分离点位置，这是因为液体距离射弹头部越近获得的动能越大。射弹的运动速度随着入水的深入而衰减，则传递给流体的动能也随着减少，因此分离点处的速度也逐渐减小。自由液面位置的喷溅在射弹入水撞击阶段中获得的速度最大，之后在重力及空气阻力的作用下喷溅的速度逐渐减小。

图11 空泡壁面速度矢量云图Fig.11 Contour of vectors on the cavity wall

3.2 流体可压缩性的影响分析

为了研究流体可压缩性对入水初期阶段流场的影响，本节分别进行了考虑气、汽、液三相流体可压缩性以及忽略所有流体可压缩性的射弹入水数值模拟，两种计算工况的初始入水速度均为400 m/s，且其他设置相同。

射弹在空气中运动时的压力分布云图，如图12所示，由图12可知高压区主要分布在射弹的头部，且压力峰值位于轴线附近；而低压区主要集中在头部边缘后侧。其中图12(a)所示的高压区压力云图为向弹体后方延伸是弓形，且压力值沿射弹头部前方以及弹体后方递减，这是由于射弹的运动速度大于传播速度为声速的压力波，导致了压力波的多重叠加，从而在射弹头部前端形成了强压缩波阵面，即激波。图12(b)中射弹头部前方形成了一个球状的高压区且沿着周向扩散，流场中的压力值随着距射弹头部的距离增加而减小，由于忽略了流体的压缩性，此时的压力没有因为运动速度超过声速而发生叠加，也就不会有激波形成。

图12 在空气中运动时的压力云图Fig.12 Contour of pressure when projectile travels in the air

图13为射弹在运动至原未扰动自由液面高度时的相图。当射弹到达自由液面高度时，忽略流体可压缩性的射弹头部直接与水面接触并由于入水冲击溅起了水花。然而，当考虑流体可压缩性时，射弹头部并没有接触到水面，而是使头部以下的自由液面发生了凹陷，由于射弹的头部为平面，则在头部与自由液面之间的空气来不及逃逸出去，从而形成了一层空气垫，另外射弹头部存在激波，激波到达自由液面也造成了气液界面的RM(Richtmyer-Meshkov)不稳定性扰动，通过以上两者的共同作用射弹头部附近的自由液面发生了凹陷。

图13 射弹运动至原未扰动自由液面高度时刻水相云图Fig.13 Contour of water volume fraction when projectile moves to the undisturbed free surface

图14给出了射弹在运动至原未扰动自由液面高度时的流线图。从图14中可以看出流体可压缩性对射弹头部流线的影响，当没有考虑流体可压缩性时，头部与水面之间的流线是沿着自由液面向两侧延展，而考虑流体可压缩性时此处的流线是指向自由液面的，这也就说明此时的空气没有像图14(b)一样得以逃逸。

另外，由图14(a)可以看出流线在射弹的头部边缘、肩部和尾部出现旋涡及明显分离线，即由于射弹的入水速度大于声速从而使射弹周围产生了激波；而图14(b)中则由于没有考虑流体的可压缩性并不存在激波，仅在射弹的头部边缘、肩部及尾部产生旋涡。

图14 射弹运动至原未扰动自由液面高度时刻流线图Fig.14 Streamlines of the fluid domain when projectile moves to the undisturbed free surface

3.3 入水高度的影响分析

为研究射弹入水高度对流场结构的影响规律，本节分别进行了射弹头部距离自由液面2 mm,6 mm,10 mm,20 mm,30 mm以及40 mm高度的入水过程数值计算，初始入水速度均为400 m/s。

图15为射弹运动至原本未扰动自由液面高度时的空气相密度云图，其中H为入水高度。从图15中可以看出当入水高度为2 mm时，由于距离自由液面较近，射弹运动至液面前激波还没有得到充分发展，尽管此时密度的峰值较大，但也仅存在头部附近极小的部分区域，而在头部与液面之间的空气密度较小，即靠近液面附近的空气没有得到较大程度的压缩；入水高度为6 mm时，头部已产生弓形激波，但从波形对比可以看出其激波的发展程度远小于较大入水高度的激波，且靠近液面附近的高密度区较小，因此可以认为空气垫效应及激波对液面的作用时间较短；而后随着入水高度的增加，激波逐渐发展完全，并且高密度区分别位于头部附近与液面附近，因此可以认为在射弹运动至原自由液面高度时空气垫效应及激波造成的RM不稳定性对自由液面已产生一段时间的作用；当入水高度大于10 mm之后，激波的形状以及射弹头部前方的空气密度分布基本一致。

图15 射弹运动至原未扰动自由液面高度时刻空气相密度云图Fig.15 Contour of the air density when projectile moves to the undisturbed free surface

射弹运动至原未扰动自由液面高度时造成的液面凹陷深度随入水高度的变化曲线，如图16所示，其中h为凹陷深度。由图16可知，当入水高度为2 mm和6 mm时凹陷深度较大，之后凹陷深度随着入水高度的增加而减小，但当入水高度继续增加时凹陷深度则变化不大。这是由于当入水高度较小时空气垫效应以及激波造成的RM不稳定性对自由液面作用的时间较短，即液面刚进入扰动阶段。对于入水高度较高的工况，正如图15所示，液面已受到一段时间的扰动了，此时液面高度有所恢复，因此相对于2 mm和6 mm入水高度的工况其凹陷较浅。

图16 水面凹陷深度随入水高度变化曲线Fig.16 Free surface depression versus entry height

通过对头部中点的压力监测可以得到压力随射弹运动时间变化的曲线，如图17所示，其中t为运动时间，t<0表示运动体在空气中运动。从全局图中可以看出，射弹在空气中运动时其头部受到的压力相对较小，由于水的密度远大于空气的密度，因此造成在触水瞬间压力在极短的时间内迅速增大，之后随着空泡的形成压力开始下降并逐渐趋于平稳。由图17左上方的局部放大图可以看出，入水高度越大，触水时刻越早；当入水高度大于10 mm时，三种工况的触水时刻基本相同。这与图16中的规律相同，即水面凹陷越深需要运动的时间越长。由图17左下方的局部放大图可知，入水高度越大，压力下降得越早。结合两张局部放大图可知，入水高度越大，整个入水撞击阶段发生得越早。

图17 头部压力随时间变化曲线Fig.17 Pressure of the nose versus time

图18为射弹入水撞击阶段头部中点的压力最大值随入水高度的变化曲线，压力值随着入水高度的增加而减小。由图15的分析可知，入水高度越大，头部与液面之间的空气垫效应越明显，因此对射弹头部入水冲击的缓冲作用越大，相应的压力峰值越小。同样地，对于入水高度大于10 mm的工况，压力峰值变化不大。

图18 头部最大压力随入水高度变化曲线Fig.18 The maximum pressure versus entry height

图19给出了射弹触水后头部附近水的密度峰值随入水高度的变化曲线，其变化趋势由式(11)可知，当考虑流体的可压缩性时，水的密度是与压力有关的函数，因此水的密度峰值随入水高度的变化规律与图18基本一致。

图19 入水后水的密度最大值随入水高度变化曲线Fig.19 The maximum density of water versus entry height

4 结 论

本文针对射弹高速垂直入水问题，考虑气、汽、液三相流体的可压缩性，对入水撞击阶段、流动形成阶段及开空泡阶段开展了数值仿真，并对入水初期的流场结构进行了研究分析，具体结论如下：

(1)由于考虑了气、汽、液三相流体的可压缩性，射弹在以跨声速运动时头部会产生弓形脱体激波，且随着射弹浸入水中，弓形激波的激波斜角逐渐减小，再压缩激波不断后退。

(2)随着入水深度的增加，空泡不断扩张，且射弹头部附近的空泡扩张速度较大，越接近自由液面空泡扩张速度越小；空泡口喷溅的速度则随着入水深度的增加而减小。

(3)当考虑了各相流体的可压缩性时，射弹运动至原未扰动自由液面高度处并没有接触到水面，而是在空气垫效应及激波对液面RM不稳定的作用下使液面发生凹陷。

(4)当入水高度为2 mm与6 mm时，头部弓形激波在入水前发展不够充分，从而导致射弹运动至原未扰动自由液面高度处引起的液面凹陷较大，入水撞击阶段发生得较晚，且入水拍击压力较大；随着入水高度的增加，头部前端的弓形激波发展越充分，入水撞击阶段发生越早，空气垫效应越明显，从而使入水拍击压力越小；但当入水高度进一步增加时，入水高度对入水流场结构的影响很小。