二级轻气炮内弹道过程数学建模及数值仿真

庄 宇，陆 欣

(南京理工大学 能源与动力工程学院，南京 210094)

【装备理论与装备技术】

二级轻气炮内弹道过程数学建模及数值仿真

庄 宇，陆 欣

(南京理工大学 能源与动力工程学院，南京 210094)

根据二级轻气炮的发射特点，采用经典内弹道模型描述药室里火药燃烧状况和活塞运动，同时采用一维非定常可压缩流动模型描述轻气室里的气体流动状态和弹丸运动，并通过活塞的运动状态将两者耦合，从而建立起二级轻气炮发射过程的内弹道数学模型；以某30 mm/120 mm轻气炮为基本计算模型，运用四阶Runge-Kutta法求解药室方程，运用二阶MacCormack格式求解轻气室方程，通过两部分的交替计算，实现二级轻气炮内弹道过程的数值仿真，为二级轻气炮的参数设计和发射性能的提高提供了理论依据。

二级轻气炮；数学建模；数值仿真；发射性能

二战后，人们对于超高速发射的研究越来越感兴趣。对传统火炮而言，由于火药气体分子量较大，滞止声速较小，弹丸获得的最大速度受到限制，很难达到3 km/s以上。轻气炮是一种利用热轻质气体(如氢和氦)膨胀做功的方式推动弹丸，从而增大气体逃逸速度，减小气体声惯性，使之获得极高速度的发射装置。在1946年世界上第一门轻气炮诞生后的几十年里，这种高温低分子量气体的发射系统在实验室中进行的高速发射研究已经取得了辉煌的成果[1]。近年来，随着轻气炮发射技术不断发展，已经可以发射各种形状的弹丸，而且弹丸的质量、尺寸和材料都有较宽的范围。同时，轻气炮应用领域也不断拓宽，已经从军工和航天扩展到了气动力、超高速碰撞、材料力学性能等方面[2-4]。

目前，国内外最为常见的是二级轻气炮。由于二级轻气炮涉及火药气体、活塞、轻质气体以及弹丸这四者之间的相互作用，发射过程较为复杂。而结构尺寸与装填条件的相互影响和共同作用也使得每门二级轻气炮都存与自身最匹配的发射参数[5]。多年来，国内外各研究单位对此进行了大量的试验和研究。美国的Douglas公司、GM通用汽车公司、NOL海军军械研究室等单位开展的超高速撞击试验为美国的军事及民用方面的发展做出了巨大贡献。其中通用汽车公司(GM)甚至将0.053 3 g的弹丸加速到10.80 km/s的惊人速度[6]。但如果仅仅通过反复试验来确定每门炮的最佳发射参数，显然事倍功半。因此，从减少实验开支等看，辅以必要的二级轻气炮的数值仿真研究显得尤为重要。

为此，针对国内外最常见、应用最为广泛的一种二级轻气炮，运用经典内弹道理论和气体动力学模型，建立二级轻气炮内弹道过程的数学模型，并提出一种行之有效的求解方法。

1 二级轻气炮数学模型

图1是一个常见的二级轻气炮的结构示意图。它的第一级工质是火药气体，第二级工质是轻质气体，一二级之间通过一个自由活塞进行能量传递。其发射过程始于固体火药的点火和燃烧，火药燃烧产生的气体驱动活塞运动，从而急剧压缩泵管中的轻质气体达到高压高温，直到用于隔离弹丸和轻气的膜片破裂，使得弹丸从静止开始加速穿过发射管[7]。需要注意的是，弹丸运动后，活塞会继续压缩轻质气体，直至进入高压段开始变形后，由于受到锥形壁摩擦和挤压以及前端气体压力所形成的阻滞力而不断减速至运动停止(如二级轻气炮参数匹配得当，可使活塞恰好停在锥段附近)。

图1 二级轻气炮结构示意图

考虑上述发射过程，药室和泵管的求解方程应该单独建立并通过活塞的运动耦合，从而建立起二级轻气炮的数学模型。本文用一组常微分方程，即经典内弹道模型描述药室内火药燃烧和活塞的运动规律，同时用一维非定常可压缩流体模型描述泵管内轻气的压缩过程和弹丸的运动规律。鉴于这样的分析，提出如下假设：火药燃烧服从几何燃烧定律；轻气视为恒定比热容的无黏气体，且为一维非定常变截面流动；活塞运动终了时，在缓冲限制器的作用下速度为零；热损失通过火药力f和绝热指数k修正，不做直接计算；其他假设等同于经典内弹道假设。

