旋转射弹高速倾斜入水多相流场与弹道数值模拟*

朱 珠，罗 松，卢丙举，于 勇

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

高速运动体入水过程中，运动体头部会发生相变产生空化，空化蒸汽与空气排开水质点形成入水空泡。空泡的存在使得弹体与水的接触面积大大减小，可以达到减阻的效果，如何形成包裹全弹身的超空泡是超空泡减阻技术的关键。同时，超空泡形态和变化规律是弹体以正确的姿态和较高的末动能达到目标地点的关键因素。对高速弹体入水空化问题进行研究，具有重要意义。

入水问题源于Worthington 等[1-2]在20 世纪初对于球体入水的大量实验，他们在实验中捕捉到了入水溅射和空泡发展。随后，由于水下鱼雷失稳、弹跳等实际问题的出现使得入水问题的研究得到重视，入水问题的研究对象从球体转向了不同头型的圆柱回转体，船舶和实际的一些弹型。俄罗斯和乌克兰学者对于高速射弹的空泡特性研究开展了大量的基础性工作[3-5]。Logvinovich[6]基于独立膨胀原理对入水空泡形态的求解成为之后研究空泡形态和空泡壁发展的理论基础。施红辉[7-8]、宋武超[9]等实验研究了细长体低速入水时入水角度和入水速度等因素对入水空泡形态的影响。马庆鹏[10]对圆柱体垂直入水问题进行了数值模拟研究，详细分析了入水空化过程的流动机理。何春涛[11]详细研究了超空泡射弹的几何参数对空泡形态和阻力特性的影响，为超空泡射弹的设计提供了参考。随后，孙健[12]、宋武超等[13]对小型回转体倾斜和垂直入水问题进行了三维数值计算，利用六自由度算法和动网格对回转体入水弹道进行了求解。齐亚飞[14]利用实验和数值计算对不同头型的回转圆柱射弹高速入水时的弹道和空泡特性进行了研究，得到了不同头型不同初速度下的弹体衰减规律以及不同头型的入水稳定性规律，同时在实验中确定了临界入水角的存在。赵成功等[15-16]研究了细长体倾斜入水的弹道特性和空泡闭合脱落特性，同时研究了弹体的尾拍运动对弹体稳定性的影响。李佳川[17]利用数值计算模拟了高速射弹入水过程中的初始扰动角速度、弹体密度和入水角度对弹道稳定性的影响。

目前对于入水问题的研究，低速入水问题研究的比较多，高速入水问题由于条件的限制，实验研究和数值计算都很少，特别是弹体高速旋转倾斜入水空化和弹道稳定性问题的研究，有待进一步深入。本文使用基于某一型号舰载弹体改进后的弹型作为研究对象，利用Fluent 平台进行计算，考虑自然空化，利用six DOF 算法计算弹体姿态轨迹，研究不同入水角度下高速旋转射弹入水过程中的空化特性、弹道特性和流体动力特性的变化规律。

1 数值计算方法

1.1 流场求解与弹道的耦合

计算采用ANSYS Fluent 流体力学计算软件，采用有限体积法求解N-S 方程组，多相流场的描述采用均质多相流模型中的VOF 模型[18]，湍流流场的描述采用SSTk-ω 湍流模型[19]。同时，本文选取基于Rayleigh-Plesset 气泡方程[20]建立的Schnerr and Sauer 空化模型来描述高速射弹入水过程中的自然空化现象。计算过程为瞬态计算，考虑重力的作用。弹体弹道的求解采用刚体运动的动力学方程和运动学方程[17]。在每个时间步，Fluent 会计算出弹体受到的流体动力及力矩，根据动力学方程求得弹体的平移加速度和旋转角加速度，之后积分得到弹体的速度和角速度，进一步计算得到弹体新的质心位置和欧拉角，更新网格之后进行下一时间步的求解，直至求解结束。

1.2 计算模型与边界条件

