王 聪，何春涛，权晓波，魏英杰

(1.哈尔滨工业大学航天学院，150001哈尔滨;2.中国运载火箭技术研究院总体设计部，100076北京)

物体高速入水时，在水面打开的1个开口空腔，称为入水空泡.入水空泡的生成、发展以及溃灭过程，涉及流固耦合、介质突变、自然空化以及强非线性等，是1个非常复杂的过程，尤其会对空投鱼雷、导弹、超空泡射弹等入水瞬间的流体动力特性和水下弹道的稳定性产生较大的影响.对入水问题的研究最早可以追溯到十九世纪末A.M.Worthington和R.S.Cole［1］所开展的关于液滴和球体等物体入水时引起的自由液面喷溅和入水空泡的研究等.1929年Von K［2］在不考虑重力、水的粘性力以及空气惯性力的基础之上完成了飞机水面降落冲击载荷的研究.上世纪50年代Albert May等［3-7］在美国海军的支持下开展了大量关于入水问题的基础研究，主要研究内容包括垂直和倾斜入水空泡、流体动力、入水撞击和入水弹道等.1968年J.G.Waugh［8］等人在美国水下海军作战中心的支持下基于弗劳德数和空气密度等参数的相似性，完成了不同头型的Mk25空投鱼雷缩比模型倾斜入水的实验研究.Albert May［7］和David Gilbarg等［9］提出将空气压强与标准大气压强的比值，或密度与标准大气压强下空气密度的比值作为入水空泡闭合前的空泡动力相似参数，并对不同压强和密度条件下的球体低速入水空泡形态进行了试验研究.研究表明空气压强对入水空泡的运动特性和空泡闭合方式等具有不可忽略的影响，因此空气压强常被作为最主要的入水空泡流动特性的相似参数之一.根据以上情况，本文将采用数值研究的方法对上述特性开展进一步的研究.随着测试设备和手段的发展，越来越多的研究人员对入水问题的机理开展进一步的研究.Jeffrey M.Aristoff等［10］利用PIV等对低密度的球体垂直入水过程中球体和空泡壁面速度变化规律进行了试验研究;Cyril Duez等［11］不同亲水性的球体垂直入水过程所产生的喷溅进行了对照试验研究;Hongmei Yan等［12］对低Re数(o(10))条件下球体垂直入水空泡动力学进行了试验研究.

综上所述，过去对入水问题的研究大多都是通过试验的方法，近年来人们(在大量宝贵的实验数据基础之上)开始发展入水空泡的数值计算方法［13-14］，由于数值计算方法大多是半经验型的，因此，对于一些复杂影响因素难以考虑周全.本文在RANS方程基础之上引入VOF多相流模型，并通过引入水和水蒸汽之间的输运关系，在考虑自由液面的同时，还考虑了由入水所引起的自然空化现象，对带150°锥角圆柱体匀速垂直入水空泡形态与空气域压强之间的变化规律进行了数值研究.在计算过程中采用动网格技术实现圆柱体入水过程的绝对运动与流场运动之间的耦合求解过程.本文假设空气域为理想气体，空气密度随空气域压强变化而发生相应的变化.分别对1/8～4倍标准大气压强(P0=101 325 Pa)条件下，150°锥角头型圆柱体以10～50 m/s速度匀速垂直入水的空泡进行了数值研究.

1 控制方程

1.1 连续性方程和动量方程

混合相连续性方程为

空气相连续性方程为

水蒸汽连续性方程为

混合相动量方程为

式中:ui表示混合物速度分量(m/s);ρm为混合物密度(kg/m3)，表达式如下:

μm是混合物动力粘性系数(m/s)，

下标m、l、g和v分别表示混合相、水、空气和水蒸汽;ag、av和1-ag-av分别为空气、水蒸汽和水的体积分数;Re、Rc表示水蒸汽产生和凝结过程的源项.

1.2 空化模型

高速入水物体在入水过程中产生空泡，仅有物体头部与水接触时，在物体头部边缘分离处会出现较大的压强梯度变化，当压强降低到饱和蒸汽压强之下时，可能会出现自然空化现象.自然空化是1个非常复杂的相变过程.本文通过式(3)中的源项来表示该过程(即Zwart-Gerber-Belamri空化模型)，表达式如下:

式中:RB=10-6m为气泡半径;anuc=5×10-4为不可凝结气体体积分数;Fvap、Fcond为经验常数，分别取50和0.001.

1.3 湍流模型

本文计算中采用标准k-ε湍流模型，湍动能k和湍能耗散率ε的输运方程分别为

