马庆鹏，魏英杰，王聪，赵成功

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

运动体高速入水过程中，在其运动轨迹上将形成一个由空气和空化产生的水蒸气构成的空腔，称为入水空泡.随着运动体不断向下运动，入水空泡不断拉长，其径向也将呈现出先扩张、后收缩的趋势.与此同时，由于运动体速度较高，在撞击水面以后，其头部将承受较大的冲击载荷并随着速度的衰减不断降低.入水空泡的形态及入水冲击载荷都与运动体的尺寸、头型有着密切联系，对于高速入水问题的研究，多采用实验与理论方法相结合的方法开展.文献［1］以口径12.7 mm及14.5 mm穿甲弹为主要研究对象，开展了以800～1 070 m/s初始速度入水的实验，研究了高速入水空泡发展、溃灭的发展规律及流体载荷，并从能量守恒角度给出了预测空泡半径的公式.文献［2］在文献［1］基础上提出了空泡深闭合的理论分析方法，对初始速度为500～1 500 m/s的球体入水空泡发展过程开展了相关研究，并与基于LSDYNA的数值模拟结果进行对比，结果吻合较好.

文献［3］基于 Besant(1859)-Rayleigh(1917)空泡压力平衡理论，将空泡流场看做有势流，推导得到空泡半径随时间变化的方程.对于球体和圆柱体入水空泡发展的理论求解结果与实验结果进行对比，得到了较好的一致性.文献［4］将文献［5］提出的独立扩张原理对小型运动体倾斜入水过程的空泡形态进行了预测和分析，并将理论分析结果和实验结果进行对比，二者具有较高的一致性.文献［6-7］利用二级轻气炮、高速摄像机等实验设备，开展了初始入水速度在100～600 m/s范围内的小型运动体水平入水实验，对比分析了运动体入水速度、深度的变化规律，并基于文献［3］的理论分析模型与实验结果进行对比，获得了较好的一致性.随着计算机技术的发展，通过CFD(computational fluid dynamics)方法对入水问题开展数值模拟，得到较为精确的空泡形态及载荷数据.文献［8］基于雷诺时均的Navier-Stokes方程，并采用k-ε湍流模型及三维动网格技术，对航行体入水过程的空泡形态及流体动力开展数值模拟研究，分析了空泡形态的发展规律及阻力系数、升力系数、俯仰力矩系数等流体动力系数，为下一步入水弹道的研究提供了输入参数.文献［9-10］等通过求解雷诺平均的Navier-Stokes方程，针对小型运动体在低速(10～50 m/s)条件下的入水问题开展了数值模拟研究，分析了入水空泡发展规律，并研究了大气压强对入水空泡面闭合的影响规律.

目前国内对入水问题的数值模拟研究多基于低速条件下(＜50 m/s)，本文基于有限体积法，求解雷诺平均的Navier-Stokes方程，对带有锥头的圆柱体以500 m/s的较高速度垂直自由入水问题开展了数值模拟研究.针对带有3种不同角度锥头的圆柱体模型，分析了入水弹道、空泡半径及流场特性与头型的关系，并通过理论方法进行了验证分析.

1 数学模型

文献［2］认为，对于初始速度＜700 m/s的入水问题，可以忽略流体可压缩性的影响，因此，本文数值模拟假设流体不可压缩，忽略入水过程中由于流体黏性产生的热耗散，考虑运动体高速入水引起的空化效应.对于气、汽、液组成的多相流动，本文采用VOF法进行求解.VOF多相流模型将多相流体看作密度可变的单相介质，对流体微团中的各相构成采用体积分数来表示，任一流体微团三相体积分数满足:

式中αg、αv、αl分别为气相、水蒸汽相及水相体积分数.

混合相的连续性方程为

混合相的动量守恒方程为

式中:ui为速度分量;ρm、μm分别为混合相的密度和动力黏度;μt=ρmCμk2/ε为混合相的湍流黏度系数.ρm、μm表达式分别为

根据伯努利方程，在运动体高速穿过水域时，将在周围形成一个低压区，当压力低于饱和蒸汽压时，便会产生空化现象.本文采用Schnerr and Sauer空化模型对流动中的空化问题进行求解，该模型描述水蒸汽相体积分数的输运方程为

式中:RB=1×10-6m为瑞利方程中的气核半径;αnuc=5×10-4为不可凝结气体的体积分数;Fvap=50和Fcond=0.001为经验常数.