本文计算采用的射弹模型是如图1 所示的实心回转体，弹体模型由头部锥台段、曲线弹身段、圆柱段及后体斜台段组成，弹体总长Lt，头部空化器直径为Dn，弹体最大直径为Dt，弹体长细比为Lt/Dt=8.746。图2 展示了计算域对称面内示意图以及边界条件。计算域长20 m，高6.5 m，宽4 m，水深4 m，空气域2.5 m，重力沿y轴负方向，弹体初始质心位置位于(8 m, 0 m, 0 m)处。初始时刻弹体对称轴与水面夹角为α，弹体初速度408 m/s，对应自然空化数为0.001 177，弹体绕着其对称轴高速旋转，旋转速度为7 500 r/min。计算域x+,x-面为滑移壁面，z+,z-面为symmetry 边界条件，计算域上方为压力入口，压强与标准大气压一致，底部为压力出口，压强为对应4 m 水深的静水压强。

图 1 弹体对称面结构示意图Fig. 1 Sketch of projectile’s symmetry plane

图 2 弹体倾斜入水流域对称面示意图Fig. 2 Sketch of the fluid domain’s symmetry plane for oblique water entry of projectile

1.3 网格划分

计算中采用了Overset 网格技术，划分网格时分为部件网格和背景网格。同时导入背景网格和部件网格进行计算，计算时Fluent 会对背景网格进行挖洞操作，然后嵌套部件网格，去除Dead Cell 后进行计算。流场参数会在部件网格和背景网格重叠边界处进行插值计算。计算网格采用结构化网格划分方式，部件网格区域为一包裹弹体的圆柱，背景网格区域是2.1 节中的长方体计算域。图3 显示了30°入水角时初始时刻计算采用的进行挖补操作之后的网格在z=0 平面上的分布。

图 3 30°入水角z=0 平面网格示意图Fig. 3 Schematic of the grid in z=0 plane of 30 degrees inclined angle

采用Overset 网格技术，设置部件网格计算域与运动体同步运动，结合six DOF 算法可以实现运动体运动的六自由度运动解算。重叠网格技术计算中，只涉及部件网格和不同区域的背景网格的挖补插值操作，不会产生网格的变形，因此可以保持整个计算过程中都有质量很高的计算网格。

1.4 数值方法与设置

计算过程中，体积分数的求解采用隐式格式(Implicit)，相间交界面的构建采用Compressive 格式，计算中设置水为主相，水蒸气和空气为第二相，相间相互作用选取Schnerr and Sauer 空化模型，入水湍流流场的模拟采用选取SSTk-ω 湍流模型；时间离散采用一阶精度，空间离散采用二阶精度；瞬时速度场和压力场的耦合计算采用Coupled 算法，压力差值采用PRESTO 格式，其余项采用二阶迎风格式，计算采用的时间步长为10-6s。

1.5 数值方法验证

Lee[21]利用能量守恒原理得到了入水空泡的计算方法。忽略弹体入水产生的粘性热，考虑流体不可压，对弹体运用牛顿第二定律可以得到弹体的入水速度方程

为了验证数值计算的可靠性，首先对郭子涛实验研究[22]的平头圆柱垂直入水问题进行数值模拟计算，实验采用的平头圆柱长25.4 mm，直径为12.56 mm，重25.1 g，初速度为603 m/s。计算采用六自由度算法和重叠网格技术。图4 显示了圆柱体入水过程中的速度衰减曲线和入水位移曲线，将仿真结果与实验结果及采用以上弹道理论解获得的结果进行了对比，可以看出，数值计算结果与实验结果和理论计算结果均具有较好的一致性，说明数值计算方法是可信的。

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

2 旋转效应对射弹入水过程的影响

超空泡射弹发射后一般都具有沿自身轴线的角速度，射弹在高速旋转入水及水下运动时，其弹道稳定性和流体动力特性都会有一定的特殊性。为了得到射弹的旋转对超空泡射弹入水过程的影响规律，这里选取图1 中的弹体在初速度408 m/s，入水角度α=20°，初始旋转速度ω0分别为0 和7 500 r/min 的工况进行数值模拟，将得到的结果记性了对比。弹体的空化器直径为Dn=35 mm，质心到头部的距离是弹体总长的60.21%。

图5 显示了射弹在无初始旋转速度和高速旋转状态下的质心在计算流域内的三维轨迹变化曲线。在相同的入水条件下，射弹旋转使得弹体质心的侧向偏移增大，为了进一步的得到弹体在初始对称平面内的弹道偏转情况，图6 显示了质心轨迹在xOy平面上的投影，可以看出在高速旋转状态下弹体质心在xOy平面内偏离初始运动方向的程度大幅度减小，这说明弹体的旋转使得弹体在初始对称平面内的弹道稳定性得到了提高。