式中:μt=ρmCμk2/ε为湍流粘性系数;σk=1.0和σε=1.3分别为k和ε的Prandtl数;C3ε=tan｜v/u｜，v为平行于重力场矢量的流速分量，u为垂直于重力场矢量的速度分量;Cμ、C1ε和C2ε为经验常数，分别取0.09、1.44，和1.92;Gk=为由平均速度梯度所引起的湍动能生成项.

2 数值求解

2.1 计算域和边界条件

本文计算模型采用150°锥角头型圆柱体，直径D=10 mm，长度L=40 mm.计算域设置和模型表面网格分布如图1所示，空气域长度为60D，水域为400D，径向为80D.文献［1］中的研究表明，轴对称物体垂直入水所产生的空泡具有轴对称性质.因此，本文数值计算均采用二维轴对称计算模型.使用四边形网格，网格数为324 750.同时还采用了107 000、168 800和515 000三种不同数量的网格进行对照计算，计算结果表明网格数较大的3种计算模型的空泡形态和空泡表面闭合时间等参数之间差值小于5%.综合考虑计算结果精度和计算量等因素，最后选取324 750数量的网格作为本文计算网格.

图1 计算域和圆柱表面网格分布

2.2 数值方法

本文采用有限体积法对控制方程进行空间离散，并应用PISO方法对动量方程和连续性方程联立求解建立速度场和压力场的耦合关系.计算过程中对流项的离散采用二阶迎风格式;水蒸汽相的离散采用QUICK格式;湍流输运方程的离散采用二阶迎风格式;对时间的离散采用一阶显示格式.计算时间步长的选取根据网格尺度、运动速度等因素综合确定.并通过动网格技术实现圆柱体的运动过程与计算域内的网格更新.由于本文计算采用结构化网格，网格动态更新过程中主要是被拉伸或被压缩变形，因此，采用动态层法进行网格更新.动态层法的核心是当网格被拉伸到大于设定的理想最大尺度((1+as)hideal)时，该网格分离成两层新的网格，当网格被压缩到小于设定的理想最小尺度(achideal)时网格与其相邻的网格合并成新的网格，即

引入动态网格后，相当于网格按一定速度作对流运动，在实际计算中应减掉网格的对流效应.因而引入动网格后控制体V对变量φ的统一形式的守恒方程可写为

式中:ρm为混合物密度;为速度矢量;g为动网格的运动速度;Γ为扩散系数;Sφ为标量φ的源项;∂V表示控制体积V的边界.

3 数值计算方法对比与验证

本文通过对Albert May等［7］提出的匀速垂直入水理想空泡模型进行对比计算，验证本文所采用的数值计算方法的可信性.为了便于分析，引入如图2所示的球体垂直入水空泡柱坐标剖面示意图.

图2 垂直入水空泡坐标

坐标原点建立在未扰动水平自由水面与球体对称轴交点处，径向半径与未扰动水平自由水面重合，竖直向下为轴向正方向.在该二维轴对称坐标系下，垂直入水早期空泡形态的理想空泡模型为［7］

式中CD为入水物体的阻力系数，d为入水物体的特征尺寸.该入水空泡模型主要是用于描述入水过程中，空泡闭合之前发展较为充分时期的外部轮廓.

根据以上理想空泡模型边界条件要求，对球体以30 m/s从空气域匀速垂直入水过程进行数值计算，空气域环境压强为101 325 Pa，空气、水和水蒸汽的密度分别为1.225、998.2和0.554 kg/m3.将以上计算结果选取入水空泡充分发展的某一时刻的空泡形态，与相同时刻下由式(12)理想空泡模型的拟合结果进行对比.阻力系数取CD=0.20［15］，并用球体直径作为特征尺寸分别对空泡直径Da和长度La进行无量纲化.数值计算与理论公式计算结果对比如图3所示.

图3 球体入水空泡轮廓计算结果与经验公式对比

由图3可以看出，数值计算结果与理想空泡模型拟合结果具有较好的一致性.在空泡轮廓两端具有一定的误差，产生该误差的主要原因是由于理想空泡模型本身忽略了自由表面处和分离点处的空泡轮廓的不规则性.通过以上对球体垂直入水空泡形态的计算方法的验证，说明了本文所采用的数值计算方法是合理可信的.并在此基础上展开空气域压强(或密度)变化对带150°锥角头型圆柱体匀速垂直入水空泡形态的影响规律的数值计算研究.

4 空气压强对空泡影响分析

