超空泡射弹运动过程模拟仿真及有限元分析

张丁雄，王义冬，刘天放

(1.海军装备技术研究所，北京 102442;2.海军航空工程学院，山东 烟台 264000)

超空泡射弹在水下航行时被超空泡包裹着，弹丸尾部与空泡壁发生剧烈的碰撞，极大的轴向冲击载荷和剧烈的尾拍使得超空泡射弹的结构响应与传统的航行体有很大的不同，作用在空化器上的力和周期性作用在弹丸尾部的力有可能使弹丸结构发生强度破坏［1］，因此研究超空泡射弹尾拍飞行时的应力变化情况具有重要的意义。

由于弹丸在水中运动速度非常快，直接对作用在弹丸上的力进行测量非常困难，本文利用Matlab 对弹丸运动过程进行模拟仿真，分析弹丸在水中的受力情况及运动速度，利用仿真得到的尾拍载荷，对做自由飞行的超空泡射弹在尾拍过程中弹丸的各个节点应力特性进行有限元仿真。

1 弹丸水中运动仿真模型

弹丸在水中运动时由于空化器的作用形成了一个包裹弹丸的超空泡，本文采用的空化器为圆盘空化器，且空化器平面与弹丸的轴线垂直，使用的模型如图1 所示［2］。

图1 弹丸模型

已有实验证明当超空泡射弹的运动速度超过8 m/s 后，重力对弹丸的作用效果很小［3］，本文不计作用在超空泡航行体上的重力，这种近似是可以接受的，故假设运动体的质心近似于沿水平直线前进。同时，由于常温下水的饱和蒸汽压相对较低(只有2 350 Pa)，泡内的气体也比较稀薄，为了简化分析，对空泡内的气体摩擦阻力也进行了忽略［4］。

超空泡射弹做尾拍飞行时速度非常大，使弹体前端面上的流体阻力和升力非常巨大，直接导致弹体前端的上、下运动变得非常困难，在弹体同时做前进和转动的刚体运动时，弹丸空化器部分在纵向平面内(偏离弹体水平前进的方向)的位移相对于其他部分弹体位置是极小的，或者说，其转动幅度极小，所以，近似地认为弹体是绕前端面转动的。另外，国外的有关文献中有如下的结论:尾拍冲击会使周围的液体受扰动并发生液体振荡，当振荡波的波长比弹体尺寸大得多时，那么就可以将弹体的前端点近似看作弹体的转动中心［5］，并且AHSUM 的系列射弹实验也证明了这一假设的合理性。

为了简化分析，将弹体模型的空化器中心点看作静止不动［6］，在初始扰动角速度ω 的作用下，弹体绕该空化器中心转动，尾部不断与空泡壁发生碰撞，其中空泡周围的水以变化的速度往弹丸尾部方向流动。即射弹飞行时，弹体与空泡一同前进，由于运动具有相对性，因此可以假设弹体和空泡不做前进运动，而利用空泡周围的水介质向相反方向的运动来代替，将弹丸看作只做一个绕空化器中心转动的运动，其中假设弹丸的头部有一个铰链，弹丸没有了向前的平动速度，只有相同大小，方向相反的来流速度。

本文弹体变形只考虑弹性变形，而不考虑塑性变形，对弹丸进行受力分析，其受力如图2 所示。弹体上的载荷有作用在空化器上的阻力Fn，它的方向为垂直于空化器表面;弹体尾部撞击空泡壁时受到水对弹丸的撞击力Fsy1，方向为垂直于弹丸的轴线;因为有作用在空化器上的阻力的存在，弹丸做减速运动，根据达朗贝尔原理可以知道弹丸受到一个惯性力FT，其中FT=Fn，方向与Fn相反，在这3 个力的作用下弹丸绕空化器中心转动［7］。

图2 弹体上的主要载荷

2 数值仿真

本文采用欧拉算法对弹丸的运动过程进行模拟分析［8］，欧拉法属于简单的数值方法，其基本思路:做等距分割，利用数值微分代替倒数项，建立差分方程。

1)第一阶段(与空泡壁发生撞击)

动力学方程为:

欧拉法可将其写为:

2)第二阶段(未与空泡壁发生撞击)

弹丸方程有:

将其写为欧拉法形式为:

本文提出了一种针对空泡壁为椭圆边的改进方案，过原点的弹体直线方程与椭圆空泡壁方程联立有:

将其化为欧拉算法可以得到

其中:ξ，( )η 为弹体与空泡壁的交点坐标;l(i +1)为弹丸在(i+1)时刻空泡的长轴长度;Dk(i+1)为弹丸在(i +1)时刻空泡的短轴长度;θi+1为弹丸在(i +1)时刻弹丸的弹丸与炮口发射方向夹角。

故浸湿长度为

其中L 为弹丸的长度。

3 仿真结果

利用MATLAB 对改进前后弹丸运动过程进行仿真，仿真结果如图3 所示。

图3 改进前后各个物理量对比

