高速旋转射弹波浪入水多相流场与弹道特征数值仿真研究

李治涛， 赵世平， 卢丙举， 于 勇

(1. 中国船舶集团有限公司第七一三研究所，郑州 450015；2. 河南省水下智能装备重点实验室，郑州 450015； 3. 北京理工大学 宇航学院，北京 100081)

射弹高速入水会产生空化现象，在真实的海况条件下，还会受到波浪等因素的影响。射弹在波浪中入水，除受浮力、自身重力、砰击载荷外，还要考虑波浪载荷对射弹的影响。波浪载荷对射弹入水过程结构的安全性以及弹体入水弹道的稳定性均有非常重要的影响。

储慧林[1]基于流体体积函数(volume of fluid, VOF)法建立了二维数值水槽，研究了波浪相位、波高对鱼雷入水过程砰击载荷、弹道和速度的影响。王平等[2]基于VOF法，建立了楔形体自由入水模型，研究了入水时楔形体所受的载荷、自由液面的变化以及物面压强分布等。杨衡等[3]采用非线性双渐进法，研究了三维圆柱体在波浪、洋流下入水运动速度及轨迹。王文华等[4]采用二维计算流体动力学(computational fluid dynamics,CFD)方法，数值模拟了圆柱在规则波中的入水过程，研究了波浪参数(波高、周期和入水相位)对圆柱入水动力学特性的影响。金禹彤等[5]用CFD方法，数值模拟了楔形体波浪入水过程，研究了规则波波形位置对楔形体受力特性和运动姿态变化的影响。邹丽等[6]通过具有造波功能的水槽，开展了两种曲面楔形体分别在静水和波浪中入水砰击试验，记录了入水过程楔形体的加速度和压力变化，用高速摄像系统捕捉砰击过程的射流飞溅现象。侯昭[7]采用大涡模拟方法，数值模拟了圆柱低速波浪入水，研究了入水空泡演变机理以及漩涡特征等。郭张霞等[8]用数值方法对回转体波浪入水进行了研究，波浪对回转体的空泡特征、头部压力以及诱导气相速度产生影响。朱仁庆等[9]用数值方法对不同刚度三维楔形体波浪入水进行了研究，弹性效应在一定程度上会减缓砰击应力。刘邦威[10]基于势流理论和边界元法开发了物体三自由入水过程的全非线性数值模拟程序，对三维圆锥波浪入水进行了研究，分析了波高对物体入水过程的影响。史崇镔等[11]用CFD法对跨介质航行器波浪入水过程进行了数值模拟，分析了有无波浪和波高对航行器入水运动特征的影响。赵庆新等[12]用CFD法对水下运载器波浪入水过程进行了数值模拟，研究了波浪参数和风速等对水下运载器受力的影响。陈晨等[13]用能够考虑气、汽、液三相可压缩性数值方法对射弹跨声速入水初期多相流场进行了数值模拟，并分析了入水高度的影响。

对波浪入水问题，目前研究较多的是楔形体、圆柱、航行器、圆锥等，高速旋转射弹波浪入水研究较少。本文用定义输入边界造波法、重叠网格技术和6自由度算法对某一高速旋转射弹波浪入水过程进行数值模拟，研究了波峰、波节1、波谷、波节2等波面不同位置入水弹体空化形态、弹道特性和流体动力载荷特性的变化规律，研究结果对复杂海况下射弹入水问题具有一定的参考价值。

1 数值计算方法

本文采用均质多相流理论中的流体体积函数模型对空气、水蒸气、水组成的多相流动系统进行描述，流动被看作单一介质的流动系统，不同流动组分共享一套动量方程，通过计算得到单元内各相流体介质的体积分数，从而确定不同相在流动系统中的分布。

1.1 控制方程组

描述混合介质流动的N-S方程组为

(1)

式中：v为速度分量；ρm=αlρl+αvρv+αgρg和μm=αlμl+αvμv+αgμg分别为混合介质密度和动力黏度；ρl，ρv和ρg分别为水、水蒸气和空气的密度；μl，μv和μg分别为水、水蒸气和空气的动力黏性系数；μt为湍流黏性系数；αl，αv和αg分别为水、水蒸气和空气相的体积分数。对于任意控制体均有如下关系成立

αl+αv+αg=1

(2)

体积分数的变化满足如下输运方程

(3)

数值计算中采用剪切应力输运k-ω湍流模型(shear-stress-transportk-ω, SSTk-ω)。该模型考虑了湍流剪切应力效应，对很多湍流流动问题的求解具有很高的精度和可信度。其数学方程为

(4)