以下对带150°锥角头型圆柱体分别以10、20、30、40和50 m/s的速度匀速垂直入水的空泡生成过程、自然空化以及不同空气压强对空泡闭合影响规律等进行数值计算，并对结果进行分析.

4.1 自然空化对空泡形态的影响分析

图4给出了圆柱体以30 m/s匀速垂直入水过程中不同时刻空泡的形态轮廓图.零时刻圆柱头部与自由水面接触.

由图4可以看出，随着入水时间的推移，圆柱体逐渐打开1个与空气相同的空腔，即入水空泡.空泡的生成过程实质上是1个圆柱体与水之间的能量转换过程.圆柱体撞击自由水面将动能传递给水体，并通过阻力作用持续地将动能传递给水体，动能将驱使水体运动，同时水体又受到由重力所产生的静压的阻碍作用.空泡壁面附近的水介质在动能和势能的共同作用下运动.由于圆柱体撞击自由液面时，在自由液面附近的水介质会产生1个绕圆柱体肩部的轴向速度分量，使得部分水介质以该速度相上运动，形成喷溅.并由于空泡壁面的运动，使得自由液面附近的水面上升，并超过了未扰动自由液面高度.在空泡入口处，由于喷溅的作用，使得入口逐渐变小以至于空泡表面闭合，阻止了外界空气进入到空泡内部，而由于惯性的作用空泡体积继续增加，使得空泡内部压强迅速降低，分布变得更加复杂.图5为标准大气压下，圆柱体以50 m/s匀速垂直入水空泡闭合瞬间肩部压力梯度云图.

图4 圆柱体以30 m/s匀速垂直入水时的空泡云图

图5 圆柱体以50 m/s入水圆柱肩部压力云图

由图5可以看出，由于头部绕流和肩部锐边使得流体在流过锐边时产生较大的顺压梯度.根据伯努利方程可知，当dP/ds=0时，压强降低到最低值.根据柯乃普［16］对空化的定义，如果压强值低于水的饱和蒸汽压强(本文取3 510 Pa)，并在时间和空间上持续一定长度，则会产生自然空化.由于此时空泡处于刚刚闭合或快要闭合状态，空泡内部压强接近空气压强，远远高于水的饱和蒸汽压强，其将阻碍空化的发展.上述压强梯度变化所对应的肩部自然空化的水蒸汽体积分数如图6所示.

图6 以50 m/s速度入水的圆柱肩部水蒸汽体积分数云图

由图6可以看出，只是在肩部后小范围内出现局部的自然空化，且水蒸汽的含量较低.通过对10、20、30和40 m/s速度垂直入水的计算结果的自然空化状态对比分析表明，随着入水速度的增加，空化区域呈增大趋势.若以10 m/s和20 m/s速度垂直入水，在圆柱体肩部几乎不发生空化.在10～50 m/s速度范围内，带150°锥角头型圆柱体匀速垂直入水，自然空化的范围较小，只是出现在肩部，对入水空泡形态的影响较小.

4.2 空气压强对空泡表面闭合影响分析

图7 给出了圆柱体在不同速度条件下匀速垂直入水时空泡无量纲闭合时间Ts/(D/V)与空气压强之间的关系曲线.

图7 入水空泡闭合时间与空气压强之间的关系曲线

其中，Ts为空泡表面闭合时间，D为圆柱体直径，V为圆柱体入水速度;P0为标准大气压强(本文取101325 Pa)，P为实际空气压强.计算压强分别为4P0、2P0、P0、(1/2)P0、(1/4)P0和(1/8)P0.由图7可以看出，不同速度条件下入水空泡无量纲闭合时间随压强的变化趋势基本一致，当空气压强高于1个标准大气压时，随着压强的升高，空泡在自由液面上闭合即表面闭合时间缩短，基本呈线性趋势变化;当压强低于1个标准大气压强时，随着压强的降低，表面闭合时间迅速增加，呈非线性趋势变化.同时，可以看出在相同压强条件下不同速度入水空泡表面闭合时间趋于常数，当空泡闭合方式由表面闭合转变为在自由液面之下某一点闭合即深闭合后，空泡闭合时间变化规律变得较为复杂.图中P/P0=1/4时，10 m/s速度入水空泡闭合方式由表面闭合转变成为深闭合.随着空气压强的降低，低速度入水空泡最先出现空泡闭合方式的转变，即由表面闭合转变成深闭合.图8给出了圆柱体以10 m/s速度分别在2P0、P0、(1/2)P0和(1/4)P0时，空泡第一次闭合瞬间空泡轮廓图.

图8 不同压强条件下入水空泡闭合瞬间轮廓图