图3 中可以看出，改进后弹丸的衰减速度比原来更快，头部的空化器阻力刚开始没有发生变化，但是随后受到的力比改进前更小了;尾拍力比原来增大。造成这种不同的主要原因是由于改进前将空泡壁看成是一条直线，撞击点的空泡直径看成是弹丸尾部的对应的空泡处直径，且撞击点的纵坐标不再发生变化，但是在真实的情况中，由于空泡壁是椭圆形的，弹丸与空泡壁相交的位置纵坐标是发生变化的，随着撞击的深入，空泡壁是一个弧形，弹丸与空泡发生撞击的分界点处对应的空泡直径实际上是变小的，故实际的浸入长度比未改进前要大，所以弹丸发生撞击时将受到更大的撞击力(即尾拍力)，弹丸速度衰减得更快，故空化器阻力也变小。

弹丸在发射的过程中，由于种种原因，出炮口时初始扰动存在差异，本文假设初始扰度为100 rad/s。针对弹体服役过程，作用于模型中的两个力由前面Matlab 中仿真给出，然后通过ABAQUS 求解器进行完全依赖时间的动态分析［9］，铰链位置采用5 个自由度进行约束，放开一个自由度的方案，在避免约束不足的前提下模拟铰链处的衔接，另外在杆件尾端加载两个集中力和一个初始角速度，在0.25 s 的时间内，加载5002 组集中力，通过ABAQUS 求解器进行分析。仿真得到一组连续变化的应力分布图，图4 为其中4 个时刻的应力分布。

图4 扰动为100 rad/s 时4 个时刻应力分布

图5 为所有节点的应力变化图。图5 中可以看出，弹丸的最大应力在在运动后不久达到，具体时间从ABAQUS 中可知为t=0.021 85 s，对应的单元编号为4 647 和4 648，最大值为σr=2 ×104MPa，此时应力分布如图6 所示。

变换视角，对应力最大值区域进行局部放大，如图7 所示。图7 中可以很直观地观察到弹丸在水中运动时，应力最大值发生在空化器后的一小段区域内，在整个过程中，此处最容易发生强度破坏现象，在图7 中，尾部的应力分布不是很明显，这是由于应力的上下限设置不同，对图片重新处理，得到图8 效果图。

图5 所有节点应力-时间变化

图6 应力最大时刻的应力分布

图7 应力最大时刻的不同视角应力分布

图8 应力上下限重新设置后应力分布

图8 中可以看出，弹丸的尾部应力也比较大。由于本文建模的时候把尾拍力施加在尾部的刚体点上，刚体与尾部一段区域耦合，故在与水撞击的时候，与水接触的那部分弹体的应力分布存在一定的误差;本文由于只考虑弹丸发生了弹性变形而没有发生塑性变形，导致弯曲过大时弹丸依旧没有发生强度破坏，所以仿真弹丸弯曲程度有可能偏离实际情况，存在一定的误差。从仿真结果动态视频中可以看出，弹丸前端的应力相对比较大，弹丸在发射过程中空化器附近区域容易发生强度破坏，导致弹丸变形。

4 结束语

本文利用Matlab 对弹丸的运动的外弹道方程进行模拟仿真，并对比了改进前后的仿真结果，发现改进后的仿真结果更加准确，将仿真得出的空化器阻力、尾拍力和初始角速度加载到弹丸水中运动的模型中进行有限元行仿真，使用ABAQUS 有限元仿真软件，仿真得到弹丸在运动过程中应力分布随时间的变化图和每个节点应力随时间的变化曲线图，可以很直观地看出最大应力随时间的变化关系，从ABAQUS软件中导出最大应力值大小，并且得到最大应力产生的时间和位置，求出此时刻的应力分布图，可以看出最容易发生强度破坏的区域。该区域的强度和刚度都是在设计弹丸时应该重点考虑的部位。

［1］Hrubes J D.High-speed imaging of supercavitating underwater projectiles［J］.Experiments in Fluids，2001(30):57-64.

［2］王献孚.空化泡和超空化泡流动理论及应用［M］.北京:国防工业出版社，2009.

［3］May A.Water Entry and Cavity-Running Behavior of Missiles［R］.［S.l.］:［s.n.］，1975.

［4］刘玉秋，张嘉钟，于开平.非流线型水下航行体阻力测量实验中的几个问题［J］.船舶工程，2006，28(4):39-42.

［5］Jenkins A，Evans L T. Sea mine neutralization using the AN/AWS-2 rapid airborne mine clearance system［J］.IEEE Aerospace Conference Proceedings，2004:2999-3004.

［6］孔德仁，狄长安，朱蕴璞，等.枪弹水中运动阻力软测量方法研究［J］.仪器仪表学报，2006，23(3):286-287.

［7］杨传武，王安稳.超空泡水下航行体振动特性分析［J］.海军工程大学学报，2008，20(4):30-35.

［8］Albert M.Water entry and the cavity-running behavior of missiles［R］.［S.l.］:［s.n.］，1975.

［9］石亦平，周玉蓉.ABAQUS 有限元分析实例详解［M］.北京:机械工业出版社，2008.