式中：k为湍动能；ω为湍动能耗散率；υt=μt/ρ为湍流动力黏性系数；σk，γ，σω和σω2为湍流模型常数。

物体在水下高速运动时会引起空化现象，需要引入空化模型，用于描述水相与水蒸气相之间的质量输运。本文选用Schnerr-Sauer空化模型，其表达式为

(5)

式中：RB=1×10-6为空化气核半径；αnuc=5×10-4为不可凝结气体的体积分数；Fvap=50和Fcond=0.001为经验常数；饱和蒸汽压pv取3 540 Pa。式(5)等式右侧第一项是水蒸气的生成项，第二项是水蒸气的凝结项。

1.2 重叠网格技术

弹体入水是一个极强的非定常、非线性问题。弹体入水运动过程中弹体与水域接触面的边界位置，需要根据之前弹体边界所受的载荷计算获得。弹体每个时刻在水域中的位置均不同，用重叠网格技术可以模拟弹体实时运动引起的流域形状随时间变化。

重叠网格法将弹体入水过程复杂的流动区域分成流场区域和弹体周围区域等两个简单的子区域。各子区域中的网格单独生成，彼此存在着重叠或嵌套关系，流场信息通过插值在重叠区边界进行匹配和耦合，各子区域无须进行繁琐的拓扑分区，降低了网格生成难度[14-15]。本文在后续计算中流场区域即背景网格，弹体周围区域即部件网格，二者均采用结构网格。

1.3 6自由度算法

计算过程中部件网格区域随弹体同步运动，部件和背景网格之间的挖补插值操作会在每一个时间步进行，弹体运动采用6自由度算法计算。

采用地面坐标系和弹体坐标系来描述弹体的运动规律。地面坐标系OXYZ与CFD计算采用的坐标系一致如图1所示。其中：OX轴在水平面内；OY轴在纵平面内，方向垂直向上，重力为-Y方向；OZ轴满足右手坐标系。弹体坐标系原点位于弹体质心，坐标系OX1轴与弹体轴线相一致，指向弹体尾端；OY1轴垂直于弹体横向平面，方向垂直向上；OZ1满足右手坐标系。

图1 地面坐标系和弹体坐标系Fig.1 Ground coordinate system and body coordinate system

射弹质心处的速度为V，根据动量定理在地面坐标系中有

(6)

式中，F为作用于射弹上的合外力。依据惯性系(地面坐标系)与动坐标系(弹体坐标系)中速度导数关系有

(7)

将式(7)代入式(6)，得

(8)

式中，Fx，Fy和Fz为作用于弹体的外力矢量在弹体坐标系中的分量。

惯性系下弹体旋转运动的动量矩方程为

(9)

式中：M为合外力对弹体质心的力矩；J为射弹的转动惯量。

根据惯性系(地面坐标系)中与动坐标系(弹体坐标系)中动量矩导数的关系有

(10)

将式(10)代入式(9)，得

(11)

式中，Mx，My和Mz为外力对射弹质心力矩M在弹体坐标系中的分量。

2 数值方法验证

用Schnerr-Sauer空化模型、SSTk-ω湍流模型和重叠网格技术，建立了高速射弹跨介质入水空化多相流场数值仿真模型。为了验证该空化多相流场数值仿真模型的可靠性和准确性，首先对其进行了验证。用文献[16]中的平头圆柱高速入水问题作为验证算例。

2.1 平头圆柱模型

采用郭子涛研究中所用的平头圆柱，其参数如表1所示。

表1 平头圆柱参数Tab.1 Parameters of flat-headed cylinder

2.2 计算网格及边界条件

部件网格为包含平头圆柱的圆柱形区域，平头圆柱位于该计算域的中心，计算域直径为15倍平头圆柱长度即0.381 m，高度为11倍平头圆柱长度即0.279 4 m，网格数目为310 374，如图2(a)所示；背景网格区域为长方体，长为0.76 m，宽为0.76 m，高为1.20 m，网格数目为1 110 262，其中水域高为1.00 m，空气域高为0.20 m；背景网格与部件网格初始化后在Z=0截面分布，如图2(b)所示，Z轴满足右手坐标系。

边界条件设置如图3所示。计算域X轴最大值与最小值的边界面、Z轴最大值与最小值的边界面均设为对称边界条件；Y轴最大值边界面设置为压力出口，压强大小与标准大气压一致；Y轴最小值边界面设置为压力出口，压强为对应水深的静水压力；平头圆柱边界设置为Wall；部件网格外边界设置为Overset边界；计算时间步长为1×10-6s。