图 5 不同旋转速度下弹体质心轨迹图Fig. 5 Centroid trajectory of projectile at different rotation velocities

图 6 不同旋转速度下弹体质心轨迹xOy 平面投影曲线Fig. 6 Projection of centroid trajectory of projectile on xOy plane at different rotation velocities

图7 显示了不同旋转速度下的弹体质心速度的变化曲线，可以看出，高速旋转状态下的弹体在入水过程中的质心速度衰减明显小于无初始旋转的弹体。弹体高速旋转使得弹体的阻力减小，末动能可以得到很大的提高。图8 显示了两种初始旋转状态下的弹体入水过程中的偏航角和俯仰角变化曲线，从图中可以看出，弹体旋转状态下入水俯仰角变化大幅变小，但偏航角变化明显大于无初始旋转速度入水工况。这说明弹体的旋转对弹体水下运动的姿态有着重要的影响，特别是在弹体表面有部分沾湿时，弹体的旋转使得弹体的俯仰姿态的稳定性得到提高，但是偏航姿态的稳定性减弱。这是由于弹体旋转过程中，由于弹体表面的不对称沾湿产生不对称的压力分布，从而产生侧向的作用力，这作用力力会产生较大的偏航力矩。而旋转状态使得弹体的沾湿区域相对于无旋转或近似无旋转状态有所变小，这使得弹体的俯仰力矩减小。这将在流体动力特性曲线中得到证实。

图 7 不同旋转速度下弹体质心速度变化曲线Fig. 7 Velocity curves of projectile at different rotation velocities

图 8 不同旋转速度下弹体姿态角变化曲线Fig. 8 Attitudes change of projectile at different rotation velocities

图9 和图10 分别显示了两种初始旋转状态下的弹体阻力系数和升力系数的变化规律，从图中可以看出在入水运动初期，两种状态的弹体阻力系数和升力系数基本没有差别，此时弹体处于超空泡航行阶段，弹体仅头部沾湿，其余部分被空化蒸汽包裹；在入水运动后期，高速旋转的弹体阻力系数和升力系数均小于无旋转弹体，此时弹体前弹身下表面部分沾湿，弹体高速旋转改变了弹体沾湿状态时的受力。

图 9 不同旋转速度下弹体阻力系数(Cd)变化曲线Fig. 9 Drag coefficients (Cd) evolution of projectile at different rotation velocities

图 10 不同旋转速度下弹体升力(Cl)系数变化曲线Fig. 10 Lift coefficients (Cl) evolution of projectile at different rotation velocities

图11 和图12 显示了两种初始旋转状态下的弹体偏航力矩系数和俯仰力矩系数的变化曲线，从图中可以看出，在入水运动前期，弹体处于超空泡航行阶段时，旋转对弹体的偏航力矩系数和俯仰力矩系数影响很小；在入水运动后期，弹体表面沾湿时，高速旋转状态的弹体受到的偏航力矩系数迅速增大，无旋转弹体的偏航力矩系数一直近似为0，高速旋转状态的弹体俯仰力矩系数的量级远小于无旋转状态的弹体。

图 11 不同旋转速度下弹体偏航力矩系数变化曲线Fig. 11 Yawing moment coefficients evolution of projectile at different rotation velocities

图 12 不同旋转速度下弹体俯仰力矩系数变化曲线Fig. 12 Pitching moment coefficients evolution of projectile at different rotation velocities

3 不同角度对射弹入水过程的影响

为了进一步得到入水角度对射弹入水过程影响的规律，对射弹在初速度408 m/s，旋转速度7 500 r/min，入水角度α 分别为20°，25°，30°情形下入水过程中的超空泡形态、弹道稳定性、流体动力特性进行了分析。

3.1 空泡形态特性

