1 引言
高速入水过程中固、液、气的耦合作用、强大的入水冲击力、液面喷溅的干扰以及自然空化现象等等都使入水过程十分复杂。高速入水在空投鱼雷、跨介质导弹、超空泡射弹等武器中都有广泛的应用，因此近年来引起了大量的关注。
实验研究方面，侯宇等[1]利用高速摄像技术开展了射弹入水试验，分析了弹体不同侧滑角入水冲击过程的喷溅演变、压力波特征等，结果表明侧滑角对弹道轨迹、喷溅形态等都有不同程度的影响。Nguyen等[2-3]利用重叠网格技术，对不同密度的斜楔和圆柱体的入水问题进行了数值模拟并对不同头型弹体水平入水的空化流动做了实验研究，对部分空化和超空化流场做了对比分析。Breton等[4]对类似圆盘状的物体出入水做了实验研究，将实验期间测得的水动力演变及润湿表面半径演变与理论结果进行了比较。Watson等[5]通过试验，研究了半涂层球体的入水问题，结果表明，半涂层球体在入水过程中会经历不对称空腔和水平力，导致入水轨迹的横向偏移。Cherdantsev等[6]通过试验来研究存在空气缓冲情况下的球体低速入水，重点关注了固液接触前自由表面的变形，发现了接触延迟现象。Sui等[7]做了弹丸的高速斜入水实验，分析了载荷与空腔行为间的关系。
数值方法研究高速入水过程也取得了很大的进展。梁景奇等[8]利用动网格移动计算域技术研究了射弹以不同攻角垂直入水的特性，发现了空泡的不对称性。Michael等[9]利用数值模拟分析了空气射流对通气超空泡的影响，结果表明，空气射流会使空泡变得更加细长。Akbari等[10]采用数值模拟的方法，研究了小角度入水圆柱弹的稳定性，结果表明，所研究的弹丸最小稳定入水角为6°。曹雪洁等[11]利用二维轴对称模型分析了水的可压缩性对流场特性、空泡形态、流体动力特性的影响。陈晨等[12]采用二维轴对称模型，同时考虑了气、汽、液三相流体可压缩的情况，对射弹跨声速入水过程进行了数值模拟。李国良等[13]建立了一套考虑液体可压缩性的高速入水数值计算模型，开展了旋成体以50～800 m/s速度入水的数值计算，发现了考虑压缩性会延缓最大载荷的出现时间，同时压缩性对速度、加速度曲线也有影响。
目前对于入水问题的数值研究大多忽略了液体可压缩性的影响。然而，对于高速入水特别是跨音速和超音速入水过程，流体可压缩性的影响是不可忽略的。本文基于计算流体力学软件Fluent，建立超空泡射弹垂直入水的三维六自由度数值计算模型，对比分析是否考虑流体可压缩性时的流场特性，将流体可压缩性对空泡尺寸、流场压力、流场密度等物理参数的影响进行定量分析。
2 数值模型与理论方法
2.1 数值模型
本文研究的内容为跨介质流动，流场中存在空气、液态水、水蒸气共三相。应用VOF多相流模型，并对RANS(reynolds-averaged navier-stokes)方程进行求解，相间界面是通过求解各次相的连续性方程来确定的，对于第q相有：
(1)
式中：αq为第q相的的体积分数； ρq为第q相的密度；
为第q相速度；
为单位时间内从第j相传递到第q相的质量，
同理。通过式(1)，即可求出控制体内各次相的体积分数，而主相的体积分数则通过下式计算：
(2)
动量方程为：
(3)
式中： ρ为密度； p为压力； μ为动力粘度；
为速度；
为重力加速度；
为其余外力。
由于考虑了流体的可压缩性，故需求解能量方程，能量方程的形式如下：
(4)
式中：E为总能量；ke为有效热导率；hw,q为q相中物质w的焓；
为q相中物质w的扩散通量；τe为粘性剪切率；等号右边前三项分别表示由热传导、物质扩散、粘性耗散引起的能量传递；S为源项，包括辐射以及其他体积热源。
由于雷诺时均法是将瞬时速度处理为平均速度与脉动速度的叠加，故RANS方程组中多出了雷诺应力这一未知项，需要使用湍流模型对方程组进行封闭，Realizable k-ε模型适用于快速应变情况，可较好的求解复杂流动，且计算开销较小，因此选取Realizable k-ε模型用于后文的计算。
当水的压力减少到饱和蒸汽压以下时，液态水便会发生汽化形成水蒸气，这种现象即称为空化现象，高速入水过程中包含着空化现象的发生。本文选取Schnerr-Sauer空化模型来求解此物理现象：
(5)
式中：R为空化源项；Fvap和Fcond分别为蒸发和冷凝的经验校准系数，默认取值分别为1和0.2；ρv为水蒸气密度，ρl为液态水密度；α为水蒸气相的体积分数；Qb为气泡半径；ph为当前环境压力；pv为水的饱和蒸汽压。
由于需要比较可压缩性与不可压缩性对流场特性的影响，当考虑各相的可压缩性时，需要引入状态方程。入水前，射弹只受到空气的作用力，顾需要考虑空气的可压缩性，本文将空气视作理想气体，有理想气体状态方程：
(6)
式中：R为气体常数；T为温度。
入水之后，由于液态水对弹的作用力远远大于空气及水蒸气，此时可忽略气体的可压缩性，仅考虑液态水的可压缩性。对于液态水，则引入Tait液体状态方程来描述液体密度和压力的关系，Tait方程的原型是P.G.Tait[14]于1888年通过大量的实验数据及经验推导出来的公式，其简易形式如下：
(7)
式中： ρ为当前密度；ρ0为参考密度；n为密度指数，对于液态水，取n=7.15；K0为参考体积模量；K为当前体积模量，有:
K=K0+n(p-p0)
(8)
式中：p为绝对压力；p0为参考压力。
激波的产生与扰动的传播速度有关，扰动在水中的传播速度为水中声速，按下式计算：
(9)
2.2 理论方法
关于空泡的形态，前人已经总结出了多种不同形式和适用范围的经验公式，本文基于细长体理论[15]对空泡形态进行预测，用于后续的数值方法验证，以增加数值模拟的可信度。
细长体理论是针对圆盘空化器的考虑可压缩性的预测空泡形态的方法，其计算公式为：
(10)
式中：
为空泡半径；Cd为阻力系数，
其中
Rn为空化器半径；
为空化数，其计算方法为：
(11)
式中：p∞为无穷远处来流压力。
3 模型介绍与方法验证
为尽可能接近真实的物理场和射弹运动过程，采用三维计算模型，利用6自由度运动方程求解射弹的入水运动。
3.1 模型介绍
射弹模型如图1所示，射弹全长L=200 mm，前锥段长65 mm，圆柱段直径D=12 mm，头部空化器直径=4 mm。
图1 射弹模型示意图
Fig.1 Projectile model
据黄闯等[16]研究，过小的流场计算域尺寸会导致壁面阻塞效应，对超空泡形态具有一定的影响。为减轻壁面阻塞效应所带来的影响，本文构建的外流场形状为长宽相等的长方体，流场的长、宽尺寸均为75倍最大弹径；流场高度为10倍弹长，其中，空气区域高度为2.4倍弹长，水相区域高度为7.6倍弹长，具体的流场计算域如图2所示。初始时刻(t=0时)，射弹距水面高度为0.75倍弹长；流场顶面采用压力入口，侧面与底面均采用固壁边界。
图2 计算域示意图
Fig.2 Computational domain
利用结构化重叠网格技术对上述计算域进行网格划分，射弹近壁面处、水面区域及运动路径的网格进行加密处理，最终网格数量在400万左右，网格划分情况如图3所示。采用与地面固连的绝对坐标系，沿初始时刻弹轴向上方向为y轴正方向，x-z面与水平面平行，坐标原点位于初始时刻射弹头部空化器的中心。
图3 网格示意图
Fig.3 Schematic diagram of grids
3.2 数值方法验证
利用张伟等[17]关于弹体高速入水实验研究所获得的实验数据与数值计算结果进行对比，以验证数值方法在较低马赫数(0.2马赫)工况下的可靠性。试验所用弹型为圆柱弹，直径12.65 mm，长25.4 mm，以初速118.8 m/s水平入水，实验获得了空泡轮廓尺寸以及速度随时间变化的数据。图4为数值模拟的速度衰减曲线与实验数据曲线，其中v表示射弹运动的速度，可以看到，数值模拟的结果与实验数据吻合较好；图5为数值模拟的空泡轮廓尺寸与实验数据曲线，由图5可以看到，两者的形态及尺寸基本一致。
图4 实验与模拟速度衰减曲线
Fig.4 Comparison of velocities between experiment and simulation
图5 实验与模拟空泡轮廓曲线
Fig.5 Comparison of experimental and simulated cavitation profiles
关于弹体跨音速及超音速入水的试验研究较为缺乏，难以获得相关数据，因此利用1.2节所提到的理论公式计算高马赫数工况(1马赫、1.3马赫)，验证数值计算结果。利用本文构建的考虑可压缩性的数值模拟方法对初始速度为 1 480 m/s的入水工况进行计算，所得的空泡形态尺寸与细长体理论计算结果如图6所示，由图中曲线可以看出空泡的形状轮廓十分吻合。
图6 理论与模拟空泡轮廓曲线
Fig.6 Comparison of theoretical and simulated cavitation profiles
综上，本文所构建的用于求解入水问题的数值模拟方法较为可信，可用于后续的计算研究。
4 结果与分析
基于上述数值计算模型，分别对初始时刻速度为300 m/s、1 480 m/s、1 921 m/s的射弹入水问题进行求解，3种工况分别对应水中的亚音速(0.2马赫)、临界音速(1马赫)以及超音速(1.3马赫)。分别采用不可压缩流体介质与可压缩流体介质模型对上述3种工况进行求解。下文凡是取某一特定时刻的流场作为研究对象的，均为射弹在水中运动距离为2L时的流场。通过对结果的分析比对，研究可压缩性对入水过程流体动力的影响。
4.1 可压缩性对空泡形态的影响
分别比较各马赫数下是否考虑可压缩性所得到的空泡轮廓如图7所示。
图7 可压缩、不可压缩3D空泡轮廓示意图
Fig.7 Comparison of compressible and incompressible cavitation profiles
从图7中可以看出，速度较低(0.2马赫)时，是否考虑可压缩性所得到的空泡轮廓尺寸相差不大，可压缩工况中弹尾处的空泡直径比不可压缩的仅小了约0.8 mm；当速度达到1马赫时，考虑可压缩性的空泡尺寸明显小于不考虑可压缩性时所得到的空泡轮廓尺寸，在弹尾处两者直径差值大约达到了2.35 mm；随着马赫数的继续增加，到1.3马赫时，这种压缩效应表现的更为明显，两者差值已达到了7 mm。这是因为1马赫、1.3马赫工况中考虑可压缩性的流场中出现了激波，如3.2节的图8所示，这会使高压区发生后掠，进而使空泡侧面液体流场压力增加，最终使1马赫、1.3马赫中可压缩工况的空泡径向尺寸明显变小。
4.2 可压缩性对流场特性的影响
图8为射弹头部附近流场的压力分布云图，每组图左侧为弹头附近压力云图，云图中白色部分代表弹体，右侧为沿射弹轴线(云图中的红色竖线)的流场压力曲线。从图8(a)中可以看出，当射弹以0.2马赫数入水时，是否考虑可压缩性对流场的压力分布影响不大，头部最大压力仅相差0.1 MPa。从图8(b)和图8(c)中可以看出，当射弹以1马赫及1.3马赫数入水时，是否考虑可压缩性计算得到的流场压力分布差异较大。不考虑可压缩性时，压力虽然在接近弹头处有较大的上升，但仍然是连续的，说明此时的流场只是受到弹体运动所产生的压缩波的影响，并没有形成激波；当考虑可压缩性时，弹头前方压力发生了极大的突变，说明此时流场中出现了脱体激波，这可以从压力云图中更直观的看出。
图8 弹头附近流场压力云图及其流场压力曲线
Fig.8 Pressure contour and its profile near the warhead
更进一步的，利用曲线图中弹头及激波面的位置坐标得到，1马赫时激波面位于沿弹轴方向弹头前方约34.9 mm处；1.3马赫时激波面位于弹头前方约9.6 mm处。说明随着马赫数的增加，激波有向弹体靠近的趋势。
对于可压缩性而言，一项比较重要的流场参数是密度，不考虑可压缩性的流场密度是恒定值，考虑压缩性时，流场的最大密度出现在射弹的头部中心紧前方位置处。如图9所示，随着马赫数的增加，流场最大密度也逐渐增加，从图还可看出：对于本文所研究的弹型，入水后流场的最大密度与射弹马赫数近似成正比关系，还可得到近似预测流场最大密度的数学关系式：
ρm=288.8Ma+956.86
(12)
式中： ρm为射弹入水后流场最大密度； Ma为射弹当前马赫数。
为验证该公式的准确性，增加了0.6及1.5马赫数的工况以额外获取2个最大流场密度与马赫数之间关系的数据，现将预测公式与数值模拟的结果画于图10。由图10可见，数值模拟的结果与预测公式较为吻合，最大误差不超过0.56%，说明本文提出的流场最大密度预测公式，对所讨论的弹型，在0.2～1.5马赫数范围内是适用的。
4.3 可压缩性对射弹动力特性的影响
图11为不同马赫数下是否考虑流体压缩性的阻力系数曲线，其中阻力系数定义为
为射弹所受总阻力，v为射弹运动速度，S取射弹头部面积为参考面积。如图11所示，亚音速入水时是否考虑压缩性的阻力系数相差不大；而达到音速及超音速时考虑可压缩性的阻力系数明显大于不考虑可压缩性时的阻力系数，且随着马赫数的提高，这种差异愈加明显。这是因为不可压缩计算模型认为波的传递速度趋近于无穷大，射弹即使达到音速后也远远不能追赶上压力波，因此压力波无法在射弹前方叠加形成激波，故在达到音速以后，其阻力系数明显小于可压缩模型的计算结果；在本文的计算范围内马赫数越高则激波强度越高，因此随着马赫数的增加，是否考虑压缩性的阻力系数差异愈加明显。
图12为入水过程射弹头部最大压力曲线。可以看出，在0.2马赫时，压缩性对最大压力几乎没有影响。速度超过1马赫后，考虑可压缩性的头部最大压力更大，且随着射弹速度的增加，两者间的差距也逐渐加大，在1.3马赫时，两者相差19.5%。这是因为在不可压缩模型计算中压力波的传播速度远大于可压缩模型；在可压缩模型中，亚音速入水时压力波不会因为追赶叠加形成激波，故此时是否考虑压缩性对头部最大压力影响不大；而在音速或超音速入水过程中，射弹速度等于或大于水中压力波传播速度，压力波叠加在弹头形成激波，从而可压缩模型计算得到的射弹头部压力更大。
图9 不同马赫数下流场最大密度散点示意图
Fig.9 Scatter diagram of maximum density of the flow field at different Mach numbers
图10 流场最大密度预测公式与数值模拟结果示意图
Fig.10 Comparison between prediction formula and simulation results
图11 不同马赫数可压/不可压阻力系数曲线
Fig.11 Comparison of compressible and incompressible drag coefficients at different Mach numbers
图12 不同马赫数下弹头部最大压力曲线
Fig.12 Maximum pressure of warhead at different Mach numbers
射弹头部的最大压力对弹体的结构设计、材料选用等具有重要的参考价值，因此针对本文所研究的弹型及所讨论的速度范围，对考虑流体可压缩性的、射弹水中航行时头部最大压力与射弹航行马赫数的关系进行推导如下，将式(7)进行简单的变形可得：
(13)
取ρ0=998.2；p0=101 325；K0=2.2×109，代入式(13)得：
(14)
联立式(12)与式(14)即可求得入水后射弹头部最大压力与射弹马赫数的数学关系式：
(15)
式(15)即为考虑流体可压缩性情况下的射弹水中航行时头部最大压力与马赫数之间的数学关系式，现利用本文数值计算的结果对该公式的准确性进行验证。
如图13所示，式(15)的预测结果与数值模拟的结果较为吻合，在本文所验证的范围内最大误差不超过5.04%，利用此公式可减少数值模拟所消耗的计算成本，对该弹型的后续研究有指导意义。
图13 射弹头部最大压力预测公式与数值模拟结果曲线
Fig.13 Comparison of prediction formula and simulation results
5 结论
1) 同马赫数下，考虑流体可压缩性后的空泡轮廓尺寸减小；射弹阻力系数增大。
2) 提出了预测入水后流场最大密度的数学关系式，并推导了在考虑流体可压缩性的情况下，水中航行时射弹头部最大压力与射弹马赫数之间的关系式，可为本文弹型的后续研究提供参考。
3) 考虑可压缩性时流场中出现激波，且随着马赫数的增加，激波逐渐向弹体靠近；马赫数越高，可压缩性对射弹入水流场的影响越明显。
[1] 侯宇，黄振贵，郭则庆，等.超空泡射弹小入水角高速斜入水试验研究[J].兵工学报,2020,41(02):332 -341.
HOU Y,HUANG Z G,GUO Z Q,et al.Experim-ental Investigation on shallow-angle oblique water-entry of a high-speed supercavitating projectile[J].
Acta Armamentarii,2020,41(02):332-341.
[2] NGUYEN V,PHAN T,PARK W.Modeling and numerical simulation of ricochet and penetration of water entry bodies using an efficient free surface model[J].International Journal of Mechanical Sciences,2020,182.
[3] NGUYEN T T,HAI D N,THAI N Q,et al.Experimental measurements of the cavitating flow after horizontal water entry[J].Fluid Dynamics Research,2017,49(05).
[4] BRETON T,TASSIN A,JACQUES N.Experimental investigation of the water entry and/or exit of axisymmetric bodies[J].Journal of Fluid Mechanics,2020,901.
[5] WATSON D A,BOM J M,WEINBERG M P,et al.Water entry dynamics of spheres with heterogen-eous wetting properties[J].Physical Review Fluids,2021,6(04):44003.
[6] CHERDANTSEV A V,GAVRILOV N V,ERMAN-YUK E V.Study of initial stage of entry of a solid sphere into shallow liquid with Synthetic Schlieren technique[J].Experimental Thermal and Fluid Science,2021,125:110375.
[7] SUI Y T,ZHANG A M,MING F R,et al.Experi-mental investigation of oblique water entry of highs-peed truncated cone projectiles:Cavity dynamics and impact load[J].Journal of Fluids and Structures,2021,104.
[8] 梁景奇，王瑞，徐保成，等.攻角对高速射弹入水过程影响研究[J].兵器装备工程学报，2020，41(07):23-28.
LIANG J Q,WANG R,XU B C,et al.Research on influence of angle of attack on process of high-speed water-entry projectile[J].Journal of Ordnance Equipment Engineering,2020,41(07):23-28.
[9] MICHEAL P K,MICHEAL H K,Ivan N K,et al.A numerical assessment of the interaction of a super-cavitating flow with a gas jet[J].Ocean Engineering,2017,136.
[10] AKBARI M A,MOHAMMADI J,FEREIDOONI J.Stability of oblique water entry of cylindrical projectiles[J].Journal of Applied Fluid Mechanics,2021,14(1):301-314.
[11] 曹雪洁，胡俊，于勇.液体可压缩性对高速射弹超空泡流的影响[J].兵工学报，2020,41(S1):72-78.
CAO X J,HU J,YU Y.Effect of liquid compres-sibility on supercavitation flow of high-speed projectiles[J].Acta Armamentarii,2020,41 (S1):72-78.
[12] 陈晨，魏英杰，王聪，等.射弹跨声速入水初期阶段多相流场特性数值研究[J].振动与冲击，2019,38(06):46-53.
CHEN C WEI Y J,WANG C,et al.Numerical analysis of the multiphase flow characteristicsin the initial period of water-entry of a projectile at transonic speed[J].Vibration and shock,2019,38 (06):46-53.
[13] 李国良，尤天庆，孔德才，等.旋成体高速入水可压缩性影响研究[J].兵工学报.2020,41(04):720-729.
LI G L,YOU T Q,KONG D C,et al.Effect of fluid compressibility on high-speed water-entry of revolutionary body[J].Acta Armamentarii,2020,41(04):720-729.
[14] 何星.基于实验的Tait、Schmidt和Renner模型比较研究[D].杭州：浙江大学，2017.
HE X.Comparison of Tait,Schmidt and Renner approaches based on experiment[D].Hangzhou:Zhejiang University,2017.
[15] SEREBRYAKOV V,KIRSCHNER I,SCHNERR G.High speed motion in water with supercavitation for sub-,trans-,supersonic Mach Numbers[J].2009.
[16] 黄闯，罗凯，党建军，等.流域径向尺度对自然超空泡的影响规律[J].西北工业大学学报，2015,33(06):936-941.
HUANG C,LUO K,DANG J J,et al.Influence law of watershed radial scale on natural supercavitation[J].Journal of Northwestern Polytechnical University,2015,33(06):936-941.
[17] 张伟，郭子涛，肖新科，等.弹体高速入水特性实验研究[J].爆炸与冲击，2011,31(06):579-584.
ZHANG W,GUO Z T,XIAO X K,et al.Experi-mental study on the characteristics of projectile entering water at high speed[J].Explosion and Shock-Waves,2011,31(06):579-584.