图2 平头圆柱计算域网格分布Fig.2 Grid distribution in computational domain of flat-headed cylinder

图3 边界条件设置Fig.3 Boundary condition setting

2.3 入水速度与位移验证

图4分别给出了平头圆柱体跨介质入水过程中，弹体入水位移曲线和速度衰减曲线。从图4中可以看出数值计算结果与试验具有较好的一致性。验证了用Schnerr-Sauer空化模型、SSTk-ω湍流模型和重叠网格技术所建入水空化多相流场数值仿真模型的可靠性和准确性。

图4 平头圆柱入水的位移与速度变化曲线Fig.4 Displacement and velocity curves of flat-headed cylinder entry water

3 数值造波

定义输入边界造波是根据行波的解析解，给定造波边界处的值，从而实现造波的一种方法[17-22]。它适用于规则波、非规则波等多种波浪模型[23]，并且还能够适应波浪和洋流耦合，在波浪和海洋工程领域应用广泛。

3.1 入射边界条件

采用定义输入边界造波法考虑到第4章射弹入水的三维数值模拟，造波采用三维计算，忽略Z向(侧向)的波浪运动，即只考虑二维波浪情形。

对于二阶Stokes波浪，入射边界处的速度和波高满足式(12)。其X和Y方向速度与波高分别为

(12)

式中,T，k，H和d分别为波浪周期、波数、波高和静水深。

定义输入边界造波法的边界条件设置，如图5所示。用Fluent中的宏DEFINE_PROFILE(name,t,i)定义入射边界处流体的速度；用宏DEFINE_PROFILE(name,t,i)定义入射边界处流体的体积分数，达到控制入射边界处流体高度，液相状态下水的体积分数为1，空气的体积分数为0。

图5 定义输入边界造波法示意图Fig.5 Schematic diagram of wave-generation method of defining inlet boundary conditions

3.2 海绵层阻尼消波

波浪传播至水槽末端，如不进行消波处理，会造成波浪的反射。反射波浪的存在会与入射波浪形成波浪叠加，影响数值水槽的波形和流场。海绵层阻尼消波是在特定的计算区域内添加阻尼项以达到消弱或消除该区域波动。本文将计算域边界前1～2倍波长的范围设置为消波区域，在动量方程中添加阻尼项，通过Fluent中的宏DEFINE_SOURCE(mom_source,c,t, dS, equ)实现消波。动量方程中添加人工黏性为

(13)

式中：α(x)为消波系数，采用线性分布的形式；x1和x2分别为海绵层消波的起始坐标和末端坐标；β为经验常数，与波浪特征等因素相关，一般取10。

3.3 造波计算及结果

数值模拟5级浪二阶Stokes波，该波的波长、波高和周期分别为40 m，3 m和5.06 s。造波区域长160 m、宽5 m、高20 m，其中水深15 m、空气域高5 m，坐标原点位于静液面对称中心处。计算采用的网格分别在静液面上下半波高附近进行了加密，网格数目为2 381 727，如图6所示；造波采用的边界条件与图5相同，入口边界条件的速度剖面、压强分布、体积分数等参数均通过自定义函数(user-defined function, UDF)载入。计算采用隐式VOF模型、SSTk-ω湍流模型、含分离算子的隐式压力(pressure-implicit with splitting of operators, PISO)算法，时间步长为0.01 s；当有较为稳定的波形出现时，逐步降低时间步长一直到1×10-6s，以适应第4章射弹高速入水的时间步长。采用水相体积分数为0.5的位置作为波面位置，二阶Stokes波面历时演化过程如图7所示。

图6 数值水槽网格Fig.6 Numerical flume grid

图7 二阶Stokes波数值生成历时演化过程Fig.7 Evolution of second-order Stokes wave generation

仿真所得T=50 s波形与理论波形对比，如图8所示。图中X坐标从20～60 m数值计算时设置为消波区即图中灰色区域，受消波与反射波浪的影响X坐标从20～80 m内仿真所得波形与理论波形差别较大；图中X坐标从-80～0 m内数值计算所得波形与理论波形吻合较好。计算所得波形能较好地反应波浪波面变化，可用于第4章高速射弹波浪入水仿真分析。

图8 数值波形和理论波面对比Fig.8 Comparisons of numerical and theoretical wavefront

4 高速旋转射弹波浪入水仿真分析

选取高速射弹在波浪波面四个不同位置入水进行数值仿真研究。波面入水位置如图9所示，分别为波峰、波节1、波谷、波节2。