1.1 药室里的基本方程

1) 火药颗粒形状函数：

(1)

式(1)中，z和ψ分别是火药颗粒的相对燃烧厚度和相对燃烧体积，χ、λ、μ、χs、λs和Zk都是火药的形状特征量。

2) 火药颗粒燃速方程：

(2)

式(2)中，ps是药室的平均压力，u1和n分别是火药的燃速系数和燃速指数，e1的值是火药颗粒初始厚度的一半。

3) 活塞运动方程：

(3)

式(3)中，mp、vp和φ分别是活塞的质量、速度和次要功系数，ppf是作用在活塞面上的氦气压力，S是活塞的横截面积。

4) 活塞速度方程：

(4)

式(4)中，lp是活塞的行程。

5) 能量守恒方程：

(5)

式(5)中，f和ω分别是火药力和火药的装药量，θ是与火药气体相关的常数，lψ是药室自由容积缩径长。

1.2 轻气室中的控制方程

1) 连续性方程：

(6)

式(6)中，ρ和u是分别是泵管里氦气的密度和速度，A是泵管的横截面积。

2) 动量方程：

(7)

式(7)中，P是泵管中氦气的压力，D和L分别是轻气室的直径和长度，cf和k分别是Darcy阻力系数和压力损失因子[7]。其中，cf可用以下关系式表示：

(8)

式(8)中，Re为雷诺数，Λ是对于高速流动十分重要的可压缩修正因子，可用半经验公式表示

(9)

其中，Ma为当地流动的马赫数，Ω为恢复因子。

除非管道面积变化很快，或在流动中隔板阻碍相当严重，压力损失才比较明显。通常情况下，式(7)中的k都比较小，可忽略不计。

3) 能量方程：

(10)

式(10)中，e是氦气的内能，λh和Nu是氦气的导热系数和努塞尔数，T和Tw分别是氦气的温度和壁温。

4) 轻气的状态方程：

(11)

式(11)中，γ和α分别是氦气的比热容比和余容。

2 数学模型的求解

2.1 数值计算方法

药室的基本方程是由常微分方程和代数方程组成，因此可采用4阶Runge-Kutta法数值求解。经典Runge-Kutta法与药室基本方程相容，解无条件收敛和稳定[8]，因而时间步长Δt仅由轻气室所采用的差分格式的稳定性来确定。

轻气室的控制方程是由一组偏微分方程和一个代数方程组成，一般采用有限差分法进行数值求解。通过对几种常见差分格式的比较，查阅相关资料，本文决定选用2阶MacCormack 预估-校正格式来求解二级轻气炮轻气室内的一维非定常可压缩流体模型。此差分格式具有时间上的一阶精度和空间上的二阶精度，整体计算精度较高。

2.2 边界条件处理

通过活塞的运动，将药室里的基本方程和轻气室中的控制方程的求解耦合到一起。一方面药室方程组的计算结果(活塞速度vp和行程lp)可为轻气室提供左边界条件，另一方面轻气室方程组的计算结果又为活塞的运动提供阻力，二者交替进行，直至整个射击过程完成为止，从而实现二级轻气炮发射过程的数值仿真。

计算过程中，轻气室右边界速度和位移分别由弹丸的速度和行程决定：

(12)

式(12)中，S2、mj、φ2分别为发射管的横截面积、弹丸质量和弹丸阻力系数，Pd为弹底压力，取轻气室右边界压力。

轻气室密度、压力、温度等其他运动边界条件可直接由计算域内邻近节点的相应参数经过插值得到。而轻气室的固壁边界条件则需采用第二网格系中的反射法获得，具体法则如下[9-10]：

(13)

2.3 网格处理

2.3.1 活塞顶端网格处理(轻气室左端)

在某个时刻活塞运动到某一位置，以活塞所处的的第k个网格的步长作为变步长进行计算。在下一个时刻活塞运动到下一个位置，此时与前面一个网格点(即第(k+1)个网格点)的距离为Δx1。若Δx1>0.5Δx,则活塞顶端仍与第k个网格联系，并以此新的网格距离Δx1作为新步长进行计算；若Δx1<0.5Δx，则将第k个网格与第(k+1)个网格合并，同时把活塞顶端各参数赋予给新合并的(k+1)网格。新网格步长为(Δx1+Δx)，并且下一时刻的计算就从第(k+1)个网格开始。随着活塞的不断向前运动，左端的网格不断减少，直到活塞运动到泵管高压段附近，停止运动。