对于流动中的湍流现象，本文采用文献［11］提出的k-ωSST湍流模型对流体控制方程进行封闭求解.该模型综合了k-ε模型对充分发展湍流求解的优势和k-ω模型的近壁面稳定性特点，并加入了涡黏度限制方程，能够更恰当的描述湍流剪切应力的传输.

2 数值计算

2.1 计算模型及网格划分

本文针对带有3种不同角度圆锥头型的圆柱体高速入水问题开展了数值模拟研究，运动体几何外形如图1(a)所示，柱段直径D=10 mm，柱段长度L=5D，锥头锥角分别为60°、90°、127°，初始入水速度均为vp|t=0=500 m/s.同时，为排除质量的影响，3种外形的运动体质量取相同数值.计算流场外形采用圆柱体，计算域二维剖面图如图1(b)所示，其中x轴为重力方向，y轴为径向方向.由于入水速度较高，为排除壁面边界对运动体附近流场影响，计算域直径取值Dd=100D.计算域总高度为230D，其中空气域高度为32D，水域高度为198D.初始时刻，保证运动体肩部位置相同，肩部拐点距自由液面高度约为2.25D.

由于运动体及外流场外形均为轴对称体，本文采用二维轴对称模型开展计算，运动体附近网格如图2所示.3种头型条件下采用相同的网格划分方法.运动体附近为空泡形成区域，流动较为复杂，对于这一区域的网格划分，首先在运动体径向6D、轴向2D范围内使用三角形网格加密过渡，运动体壁面首层网格高度为D/400，加密区最外层网格高度为D/20，加密区外为过渡四边形网格，其中沿x轴负向为D/20高度均布，沿径向及x轴正向按一定比例均匀过渡.3 种条件下(θ=60°、90°、127°)网格数量分别为282 502、275 768、273 284.

图1 运动体及计算域示意

图2 运动体附近网格

2.2 数值计算及动网格实现

本文采用有限体积法对离散RANS方程，压力及速度场的求解采用PISO(pressure-implicit with splitting operators)算法［12］，这种算法既具有与迭代的隐式算法相同的精度，又可以取较大的时间步长;对压力场的空间离散采用PRESTO!格式实现;单元内各相体积率离散采用CICSAM格式.计算中嵌入UDF自编程语句，实现运动体沿重力方向的垂直自由运动，在运动过程中，网格的更新与生成采用动态层法实现.

动态层法的主要思想是预设合并或分裂新网格的网格高度临界值，当紧靠运动边界的网格高度低于或超出临界值时，合并或分裂动态层网格来实现网格的更新.其数学表达式为

式中:hj为运动边界的网格高度;hideal为预设的网格更新临界高度;αc、αs分别为网格层合并因子及分裂因子.

网格合并和分裂过程如图3所示，图3(a)、(b)分别为临近运动边界网格层被压缩和拉伸的过程.图3(a)中，当运动边界向上运动，网格高度hj满足式(1)时，第j层网格被压缩并与第i层合并;在图3(b)中，当运动边界向下运动，网格高度hj满足式(2)时，第j层网格就分裂为两层新的网格.

图3 动态层网格更新方法

由于空气域及产生空泡的区域流动较为复杂，为保证该区域计算结果的稳定性，运动体后部网格更新高度均为0.5 mm，即D/20，从而保证该区域网格节点分布保持不变，以提高结果的精确性.

3 结果分析

3.1 入水弹道

图4为3种头型条件下运动体入水速度vp及入水深度的变化曲线，其中，对入水深度H进行了量纲一的处理(H/D).可以看出，角度不同的锥头运动体，入水后速度衰减趋势基本一致，但是速度值有较大差别.锥角越小，阻力系数越小，速度衰减更缓慢，在相同的入水时刻入水深度也更大.随着入水深度的增加，运动体速度逐渐降低，阻力也随之减小，因此在入水初期速度衰减较快，然后逐渐放缓.

图4 不同头型运动体入水弹道曲线

3.2 空泡形态

图5给出了3种头型条件下，不同时刻的入水空泡形态，其中，对空泡半径R进行了量纲一的处理(R/ro).可以看出，由于3种头型阻力系数不同，每种时刻下，运动体的入水深度也有明显的差别，因此在接近运动体端面位置空泡形态产生了交叉.4种时刻下，入水空泡形态均表明，头部锥角越大，上方同一深度处空泡半径越大，而空泡整体轮廓发展规律则保持一致.

图6、7给出了3种头型运动体在4种不同入水深度及4种不同瞬时速度时的入水空泡形态.从图6、7可以看出，当入水深度相同时，头部锥角越大，同一深度处的空泡半径越大.当入水瞬时速度相同时，头部锥角越小，其入水深度越深，同一深度空泡半径也越大;另一方面，自由液面处空泡半径基本相同.