计算所用射弹模型如图10所示，弹体长为Lt，头部直径为Dn，弹体最大直径为Dt，弹体细长比为8.746，弹体质心到头部距离与弹体长度的比值为48.37%。弹体初速度为408 m/s，名义入水角为25°(弹体轴线与水平面的夹锐角为名义入水角)，初始旋转速度为7 500 r/min。

图9 高速旋转射弹波面不同入水位置Fig.9 Different water-entry positions of wave surface of high-speed spinning projectile

图10 弹体几何示意图Fig.10 Geometric sketch of projectile

4.1 射弹波浪入水弹体网格与计算方法

高速旋转射弹波浪入水，采用1.2节所述重叠网格技术，网格分为背景网格和部件网格。背景网格见图6；部件网格区域是包围弹体周围一定区域，使弹体周围有较好质量且足够精细的网格用于计算。采用三套不同大小流场区域的部件网格，它们的几何尺寸如图11所示，图中D为弹体最大直径。虽然三套部件网格流场区域大小不同，但它们的网格拓扑是一样的。这里仅给出第3套网格分布如图12所示，第1套和第2套网格与其类似不再附图。图12分别给出了过弹体重心轴截面的网格分布和其他两个剖面的典型网格，这样划分的结构化网格质量很高。

图11 部件网格流场流域示意图Fig.11 Schematic diagram of flow field area of component grid

图12 部件网格局部截面分布示意图Fig.12 Schematic diagram of component grid local section distribution

3套不同大小流场区域的部件网格信息，如表2所示。背景网格尺寸大小的选取应与部件网格最大网格尺寸保持基本一致。这样可以保证在背景网格和部件网格交界面上插值的精度，在提高部件(弹体)周围流场分辨率的同时降低总的网格数目。

表2 部件网格信息Tab.2 Component grid information

图13给出了3套不同部件网格，射弹以初速度为408 m/s、名义入水角7.2°和初始旋转速度为7 500 r/min平静水面工况入水，弹体质心运动轨迹在XOY面的投影。从图13中可以看出，受部件网格区域大小的影响第1套网格弹体在XOY平面比第2套和第3套网格更早缓慢抬头，第2套和第3套网格弹体运动轨迹几乎没有差别。增加部件网格区域的尺寸势必增大网格的数量，加大计算负担，所以部件网格尺寸的选取需要综合平衡计算精度和计算资源的取舍。事实上，当部件网格尺寸达到一定值后，再增加部件网格尺寸对计算结果基本不会有影响。经过上述分析，又考虑到第3套网格仅比第2套增加了5万多网格量，计算量增加不是太大，后续仿真计算中均采用第3套部件网格。

图13 3套不同部件网格，弹体质心运动轨迹在XOY面的投影Fig.13 Three sets of different component grid, projection of projectile centroid trajectory on XOY plane

将背景网格和部件网格同时导入Fluent软件中，设置弹体初始位置和姿态，弹体不动，用第3章所述定义输入边界造波法进行数值造波，在造出稳定波形后逐步减小时间步长，使其与射弹波浪入水时间步一致为1×10-6s；之后设置弹体6自由度运动，打开空化模型、三相流体等，进行射弹波浪入水仿真，直到计算结束。高速旋转射弹入水位置处的波形(波峰、波节1、波谷和波节2等)可用式(12)中的波高预示获得。

4.2 空化效应对高速射弹入水特征的影响

图14给出了是否考虑空化效应，高速射弹波峰入水Z=0平面不同相的体积分数云图。图14(a)未考虑空化效应，即计算中未开启1.1节Schnerr-Sauer空化模型，该图中深灰色为水相、白色为空气相；图14(b)考虑空化效应，即计算中开启1.1节Schnerr-Sauer空化模型，该图中浅灰色为水相、深灰色为空气相、白色为水蒸气相。通过对比图14中是否考虑空化效应可发现：未考虑空化效应射弹波峰高速入水，刚入水时弹尾携带部分空气图中白色，随着入水深度增加弹尾空气逐渐从弹尾断裂溢出，入水9 ms以后弹尾几乎没有空气相，入水后受水动力的作用弹体从6 ms后逐渐缓慢抬头(见图14(a))；考虑空化效应，射弹刚入水时弹体头部产生了空化效应，即水相(图中浅灰色)汽化为水蒸气(图中白色)，同时弹尾携带部分空气(图中深灰色)，由于空化效应在12 ms内弹体仍保持在超空泡中直线航行。

图14 是否考虑空化效应，高速射弹波峰入水Z=0平面不同相体积分数云图Fig.14 Whether cavitation effect is considered, volume fraction contour of different phases in the Z=0 plane of wave crest entering water of high-speed projectile