2.3.2 弹丸底部网格处理(轻气室右端)

当弹丸开始运动后，轻气室右边界由固壁边界变为运动边界，此时计算网格的数量要增加。

假设在某一时刻弹丸运动到某一位置处，以弹底与后面第一个格点k的距离作为变步长进行计算。在下一个时刻弹丸运动到下一个位置，此时与前面第k个网格点的距离为Δx2。若Δx2<1.5Δx，则弹底仍与第k个网格联系，并以此新的网格距离Δx2作为新步长进行计算；若Δx2>1.5Δx，则分裂出一个新网格(k+1)，并让此网格与弹底发生联系，此新网格以(Δx2-Δx)作为新步长进行计算。

相对于轻气室左端的网格处理，增加了对于新格点各流动参量的求解。对于新格点的处理，可以采用前后相邻格点线性插值的方法，来求得新网格点气体的各个状态参量。此时，所有网格点的各状态参量就可以求得，从而为下一个时间层的计算做准备。

2.4 截面滤波

特别要指出，在计算轻气室各参数过程中，由于差分格式的色散性，使得波头可能产生振荡，且随着计算的进行，振荡积累到一定程度可能使得数值计算无法继续进行，从而造成停机现象。另外，截面的变化也会使面积的导数产生间断。诸如此类问题，通常需要采用滤波处理。常用的滤波方法是Shumann滤波(主要用于截面变化处)，采用的是加权平均方法，即：

(14)

3 数值仿真结果及其分析

在此基础上，对一个30 mm/120 mm二级轻气炮进行数值仿真，从而检验上述数学模型及其计算方法的可行性。30 mm/120 mm二级轻气炮的主要参数如表1所示。

表1 30 mm/120 mm二级轻气炮主要参数

数值仿真结果如图2～图6所示。结合图2、图3和图4可知，当药室压力达到活塞起动压力后，活塞开始运动，并不断压缩轻质气体。随着火药的不断燃烧，药室内的压力开始逐渐增大。在大约4.6 ms时，药室压力增加到最大值310 MPa。随后，压力开始下降。大约在6 ms时，火药燃烧结束。此后，随着活塞的运动，火药气体空间增加，压力持续下降。

图2 药室压力随时间的变化

图3 火药相对已燃百分比随时间的变化

结合图4和图5可知，在药室压力作用下，活塞速度不断增大至最大速度728 m/s。此后，随着轻气室压力的增大，活塞速度开始减小。当活塞运动到高压锥段附近时，由于锥壁阻力的影响，速度迅速降为零，最终停止运动。

图4 活塞和弹丸的速度随时间的变化

图5 活塞速度随行程的变化

结合图4、图6和图7可以看出，约在22.5 ms时，弹丸开始启动，弹底压力开始迅速上升。在维持一段较为平稳的压力过载后，约在24.8 ms，第一个压力波到达弹底，弹底压力达到第一个最大值347 MPa。之后，由于轻气室流场中各种压力波的运动和叠加，使得弹底压力变化剧烈。随着压力波在活塞和弹丸之间不断传递反射，弹丸不断加速运动，直至在出炮口处弹丸速度达到5 906 m/s。可以预见，如果将燃烧过程、活塞质量、轻气室结构以及弹丸质量之间关系匹配得当，使得弹丸在出炮口有一个大幅值的压力波位于弹底，这对增加弹丸的速度非常有利。

图6 弹丸速度随行程的变化

图7 弹底压力随时间的变化

从图8可以看出，在初始时刻，活塞不断压缩轻质气体，导致轻气室左端压强比右端压强稍大，因为此时压缩波尚未传到右端。当接近高压段时，轻质气体分布在一个相对较小的空间里，压强急剧增大，且在空间上的分布出现极大跳动，产生了激波。随后轻质气体绝热膨胀，推动弹丸向前运动，活塞停止，气体压强逐渐减小。

图8 不同时刻轻气室压强沿轴向分布