图5 不同时刻下空泡形态对比

图6 不同入水深度下空泡形态对比

为进一步分析3种头型对入水空泡形态的影响，在相关文献基础上对入水空泡半径的影响因素进行分析.

文献［2］从能量守恒的角度给出了入水空泡半径预测公式，其基本假设条件为流体不可压缩，忽略流体黏性引起的能量损耗及重力势能，得出运动体入水过程的动能损失量等于排开水的动能及压力势能之和.式(3)给出了t时刻距自由液面距离为xb深度的空泡半径表达式.

式中:tb为运动体到达xb深度的时刻;r0为运动体特征半径;A(xb)、B(xb)为引入变量，其表达式为

图7 不同瞬时速度下空泡形态对比

将式(4)、(5)代入式(3)，可以得到

忽略重力影响，由牛顿第二定律可得

式中:A0为运动体特征面积;Cdx=Cd0+σ［13］为运动体阻力系数;为空化数;Cd0为空化数为零时的阻力系数.

将式(7)代入式(6)，式(6)右侧可改写为

对于本文研究，在自由液面处，近似取xb=0、tb=0，并考虑到，式(6)可表达为

由于vp(0)=vp|t=0为已知量，因此从式(8)可以看出，当运动体特征面积一定时，在自由液面处，影响空泡半径的主要因素为运动体阻力系数及入水时间.进一步，考虑入水前期空泡扩张阶段，此阶段入水时间较小，可略去pg(0)/(ρlN)及t2两项相对高阶小量，则式(8)可表达为

由式(9)可以得到，对于本文3种头型条件下入水问题，vp(0)可近似认为均等于初始入水速度500 m/s.因此，在同一时刻，3种头型入水过程中，自由液面处空泡半径主要取决于项，即运动体的阻力系数.

根据式(9)，对 3种头型情况下，从入水0.3～1.0 ms这8种时刻自由液面处空泡半径数值模拟结果进行分析，得到了各工况下C0=的值，如图8所示，其中各头型阻力系数Cd0由文献［14］得到，如表1所示.

图8 不同头型各时刻C0

表1 文献［14］中不同锥角在水中阻力系数

从图8可以看出，对于同一锥角运动体，不同时刻下的C0随着入水时间的增大而呈现出逐渐增大的趋势;对于θ=60°及θ=90°的情况，C0在稍微上升后又呈现出下降的趋势，通过对数值模拟结果的分析发现，在C0最大时刻附近是空泡发生了表面闭合，这表明，表面闭合对空泡半径的扩张有一定的影响.对于不同锥角运动体，当入水时间相同时，锥角较大时，C0也较大.

通过以上分析可以看出，对于同一锥角的入水问题，C0的值近似等于某一常数，对于不同的入水时刻，其大小在一定的误差范围内;对于不同的锥角，C0有一定的差异.这表明式(8)简化至式(9)时有一定的误差，对于不同锥角大小，省略部分有一定的差异.

此外，通过常数C0的表达式可知，入水后在相同时刻下，阻力系数越大，自由液面处空泡半径越大.对于本文研究的3种不同头型，锥角较大时，阻力系数也较大，因此图5中，在同一时刻，锥角为127°时空泡半径最大，90°锥角次之，60°锥角时空泡半径最小.

3.3 压力场分布

图9为3种头型运动体入水过程锥头顶点量纲一的压力曲线.从图9(a)可以看出，3种头型条件下，压力变化规律保持一致，在峰值上稍有差别.图9(b)为入水初期0.2 ms以内的局部放大曲线，可以清晰地看出，头部锥角越大，入水压力峰值也越高，θ=127°时可达到约4 000倍大气压力.不同头型条件下，入水冲击压力峰值脉宽基本相同，约为0.02 ms.在该压力峰值之后，3种条件下的压力值趋于一致，但此时，锥角较小时，压力值较高，这是由于在相同的入水时刻下，锥角较小时运动体的瞬时速度较高，头部也相应承受更大的压力.

图9 不同头型运动体入水过程锥头顶点压力变化曲线

图10给出了3种头型运动体入水速度衰减为300 m/s时锥头母线压力系数曲线.压力系数Cp表达式为

式中:p为当地压力;p∞为无穷远处环境压力;vp为运动体速度.

由图10可以看出，3种头型下，运动体锥头母线位置压力系数趋势相同，距离锥头顶端量纲一的距离(l/l0)越小，压力系数越大，在母线末端与柱体连接处，压力系数值在0附近.此外，对比3条曲线可知，在相同的速度下，锥角越大，运动体头部压力系数越大.这表明，锥头角度对运动体入水时受到的压力载荷有较大的影响.