图15给出了是否考虑空化效应，弹体波峰入水其质心运动轨迹。从图15中可以看出，未考虑空化效应由于弹体高速旋转其受到较大的侧向水动力载荷作用，使其产生了侧向位移，12 ms内弹体质心产生了约40 mm的侧向位移；考虑空化效应，12 ms内弹体质心几乎未产生侧向位移。

图16给出了是否考虑空化效应，弹体波峰入水其质心运动轨迹在XOY平面的投影。从图16中可以看出，未考虑空化效应，弹体入水一段时间后由于水动力载荷的作用其缓慢抬头；考虑空化效应，弹体在12 ms内沿弹体初始对称轴直线航行。

图15 是否考虑空化效应，弹体波峰入水其质心运动轨迹Fig.15 Whether cavitation effect is considered, centroid trajectory of projectile of wave crest entering water

图16 是否考虑空化效应，弹体波峰入水其质心轨迹在XOY平面投影Fig.16 Whether cavitation effect is considered, projection of centroid trajectory of projectile of wave crest entering water on XOY plane

仿真计算时保存弹体三个方向的受Fx，Fy和Fz，以及力矩中心在原点三个坐标方向的力矩Mx，My和Mz，输出的力和力矩是关于全局坐标系即地面坐标系下的值，需要变换到弹体速度坐标系计算其阻力、升力和侧向力，变换到弹体坐标系计算滚转力矩、偏航力矩和俯仰力矩。

过弹体质心全局坐标系与弹体坐标系的变换矩阵[24]为

(14)

式中，γ，ψ和ϑ为弹体的三个欧拉角，即绕X轴的滚转运动、绕Y轴的偏航运动和绕Z轴的俯仰运动。

通过全局坐标系下弹体质心受力FG和变换矩阵Ω，可求得弹体坐标系下的受力FB

FB=ΩFG

(15)

将弹体坐标系下的受力FB向弹体速度坐标系变换可得弹体运动过程中受到的阻力、升力和侧向力。将弹体坐标系下的载荷用0.5ρV2S进行无量纲处理可得运动过程中弹体的阻力系数、升力系数和侧向力系数，其中：ρ为水的密度；V为弹体初速；S为弹体最大横截面积。

全局坐标系下力矩中心在原点时弹体受到力矩Mo=[MxMyMz]T，全局坐标系下过弹体质心的弹体力矩为MG，过质心弹体坐标系下的力矩为MB=[MrollMyawMpitch]T，其中，Mroll，Myaw和Mpitch分别为弹体的滚转力矩、偏航力矩和俯仰力矩。通过过原点全局坐标系下的力矩可得过质心全局坐标系下的力矩，然后通过变换矩阵可得弹体坐标系下的力矩，即

MG=Mo-RoG×FG
MB=ΩMG

(16)

式中,RoG为全局坐标系中坐标原点到弹体质心矢量。将力矩用0.5ρV2SLt无量纲处理可得力矩系数。

图17～图22分别给出了是否考虑空化效应，波峰入水弹体阻力系数、升力系数、侧向力系数、滚转力矩系数、偏航力矩系数和俯仰力矩系数随时间变化曲线。从图17～图22中总体可以看出，未考虑空化效应弹体入水后直接与水接触，其受到的水动力载荷幅值整体大于考虑空化效应的仿真结果，且在数值仿真计算12 ms内未考虑空化效应载荷幅值变化较为剧烈。这是由于弹体以相同速度在水中航行和在水蒸气与空气相混合气体中航行所受的流体动力载荷幅值有非常大的差别，其在水蒸气与空气相混合气体中所受载荷较小，特别是阻力，极大的提高了射弹在水下航行速度，受到国内外水动力领域的广泛关注和重视。综合本节是否考虑空化效应对射弹高速波峰入水不同相体积分数云图、运动轨迹以及载荷特征分析可知，空化效应是射弹高速入水过程一个非常重要的物理现象，对射弹高速入水过程的运动轨迹以及载荷有非常重要的影响，故在后续射弹高速入水数值仿真过程均考虑了空化效应。

图17 是否考虑空化效应，弹体阻力系数变化曲线Fig.17 Whether cavitation effect is considered, curves of drag coefficient of projectile

图18 是否考虑空化效应，弹体升力系数变化曲线Fig.18 Whether cavitation effect is considered, curves of lift coefficient of projectile

图19 是否考虑空化效应，弹体侧向力系数变化曲线Fig.19 Whether cavitation effect is considered, curves of lateral force coefficient of projectile

图20 是否考虑空化效应，弹体滚转力矩系数变化曲线Fig.20 Whether cavitation effect is considered, curves of rolling moment coefficient of projectile