结合图2～图8，将其与文献[11]中提到的中国工程物理研究院某30 mm口径的二级轻气炮的实验数据进行比较后发现，该仿真结果较为准确地反映了二级轻气炮发射过程的实际状况。数值仿真的弹丸出炮口速度为5.91 km/s，实验值为5.56 km/s，误差为6%，在允许范围内。这也在一定程度上验证了所建立的二级轻气炮的数学模型和所采用的数值求解方法。

4 结论

本文运用经典内弹道理论和气体动力学模型建立了二级轻气炮内弹道过程的数学模型，给出了一种运用四阶Runge-Kutta法和二阶MacCormack格式分别交替求解药室方程和轻气室方程的数值计算方法，从而实现了二级轻气炮的数值仿真，较为准确地刻画了它的发射过程。由此可见，通过数值仿真来辅助二级轻气炮的实验研究是可行的，可以为二级轻气炮的试验装填参数优化和发射性能的提高提供一定的理论依据。

[1] BERGGREN R E，REYNOLDS R M.The Light-Gas-Gun Model Launcher[R].AD 713915，1970.

[2] 黄洁,梁世昌,李海燕,等.二级轻气炮发射过程内弹道数值计算研究[J].空气动力学学报，2013，31(5):657-661.

[3] 柳森，李毅，黄洁，等．弹丸超高速撞击单层和多层板结构的碎片特征研究[J]．宇航学报，2010，31(6)：1672-1677.

[4] ZENG X J,YU Z F,BU S Q,et al.Research on the RCS of Hypervelocity Model and Its Plasma Sheath[J].ACTA Aerodynamica Sinica,2010,28(6):645-649.

[5] 管小荣，徐诚.二级轻气炮发射过程数学模型和计算方法[J].南京理工大学学报，2007，31(1)：22-26.

[6] 王东方，肖伟科，庞宝君.NASA二级轻气炮设备简介[J].实验流体力学,2014,28(4):99-104.

[7] 金志明，翁春生.高等内弹道学[M].北京：高等教育出版社，2003.

[8] 郑慧娆，陈绍林，莫忠息，等.数值计算方法[M].武汉：武汉大学出版社，2002.

[9] 翁春生，王浩.计算内弹道学[M].北京：国防工业出版社，2006.

[10]黄玲智.具有脉冲扩散和扩散时滞的单种群模型研究[J].重庆工商大学学报(自然科学版),2015,32(9):5-9.

[11]冯福全.二级轻气炮内弹道仿真研究[D].南京：南京理工大学，2008.

[12]黄滔，张相炎，刘宁.燃烧轻气炮点火对内弹道性能影响的仿真分析[J].兵工自动化，2014(12):25-28.

(责任编辑周江川)

Mathematical Modeling and Numerical Simulation for Interior Ballistics Process of Two-Stage Light-Gas Gun

ZHUANG Yu, LU Xin

(School of Energy and Power Engineering,Nanjing University of Science and Technology, Nanjing 210094, China)

According to the launch features of two-stage light gas gun, classical interior ballistic model was used to describe combustion of powder in chamber and piston motion, and one dimensional unsteady compressible fluid model was used to describe flowing of gas in light-gas chamber and projectile motion. The two groups of equations were linked through piston motion to develop a mathematical model for interior ballistics launch process of two-stage light-gas gun. A 30 mm/120 mm light-gas gun was taken for the primary computing model.The equations in powder chamber were solved numerically with the fourth-order Runge-Kutta method, and the control equations in light-gas chamber were solved numerically with the second-order MacCormack scheme.The numerical simulation for interior ballistics process of two-stage light-gas gun was accomplished by alternation between these two computations. The simulation results provide theoretical references for parameter design and improvement of launching performance of two-stage light-gas gun.

two-stage light-gas gun; mathematical modeling; numerical simulation; launching performance

2016-07-27；

江苏省前瞻性联合创新资金(BY2014004-11)

庄宇(1992—)，男，硕士研究生，主要从事新型发射理论与控制技术研究。

10.11809/scbgxb2016.12.018

庄宇，陆欣.二级轻气炮内弹道过程数学建模及数值仿真[J].兵器装备工程学报，2016(12):75-79.

format:ZHUANG Yu, LU Xin.Mathematical Modeling and Numerical Simulation for Interior Ballistics Process of Two-Stage Light-Gas Gun[J].Journal of Ordnance Equipment Engineering,2016(12):75-79.

TJ012.1

A

2096-2304(2016)12-0075-05

修回日期：2016-08-30