图10 速度为300 m/s时锥头母线压力系数

3.4 速度场分布

图11给出了运动体在水下运动过程中速度衰减为300 m/s时流场速度矢量分布.首先，可以看出，在运动体肩部附近，排开水的速度最大，而空泡壁位置速度较小，且越接近自由液面，速度值越小.这表明，运动体后方空泡扩张受到周围流体静压的作用，扩张速度逐渐减小;其次，对比3种头型肩部位置流体速度值可以看出，不同锥角条件下，排开水的速度值也有较大区别，锥角较大时，排开水的速度也较高.这表明运动体与前方流体撞击分离时损失的能量因为头型的不同而有较大差异，头型阻力系数较大时，能量传递较大，周围水体获得的动能较高，向径向排开的速度也较大.

图11 300 m/s速度时运动体周围及空泡壁速度矢量分布

4 结论

1)锥角不同时，其入水后速度衰减及弹道也有较大差异，锥角越大，其速度衰减也越快，相同时间内运行距离越短.

2)入水空泡半径与头型直接相关，对于不同头型的运动体，入水空泡半径的主要影响因素为阻力系数Cd0，对于同一锥角的头型，在自由液面处，空泡半径r(0，t)近似满足)结果为一常数;对于不同的锥角，角度越大，Cd0越大，得到的C0也越大.

3)入水初期，运动体受到极强的冲击载荷，压力峰值可达数千倍大气压力，且锥角越大，压力峰值越高.

4)相同速度条件下，锥角不同，锥体母线压力系数趋势相同，但在压力系数数值上，锥角越大时，压力系数也越大.

5)运动体在水下运动时，其肩部附近排开水速度远大于后方空泡扩张速度，空泡扩张速度受周围水体影响较大;在相同速度下，锥角较大时，肩部排开水的速度也较大.

［1］ LUNDSTROM E A，FUNG W K.Fluid dynamic analysis of hydraulic ram III.AD 031644［R］.Washington，DC:Joint Technical Coordinating Group for Aircraft Survivability，1976.

［2］ LEE M，LONGORIA R G，WILSON D E.Cavity dynamics in high-speed water entry［J］.Physics of Fluids，1997，9(3):540-550.

［3］DUCLAUX V，CAILLÉ F，DUEZ C，et al.Dynamics of transient cavities［J］.Journal of Fluid Mechanics，2007，591:1-19.

［4］ TRUSCOTT T T.Cavity dynamics of water entry for spheres and ballistic projectiles［D］.Cambridge，MA:Massachusetts Institute of Technology，2009.

［5］LOGVINOVICH G V.Hydrodynamics of free-boundary flows［M］.Jerusalem:Israel Program for scientific translation，1972.

［6］ GUO Z，ZHANG W，WANG C.Experimental and theoretical study on the high-speed horizontal water entry behaviors of cylindrical projectiles［J］.Journal of Hydrodynamics，Ser B，2012，24(2):217-225.

［7］ GUO Zitao，ZHANG Wei，XIAO Xinke，et al.An investigation into horizontal water entry behaviors of projectiles with different nose shapes［J］.International Journal of Impact Engineering，2012，49:43-60.

［8］胡平超，张宇文，袁绪龙.航行器垂直入水空泡特性与流体动力研究［J］.计算机仿真，2011，28(6):5-8.

［9］王 聪，何春涛，权晓波，等.空气压强对垂直入水空泡影响的数值研究［J］.哈尔滨工业大学学报，2012，44(5):14-19.

［10］何春涛，王聪，闵景新，等.回转体匀速垂直入水早期空泡数值模拟研究［J］.工程力学，2012，29(4):237-243.

［11］MENTER F R.Multiscale model for turbulent flows［C］//Proceedings of the 24th Fluid Dynamics Conference，American Institute of Aeronautics and Astronauticsl.Orlando，USA:AIAA，1993:1-21.

［12］ISSA R I.Solution of the implicitly discretised fluid flow equations by operator-splitting［J］.Journal of computational physics，1986，62(1):40-65.

［13］SEDOV L I.Two-Dimensional problems in hydrodynamics and aerodynamics［M］.New York:John Wiley＆ Sons Inc，1965:1-427.

［14］MAY A.Water entry and the cavity-running behavior of missiles，AD A020429［R］.Silver，Spring，MD:Final Technical Report NAVSEA Hydroballistics Advisory Committee，1975.