图21 是否考虑空化效应，弹体偏航力矩系数变化曲线Fig.21 Whether cavitation effect is considered, curves of yawing moment coefficient of projectile

图22 是否考虑空化效应，弹体俯仰力矩系数变化曲线Fig.22 Whether cavitation effect is considered, curves of pitch moment coefficient of projectile

4.3 不同波面位置入水空化特性分析

图23分别给出了波峰、波节1、波谷和波节2四种工况下弹体入水过程不同相体积分数云图，图中浅灰色为水相、深灰色为空气相、白色为水蒸气相。从图23中可知：由于波面的存在，弹体不同位置入水的有效入水角(弹体入水位置处波面的切线与弹体轴向的夹锐角为有效入水角；弹体波峰、波节1、波谷与波节2入水的有效入水角分别为25.003 2°，39.662 0°，25.009 8°与10.309 4°)发生了变化，与名义入水角相比，波节2有效入水角减小了14.690 6°，而波节1有效入水角增加了14.662 0°，波节2工况弹体失稳最快。

图23 高速射弹波面不同位置入水Z=0平面不同相体积分数云图(图中浅灰色为水相,深灰色为空气相,白色为水蒸气相)Fig.23 Contour of different phases volume fraction of high-speed projectile entry water at different positions of wave surface in Z=0 plane (in the figure, light gray is water phase, dark gray is air phase and white is water-vapor phase)

图24分别给出了波峰、波节1、波谷和波节2四种工况下弹体的空泡形态演化过程，其中空泡界面取水蒸气相体积分数为0.5的等值面。从图23和图24中可以看出，弹体入水经历了撞击、空泡形成、超空泡航行和沾湿阶段；对比不同波面位置入水弹体产生的空泡形态变化过程可以发现，波节2工况弹体姿态最先发生翻转，波节1工况弹体下表面沾湿和姿态变化最慢，相同时间内弹体姿态发生很大变化，由大到小依次是波节2、波峰、波谷和波节1。

图24 高速射弹波面不同位置入水空泡形态图Fig.24 Cavitation shape map of high-speed projectile entry water at different positions of wave surface

4.4 波面不同位置入水弹道特性分析

图25和图26分别给出了波峰、波节1、波谷和波节2入水弹体质心在三维空间中的运动轨迹和弹体质心运动轨迹在XOY平面上的投影。从图25和图26中可以看出，入水初期弹体质心运动轨迹保持在弹体初始对称轴方向，随后发生偏转；波节1工况弹体质心轨迹侧向偏移和纵对称面内的偏转均很小，而其他三种工况在入水达到一定深度后弹体质心位移变化较大；四种入水工况弹体侧向偏转程度由大到小依次为波节2、波峰、波谷和波节1，图中仿真结束时波节2工况侧偏了0.321 6 m、波峰工况侧偏了0.061 6 m、波谷工况侧偏了0.006 17 m、波节1工况侧偏了0.000 638 m。

图27和图28分别给出了弹体运动过程中X向、Y向、Z向的速度衰减曲线和总速度大小衰减曲线，其中速度分量的正负与坐标轴方向一致。从图27和图28中可以看出，入水前期弹体基本处于超空泡航行，各方向速度变化基本一致，且侧向Z向速度分量近似为0；入水中期不同工况弹体沾湿出现的时间和变化趋势不同，速度变化量由大到小依次是波节2、波峰、波谷、波节1，在此阶段逐渐产生了侧向Z向速度分量；波节1工况弹体速度变化相对于其他三种工况得到了很大的延缓，同一时刻该工况弹体轨迹和速度变化量均小于其他工况。

图25 高速射弹波面不同位置入水，弹体质心运动轨迹图Fig.25 Centroid trajectory of projectile entry water at different positions of wave surface

图26 高速射弹波面不同位置入水，弹体质心运动轨迹在XOY平面投影Fig.26 Projection of centroid trajectory of projectile on XOY plane entry water at different positions of wave surface

图27 高速射弹波面不同位置入水，弹体质心在X，Y，Z方向速度分量变化曲线Fig.27 Velocity component curves of projectile centroid in X,Y,Z direction entry water at different positions of wave surface

图28 高速射弹波面不同位置入水，弹体质心速度变化曲线Fig.28 Curves of centroid velocity of projectile entry water at different positions of wave surface