图13 显示了不同入水角度不同时刻空泡形态的变化，空泡壁面取水蒸气相体积分数为0.5 的等值面，图中时间从1～19 ms，时间间隔为3 ms。从图中可以看出，弹体迎水面空泡发展快于背水侧，形成的空泡不对称，随着入水深度的增加，这种不对称性逐渐变小，但当弹体运动状态变化使弹身刺破空泡壁时，空泡的对称性被打破，空泡的形态发生较大变化。弹体自然空化剧烈，空泡内充满蒸汽。随着入水深度的增加，自然空化数增大，空泡尺寸有所减小，弹体下表面有所沾湿，弹体姿态开始变化，其中20°入水工况弹体下表面沾湿和姿态变化较明显。同时可以看出，入水角度越小，撞击水面形成的空泡越不对称，开空泡阶段的空泡形态越不对称。

图 13 弹体不同角度入水空泡形态图Fig. 13 Cavitation shapes of the projectile at different inclined water entry angles

通过分析不同入水角度弹体产生的空泡形态变化过程，我们可以得知，在本文计算工况下，弹体能够较为稳定的射入水中，但随后弹体姿态的轻微变化导致弹体发生“尾拍”现象，不同于正常的“尾拍”情形，由于空泡直径偏小，空泡壁面与弹体下表面距离较近，同时弹体曲线段曲率较大，弹体发生“尾拍”时并非是尾部段接触空泡壁，而是弹体的曲线弹身段先刺破空泡壁沾水，沾湿区域会受到极大的摩擦阻力和和压力，而这一位置位于弹体质心的前方，使得弹体受到的力和力矩是不稳定的，导致弹体轨迹和姿态发生很大变化。

3.2 弹道特性

图14 显示了不同入水角度下弹体质心在计算流域内的运动轨迹，可以看出在计算的工况下，入水角度越大，弹体质心轨迹偏离初始弹体对称轴方向的程度越小，不同入水角下弹体都会有一定程度的侧向偏移，入水角越大，侧向偏移量越小，其中20°入水工况下弹体能够稳定入水，但在入水深度达到一定程度后，弹体逐渐失稳，质心轨迹发生很大的偏转。同时，图15 显示了弹体入水时的速度衰减规律，可以看出在入水运动前期，各角度下弹体的速度衰减较为一致，此时弹体仅头部沾湿，阻力仅来源于头部，随后弹身下表面沾湿，弹体的阻力变大，弹体速度衰减变快，不同角度下弹体表面沾湿大小不一样，速度衰减产生差异，其中20°工况速度衰减明显快于25°和30°工况，这是由于20°工况下的弹体沾湿面积扩张的更快，导致阻力增加的更快。通过这里可以看出，弹体下表面的沾湿对弹体的轨迹和速度会产生很大的影响，弹体侧表面的不对称沾湿是弹体轨迹偏转和速度急剧减小的重要原因。

图16、图17 和图18 显示了弹体入水过程中的欧拉角变化，对于计算的几种工况，滚转角由于初始转动角速度的存在持续增大，受到的扰动较小，不同入水角度下滚转角的变化几乎一样，只在入水后期产生轻微差别；弹体运动初期稳定入水，弹体处于超空泡航行阶段，弹体的俯仰角和偏航角初始阶段变化很小，随后由于弹体姿态的轻微变化导致弹体前弹身下表面沾湿，使得弹体产生较大的力矩，弹体的俯仰角和偏航角随之发生很大的变化。同时可以看到，入水角度越小，弹体俯仰角和偏航角的变化会越大，这种扰动变化会随着入水深度的增加逐渐累积变大。其中对于20°工况弹体在入水初期运动稳定，但当其运动到一定深度后，弹体俯仰角发生极大变化，由正变负，弹体逐渐发生失稳，弹体运动姿态由入水姿态变为出水姿态，弹体抬头发生翻转，同时弹体产生很大的偏航角，弹体的运动姿态相对于初始状态发生很大变化。弹体的侧表面的不对称沾湿对弹体的运动姿态产生了极大的影响，是影响弹体姿态稳定性的重要因素。

图 14 不同入水角度下弹体质心轨迹图Fig. 14 Centroid trajectory of projectile at different inclined water entry angles

图 15 不同入水角度下弹体质心速度变化曲线Fig. 15 The velocity curves of projectile at different inclined water entry angles

图 16 不同入水角度下弹体滚转角(φ)变化曲线Fig. 16 Roll angle (φ) curves of projectile at different inclined water entry angles

图 17 不同入水角度下弹体偏航角(θ)变化曲线Fig. 17 Yaw angle (θ) curves of projectile at different inclined water entry angles