由图8可以看出，在压强高于(1/2)P0时空泡在表面发生闭合，当压强为(1/4)P0时空泡在表面以下某点发生闭合.并且随着空气压强的降低空泡闭合时间延后，当降低到一定程度后空泡闭合方式发生转变，即由表面闭合转变为自由液面以下某一点的深闭合.根据本文对不同速度条件下圆柱体入水计算结果分析表明，随着速度的增加，空泡发生闭合方式转变的临界压强降低.同时可以看出，在空泡闭合时刻，入水点附近自由液面的高度随着空气压强的降低而增加，空泡直径随着压强的降低也呈增加趋势.图9为在不同空气压强条件下，不同速度垂直入水空泡闭合瞬间空泡最大直径变化规律.

由图9可以看出，在同一压强条件下，以不同速度入水，其入水空泡最大无量纲直径几乎趋于某一定值;同时，随着空气压强的降低入水空泡闭合瞬间最大无量纲直径呈递增趋势，在空泡发生闭合方式转变之前，空泡无量纲直径随压强变化的趋势较为稳定，当压强降低到一定程度后，空泡无量纲直径出现较大幅度的增加，并出现不稳定现象.

图9 空泡最大直径变化规律

5 结论

本文对带150°锥角头型圆柱体匀速垂直入水空泡形态进行了数值计算研究.计算工况为:入水大气压强从(1/8)P0～4P0，入水速度从10 m/s到50 m/s.通过对以上计算结果分析得出以下结论.

1)通过对不同速度和不同压强条件下入水自然空化的计算结果分析表明，在本文所计算入水速度范围内，在圆柱体肩部出现小范围内的局部自然空化，自然空化对整个空泡形态影响较小.

2)随着空气域压强的降低，空泡表面闭合无量纲时间呈增加趋势变化;当空气压强高于1个标准大气压强时，空泡无量纲闭合时间随压强的变化趋势较为稳定;当压强低于1个标准大气压强时，空泡无量纲闭合时间变化较为复杂，呈现出较强的非线性.

3)当空气域压强降低到一定程度时，空泡闭合方式将由表面闭合向深闭合转变.

4)空泡闭合瞬间最大无量纲直径在同一空气域压强条件下趋于某一稳定值;并随着压强的降低，出现一定的不稳定性.

［1］WORTHINGTON A M，COLE R S.A study of splashes［M］.New York:Longmans Green and Company，1908:13-119.

［2］VON K.The impact on seaplane floats during landing［R］.Washington:National Advisory Committee for Aeronautics，1929:1-8.

［3］MAY A.Review of water-entry theory and date［J］.Journal of hydrodynamics，1970，4(4):140-142.

［4］MAY A，WOODHULL J.C.Drag coefficients of steel spheres entering water vertically［J］.Journal of Applied Physics，1948，19(12):1109-1121.

［5］MAY A.Effect of surface condition of a sphere on its water-entry cavity［J］.Journal of Applied Physics，1951，22:1219-1222.

［6］MAY A.Vertical entry of missile into water［J］.Journal of Applied Physics，1952，23(12):1362-1372.

［7］MAY A.Water entry and the cavity-running behavior of missiles［R］.Virginiai Nowal Surface Weapons Center，1975:1-37.

［8］WAUGH J G.Water-entry pitch modeling［J］.Journal of hydrodynamics，1968，2(2):87-92.

［9］GILBARG D，ANDERSON R A.Influence of atmospheric pressure on the phenomena accompanying the entry of spheres into water［J］.Journal of applied physics，1948，19(2):127-139.

［10］ARISTOFF J M，TRUSCOTT T T，TECHET A H，et al.The water entry of decelerating spheres［J］.Physics of fluids，2010，22(3)，032102:1-7.

［11］DUEZ C，YBERT C，CLANET C，et al.Making a splash with water repellency［J］.Nature physics，2007，3:180-183.

［12］YAN H，LIU Y，JAKUB K，et al.Cavity dynamics in water entry at low Froude numbers［J］.Journal Fluid Mech，2009，641:441-461.

［13］陈学农，何友声.平头物体三维带空泡入水的数值模拟［J］.力学学报，1990，22(2):129-138.

［14］ARISTOFF J M，BUSH J W M.Water entry of small hydrophobic spheres［J］.Journal Fluid Mech，2009，619:45-78.

［15］李玉柱，苑明顺.流体力学［M］，北京:高等教育出版社，2003:173-179.

［16］柯乃普，戴利，哈密脱.空化与空蚀［M］.水利水电科学研究院，译.北京:水利出版社，1981:1-12.