图29～图31分别给出了四种工况弹体入水运动过程中的滚转角、偏航角与俯仰角变化规律。从图29～图31中可以看出：弹体滚转角持续近似线性增长，不同工况弹体滚转角的变化在入水前期几乎一样，入水后期产生少许差别；在计算时间内弹体偏航角逐渐增大，偏航角增大量由大到小依次为波节2、波峰、波谷和波节1，图中仿真结束时刻波节2偏航角增大了36.062 2°、波峰增大了16.144 8°、波谷增大了5.148 7°、波节1增大了2.979 2°；弹体俯仰角在弹体接触水面前近似不变，接触水面后俯仰角逐渐减小，在相同时间内俯仰角减小量由大到小依次为波节2、波峰、波谷、波节1，且波节2工况弹体俯仰角变化最快；仿真结束时波节2工况弹体俯仰角减小了75.877 8°，弹体由入水姿态变为出水姿态如图23(d)和图24(d)时刻23 ms所示，波峰工况弹体俯仰角减小了17.433 9°，波谷工况弹体俯仰角减小了6.185 4°，波节1工况弹体俯仰角减小了3.005 9°。

图32～图34分别给出了四种工况弹体入水运动过程的滚转角速度、偏航角速度和俯仰角速度变化曲线。从图32～图34中可以看出，弹体滚转角速度逐渐减小，减小程度由大到小依次为波节2、波峰、波谷、波节1，仿真结束时刻波节2工况滚转角速度减小了142.586 6 rad/s，波峰工况减小了27.884 1 rad/s，波谷工况减小了5.496 3 rad/s，波节1工况减小了2.390 2 rad/s；弹体偏航角速度逐渐增大，仿真结束时刻波节2工况出现了极大值；同时可以看出弹体撞击水面产生的俯仰角速度扰动比较明显，而滚转角速度和偏航角速度几乎没有受到影响。

图29 高速射弹波面不同位置入水，弹体滚转角变化曲线Fig.29 Curves of roll angle of projectile entry water at different positions of wave surface

图30 高速射弹波面不同位置入水，弹体偏航角变化曲线Fig.30 Curves of yaw angle of projectile entry water at different positions of wave surface

图31 高速射弹波面不同位置入水，弹体俯仰角变化曲线Fig.31 Curve of pitch angle of projectile entry water at different positions of wave surface

图32 高速射弹波面不同位置入水，弹体滚转角速度变化曲线Fig.32 Curves of rolling angular velocity of projectile entry water at different positions of wave surface

图33 高速射弹波面不同位置入水，弹体偏航角速度变化曲线Fig.33 Curves of yaw angular velocity of projectile entry water at different positions of wave surface

图34 高速射弹波面不同位置入水，弹体俯仰角速度变化曲线Fig.34 Curves of pitch angular velocity of projectile entry water at different positions of wave surface

4.5 波面不同位置入水载荷特性分析

图35～图37分别给出了四种工况弹体入水运动过程的阻力系数、升力系数和侧向力系数的变化曲线。撞击水面时弹体阻力系数产生瞬间的峰值，随后趋于平缓；入水后期由于弹体表面的沾湿，弹体阻力系数迅速增大，其中波节2工况阻力系数最先增大，这与图23和图24波节2工况弹体最先失稳结论一致；撞击水面时弹体升力系数产生瞬态峰值后趋于平缓，入水后期由于弹体下表面的沾湿，弹体升力迅速增大，由于波节2工况弹体最先沾湿，故其升力系数最先增大；撞击水面时弹体侧向力系数不受影响为0，入水后期弹体侧向力系数逐渐增大，由大到小依次是波节2、波峰、波谷和波节1，这与图25所示弹体侧向位移结论是一致的。

图35 高速射弹波面不同位置入水，弹体阻力系数变化曲线Fig.35 Curves of drag coefficient of projectile entry water at different positions of wave surface

图36 高速射弹波面不同位置入水，弹体升力系数变化曲线Fig.36 Curves of lift coefficient of projectile entry water at different positions of wave surface

图37 高速射弹波面不同位置入水，弹体侧向力系数变化曲线Fig.37 Curves of lateral force coefficient of projectile entry water at different positions of wave surface

图38～图40分别给出了四种工况弹体入水运动过程的滚转力矩系数、偏航力矩系数和俯仰力矩系数的变化曲线。在弹体带空泡运动初期，弹体的力矩系数趋于平缓且量很小；弹体表面沾湿后，弹体力矩系数发生很大的变化，由于弹体旋转运动的存在，弹体的偏航力矩系数和俯仰力矩系数呈现正负交替变化。

图38 高速射弹波面不同位置入水，弹体滚转力矩系数变化曲线Fig.38 Curves of rolling moment coefficient of projectile entry water at different positions of wave surface