图 18 不同入水角度下弹体俯仰角(ψ)变化曲线Fig. 18 Pitch angle (ψ) curves of projectile at different inclined water entry angles

3.3 流体动力特性

定义弹体的流体动力系数如下：

图19 和图20 显示了弹体入水过程中阻力系数和升力系数变化情况，从图中可以看出，弹体撞击水面时，阻力系数和升力系数均在极短的时间内达到峰值，随后随着空泡的形成，弹体的阻力系数和升力系数缓慢下降，随后趋于稳定，稳定运动时的升力系数趋于0，阻力系数很小；但随后弹体沾湿，升阻力系数开始持续增大，从这里可以看出，弹体在姿态发生变化时，表面沾湿形成的升力并未成为恢复升力使弹体恢复到空泡内部运动，弹体沾湿形成的受力使弹体的轨迹变化加剧，沾湿面积增大，从而进一步增大了受力。其中对于20°工况弹体失稳后升阻力变化明显大于25°和30°工况。

图 19 不同入水角度下弹体阻力系数变化曲线Fig. 19 Drag force coefficients evolution of projectile at different inclined water entry angles

各力矩系数的变化规律如图21、图22 和图23 所示。在弹体撞击水面时刻，偏航力矩系数和俯仰力矩系数受到轻微的扰动，滚转力矩系数没有受到扰动，在弹体带空泡运动初期，弹体的力矩系数趋于平缓且量级很小，在弹体表面沾湿后，弹体的力矩系数产生较大变化和波动，其中25°，30°入水角下弹体力矩系数变化趋势较为一致，20°入水角情形下弹体的力矩系数量级较大，因此20°工况下弹体的运动姿态变化比另外两种工况快。滚转力矩系数持续增大，而偏航力矩系数和俯仰力矩系数正负交替，这是由于旋转导致的弹体坐标系方向的改变造成的。

图 20 不同入水角度下弹体升力系数变化曲线Fig. 20 Lift force coefficients evolution of projectile at different inclined water entry angles

图 21 不同入水角度下弹体滚转力矩系数变化曲线Fig. 21 Rolling moment coefficients evolution of projectile at different inclined water entry angles

图 22 不同入水角度下弹体偏航力矩系数变化曲线Fig. 22 Yawing moment coefficients evolution of projectile at different inclined water entry angles

图 23 不同入水角度下弹体俯仰力矩系数变化曲线Fig. 23 Pitching moment coefficients evolution of projectile at different inclined water entry angles

4 结 论

（1）当弹体处于超空泡航行阶段时，弹体的旋转对弹体的弹道姿态，流体动力特性无明显影响；但入水深度达到一定值时，此时形成的超空泡不能完全包裹弹体，弹体的旋转使得弹体弹道的侧向稳定性下降，但初始对称平面内的弹道稳定性得到大幅度的提高，弹体的旋转使得弹体的升阻力系数减小，弹体速度衰减变慢，弹体的旋转使得弹体的俯仰力矩系数大幅变小，姿态角变化相对于无旋转弹体大幅减小，但旋转使得弹体的偏航力矩系数增大，偏航角变化相对于无旋转弹体加剧。（2）在本文计算的倾斜入水角度下，弹体入水初期产生的空泡不对称性较小，弹体产生覆盖弹身的超空泡，入水角度越小，空泡不对称性越强；入水中后期弹体轨迹和姿态发生变化，弹体下表面刺破空泡壁面，空泡形态产生变化，其中20°工况的弹体沾湿和空泡形态变化最明显。（3）在本文的入水角度下，弹体均能稳定入水，但弹体的运动轨迹会产生不同程度的偏移，同时弹体的姿态也会发生不同程度的变化，弹体入水角度越大，弹体质心的偏移量越小；相较于另外两种入水角度，20°入水时弹体入水运动一段时间后发生严重失稳，弹体质心轨迹和弹体姿态产生的变化极大，速度衰减急剧变快。（4）三种角度入水角度下，弹体在下表面未沾湿前受到的力和力矩系数差别不大，但是当弹体下表面沾湿后，入水角度越小，弹体受到的力和力矩会越大，其中20°入水角下弹体的力和力矩系数的量级远大于另外两种入水角度下的力和力矩中时间后弹体质心轨迹和弹体姿态产生的变化极大，弹体发生反转。