图39 高速射弹波面不同位置入水，弹体偏航力矩系数变化曲线Fig.39 Curves of yawing moment coefficient of projectile entry water at different positions of wave surface

图40 高速射弹波面不同位置入水，弹体俯仰力矩系数变化曲线Fig.40 Curves of pitch moment coefficient of projectile entry water at different positions of wave surface

4.6 高速射弹入水不同阶段的载荷与弹道特征

不同学者对于运动体入水各阶段的划分不尽相同，May[25]将整个过程分为六个阶段：冲击阶段、流动形成阶段、空泡敞开阶段、空泡闭合阶段、空泡溃灭阶段以及完全沾湿阶段。本文高速射弹波面入水仅研究前五个阶段，由图23可知高速射弹波面不同位置入水四个阶段(这里将空泡闭合与溃灭合为一个阶段)时间划分如表3所示。

文献[26]指出：当t≪V/g时，重力在许多高速入水问题中通常被忽略。高速射弹波面不同位置入水，仅分析前23 ms，由0.023 s≪(409/9.80)s=41.735 s，故在高速射弹波面入水整个过程中重力因素可以被忽略。事实上在本文数值仿真过程中均考虑了重力效应。

通过图25～图40可知弹体冲击和流动形成这两个阶段，弹体所受重要载荷以及由此引起弹道变化如图41所示。这两个阶段弹体阻力急剧增大产生局部极大值，然后减小趋于稳定值如图41(a)所示，且波面不同位置入水四种工况稳定值几乎完全相同，约为44.73 kN；这两个阶段产生了冲击升力载荷然后趋于稳定值，如图41(b)所示；这两个阶段使弹体产生了初始俯仰角速度扰动如图41(e)所示。

表3 高速射弹波面不同位置入水四个阶段时间划分Tab.3 Time division of four stages of high-speed projectile entry water at different positions of wave surface ms

图41 高速射弹波面不同位置入水，冲击和流动形成阶段弹体载荷与弹道特征Fig.41 Load and ballistic characteristics of projectile entry water at different positions of wave surface in the formation stage of impact and flow

高速射弹波面不同位置入水空泡敞开阶段、空泡闭合阶段以及空泡溃灭阶段，弹体所受载荷与弹体特征详见4.4节和4.3节。

高速射弹波面不同位置入水弹体所受的波浪载荷已包含于其水动力载荷中，由式(12)二阶Stokes波的最大横向速度和垂向速度分别为3.08 m/s和3.06 m/s，远小于射弹入水速度409 m/s；同时仅分析高速射弹入水前23 ms，远小于二阶Stokes波的周期5.06 s。故高速射弹入水过程中波浪载荷对弹体运动轨迹影响不大，可以忽略。

5 结 论

本文针对某一型舰载射弹在波峰、波节1、波谷及波节2下高速旋转入水过程进行了数值模拟，分析了不同波面位置下弹体入水过程中的空泡形态、弹道特征和流体动力载荷特征的变化规律，得到如下结论：

(1) 波面不同位置入水，入水初期产生覆盖弹身的超空泡，空泡不对称，随着入水深度的增加，空泡的不对称性逐渐减小；入水中后期弹体尾部下垂，刺破空泡壁面，空泡形态产生变化，波节2工况弹体姿态最先发生翻转，波节1工况弹体下表面沾湿和姿态变化最慢，相同时间内弹体姿态发生很大变化依次是波节2、波峰、波谷和波节1。

(2) 波面不同位置入水，入水后期弹体质心均会产生不同程度的侧向Z向偏移；仿真结束时刻弹体侧向位移、弹体质心速度变化量、姿态角变化量和角速度变化量由大到小工况依次为波节2、波峰、波谷和波节1。

(3) 撞击水面时，弹体阻力和升力系数产生瞬间峰值随后趋于平缓；入水后期由于弹体表面的大量沾湿，阻力和升力系数迅速增大。入水初期弹体在空泡中航行，弹体力矩系数趋于平缓且量较小；由于弹体旋转运动的存在弹体偏航力矩和俯仰力矩系数呈现正负交替变化，仿真结束时弹体偏航力矩系数、俯仰力矩系数的量级从大到小依次是波节2、波峰、波谷和波节1。

(4) 弹体不同波面位置入水，实际上改变了弹体入水的有效入水角，虽然四种工况弹体名义入水角均为25°，但是波峰、波节1、波谷与波节2入水弹体的有效入水角分别为25.003 2°，39.662 0°，25.009 8°与10.309 4°，因此在相同仿真时间内弹体最先发生失稳的工况依次为波节2、波峰、波谷和波节1。