防护结构设计是被动防护任务的核心，是航天器空间碎片防护领域研究的重点^[1]。为更好地开展相关研究，各种超高速发射技术应运而生，但大部分发射技术都存在一定的局限性。二级轻气炮因具有众多优点从各类发射技术中脱颖而出，如：适用弹丸种类多；在低加速度和低应力条件下，弹丸能获得较高初速；对同一弹速有较好的可重复性^[2-5]。

为保证实验安全，降低气炮运行成本，常采用以压缩氮气为驱动源的二级轻气炮技术^[6]，发射过程为：通过高压气室内填充的高压气体膨胀做功，推动活塞加速运动，压缩前端轻质气体；当活塞临近锥段，它前端的轻质气体压力激增，超过活塞尾部的压力，使活塞减速；当膜片处的压力达到破裂临界值时，膜片破裂，受压的轻质气体推动弹丸在发射管内加速运动，达到预设速度^[7]。该驱动方式安全环保、不会产生有毒有害气体^[8-9]，但是压缩气体的能量密度远低于固体发射药的，如何有效提高能量利用率，获得较高弹速，是亟待解决的难题^[10-11]。

研究者们运用数值计算方法研究二级轻气炮发射的物理过程，已取得了一定的研究进展。林俊德^[12]在Seigel^[13]工作基础上，推导出一套针对高压气体驱动二级轻气炮的计算方法，以一维非定常描述气流运动，忽略活塞与管壁的摩擦，讨论了高压气室初始压力、膜片破膜压力及活塞质量对发射效果的影响，较客观地描述了发射参数和发射效果的关系。Groth等^[14]采用非定常准一维可压缩流动模型来预测泵管中轻质气体和弹丸头部低压气体运动，建立摩擦模型来描述活塞和弹丸的运动。目前气炮计算工作大多采用简化一维模型描述流场，该处理方式能较快得到计算结果，帮助归纳发射规律，但无法考虑活塞入锥过程的流场变化；锥段流场直接影响活塞减速过程及作用于弹丸尾部的初始压力，将该流场简化为一维模型会带来计算误差。

Piacesi等^[15]将气体、活塞、弹丸运动简化为沿孔径方向距离和时间的函数（即一维分析），但考虑了在变截面处气室与泵管、泵管与发射管的横截面积变化。胡天翔等^[16]在FLUENT中建立气炮二维数值模型，详细讨论了不同初始条件对气炮内弹道性能的影响。二维轴对称模型无法考虑气炮的非对称结构，如：高压气室与泵管连接处、锥段；高压气室与泵管处流场直接影响活塞加速过程，锥段流场直接影响活塞入锥及弹尾压力，造成计算误差。

随着计算机性能及数值方法的发展，采用三维有限元方法计算轻气炮发射过程成为可能。传统有限元方法计算气炮问题存在一些不足，拉格朗日法不能给出准确流场变化，欧拉法虽能描述气体流动，但对活塞、弹丸边界描述模糊，计算精度不够。而Noh^[17]提出的耦合欧拉-拉格朗日（coupled Eulerian-Lagrangian, CEL）算法，结合拉格朗日法和欧拉法的优点，不仅对材料运动分析较准确，而且允许材料发生大变形。欧拉域可以填充多种材料，对每个时间增量步计算气炮内高压气体、轻质气体体积分数，来确定气体在欧拉单元的流动轨迹。将CEL算法应用到轻气炮计算中，不仅能较精确地描述活塞、弹丸运动，还可以给出各个阶段气炮内气体分布、压力变化，为观测气炮发射过程、流场变化提供有效手段。

本文中，基于CEL算法，以高压气体驱动二级轻气炮为例，提出一种适用于多级轻气炮的分级数值模拟方法。由于CEL算法稳定时间步正相关于模型最小单元尺寸，膜片与气炮整体模型跨尺度，按照膜片几何尺寸划分全模型网格，单元数量将达千万量级，严重影响计算效率。因此，合理简化气炮模型，根据膜片破裂与否，将模型解耦为一级、二级模型，减小计算规模。目前，二级轻气炮内弹道计算没有活塞在高压段挤进变形过程的模拟，没有考虑摩擦、膜片破膜压力^[18]，导致计算弹丸终速与实验结果存在一定偏差。本文中设计正交试验，拟合确定计算重要参数——材料摩擦因数和膜片破膜压力。依据上述参数和分级模型建立数字化二级轻气炮，模拟结果与实验高度吻合。需要指出，该计算方法也适用于其他驱动形式多级轻气炮数值计算，对设计和优化多级轻气炮具有指导意义。

1.   数值模型

以14 mm口径高压气体驱动二级轻气炮为例，详细讨论气炮模型简化、分级依据及模型边界条件，提出依据CEL算法的分级数值模型。

1.1   计算方法的选择

14 mm口径高压气体驱动二级轻气炮模型如图1所示，由泵管堵头、活塞、高压气室、泵管、锥段、膜片、弹丸、发射管和靶室组成。将流场简化为一维模型，无法考虑活塞运动过程中高压气室与泵管连接处流场变化及锥段内流场变化，带来计算误差。简化二维轴对称气炮模型如图2所示，忽略非对称结构（高压气室与泵管连接结构和气孔），并假设高压气体直接作用于活塞尾部。真实发射过程是高压气体通过气室与泵管连接处气孔作用于活塞内凹处，推动活塞加速运动。

图  1  14 mm口径高压气体驱动二级轻气炮示意图

Figure  1.  Schematics of a 14-mm-caliber two-stage light gas gun based on high-pressure gas driving

下载: 全尺寸图片 幻灯片

图  2  14 mm口径高压气体驱动二级轻气炮二维模型

Figure  2.  Two-dimensional models for the 14-mm-caliber two-stage light gas gun

下载: 全尺寸图片 幻灯片

因此，为了拓展气炮数值计算方法，得到更精确的结果，在ABAQUS中采用CEL算法计算气炮三维模型。需要指出，二级轻气炮模型较大，三维方法无法像一维、二维方法较快得到计算结果，但有如下优势：(1)三维数值模型能考虑炮体非对称结构，用以研究活塞在锥段挤进变形过程、气固两相相互作用等问题；(2) CEL算法不仅较精确描述活塞、弹丸运动，而且考虑各阶段流场变化、气体分布，得到的流场更符合气炮真实发射状态；(3)欧拉场计算结果包括气体温度，可用于研究锥段内气炮烧蚀问题。

1.2   模型简化和分级

首先，简化计算模型的非必要部分：(1)省去阀门，忽略开阀过程。假设高压气体已填充，将活塞即将运动作为计算初始状态；(2)删去高压气室倒角。将高压气室简化为圆柱形，保持体积不变，便于划分高质量单元；(3)锥段内壁保持不变，其两端外径分别等同于泵管外径和发射管外径，并将锥段、泵管和发射管作为一个整体。

然后，需对简化后的计算模型划分单元。在CEL算法中，稳定时间步正相关于模型最小单元尺寸。单元过小不仅减小时间步，也会增加模型自由度，共同导致计算时间急剧增加。膜片与其他部件几何尺寸相差极大，膜片厚度仅约0.40 mm，但气炮总长可达19.36 m。要获得精确计算结果，模型中单元尺寸不可相差过大。若按照膜片几何尺寸划分单元并形成气炮全模型，单元数量将达千万量级，计算效率难以接受。

为此，从二级轻气炮工作原理出发，提出分级数值模拟方法，具体如下：(1)活塞和弹丸在膜片破裂前不相关，仅在膜片破裂后相互影响。(2)膜片破裂前，活塞在泵管内运动，挤压前端预先注入的轻质气体。气体无法透过膜片进入发射管，弹丸保持静止。此时膜片没有破裂，将其视为刚体，可得到膜片处气体压力历程。(3)测试膜片破裂压力发现，膜片破裂过程历时仅数微秒，远小于气炮发射过程历时。因此，可将膜片破裂视为瞬时完成，即膜片达到破膜压力时，膜片完全破坏。(4)膜片破裂后，受压的轻质气体进入发射管，推动弹丸加速运动。此时，活塞运动速度较低或接近停止，活塞尾部压力远低于锥段内气压，可忽略。

通过上述处理，将气炮解耦为两部分，如图3所示：一级模型为计算开始至膜片破裂，包含高压气室、活塞、泵管和锥段；二级模型计算膜片破裂至弹丸出膛，包含部分泵管、活塞、锥段、弹丸和发射管。

图  3  气炮模型

Figure  3.  Light-gas gun models

下载: 全尺寸图片 幻灯片

1.3   边界条件

对上述计算模型划分单元及设定边界条件。对炮体、活塞和弹丸，采用C3D8R型拉格朗日单元。对气体流经区域，即高压气室、泵管、锥段和发射管，建立欧拉域，采用EC3D8R型欧拉单元。经测试发现，当欧拉单元与拉格朗日单元尺寸比例保持1∶3到1∶2时，流体与固体接触效果较好。因此，气炮模型均保持该单元比例，模型最小单元尺寸为1.167 mm。

计算中每个单元都会记录各种参数，包括但不限于几何位置、速度、气体压力、密度和温度等。模型所有参数均存储在inp文件中，如：模型坐标存储在PART INSTANCE模块中，材料参数存储在MATERIALS模块中。可从该文件得到任意输出时间步气体参数、活塞和弹丸速度等。

以膜片在气炮轴向坐标x_d为分界线，从inp文件中分离出高压气室、活塞、泵管以及锥段模型参数生成一级模型，如图3(a)所示。因模型过于细长，图中隐去部分泵管。分别定义高压气体、低压气体在一级模型欧拉域中的初始区域。

对单元接触采用通用接触算法。该接触算法采用罚函数方法，自动指定主从面，充分考虑流体与固体、固体与固体的摩擦接触关系。由于本文中将气体处理为理想气体，气体与管壁碰撞完全弹性，没有动能损失，因此计算中仅考虑固体间的摩擦接触关系。对于一级模型，仅活塞与泵管产生相对位移，存在摩擦，即接触算法中摩擦因数定义为活塞与管壁的摩擦因数。需要注意，通用接触算法在同一计算中仅能定义一次且仅可定义一个摩擦因数。

完成一级模型计算后，根据膜片破膜压力，在时间-压力曲线上确定对应时刻，记为t_p。此时刻，活塞尾部对应气炮轴向坐标记为x_p，运动速度为v_p。根据x_p的数值，保留坐标大于x_p部分模型作为二级模型，包含部分泵管、锥段及发射管，如图3(b)所示。类似地，因模型过于细长，图中隐去部分发射管。将一级模型与二级模型重合部分的参数赋予二级模型，作为初始条件，包括气体密度、压力、温度以及活塞位置和速度。

二级模型活塞尾部因无气体存在，对活塞尾部施加稳态压力（压力数值为t_p时刻一级模型活塞尾部压力），近似模拟活塞尾部气体的作用。与一级模型类似，以弹丸和发射管的摩擦因数作为摩擦因数。需注意，此处存在2个误差：(1)活塞在泵管中向前运动，其尾部压力应动态变化，施加稳态压力将存在误差。一方面，活塞最大运动距离（即坐标x_p与锥段的距离）一般小于1 m，泵管全长11.82 m，活塞后部泵管内气体体积最大可膨胀8%，即活塞尾部压力变化应小于8%。另一方面，随着活塞入锥，其头部压力远超尾部压力，尾部压力误差影响更小可以忽略。(2)与一级模型不同的是，二级模型中活塞与泵管、弹丸与发射管均存在摩擦。因仅可定义一个摩擦因数，此处活塞与泵管的摩擦存在误差。考虑到二级模型中活塞运动距离较短，摩擦因数对活塞运动影响较小，且二级模型分析的主体为弹丸，活塞和发射管间的摩擦对活塞运动造成误差可以接受。

如此，二级轻气炮根据膜片破裂与否被解耦为一级模型和二级模型。类似地，该分级方法同样适用于其他多级轻气炮。

2.   参数设置

2.1   几何参数及工况参数

14 mm口径高压气体驱动二级轻气炮如图1所示，其几何参数见表1。其中泵管长度为气炮的完整泵管长度。经过分级处理后，一级模型泵管长度与此相同，为11.82 m。二级模型中，泵管保留部分是根据膜片破裂时活塞所处位置确定。因此，在不同工况下，二级模型保留的泵管长度一般不同。

表  1  14 mm口径高压气体驱动二级轻气炮几何参数

Table  1.  The geometrical parameters of the 14-mm-caliber two-stage light-gas gun

部件	长度/mm	内径/mm	入锥角度/(°)
高压气室	720	90
泵管	11820	60
发射管	6600	14
锥段	320	14～60	4.5

下载: 导出CSV 

| 显示表格

该气炮已有实验发射参数及弹丸终速见表2。3组实验高压气室初始压力、泵管初始压力及弹丸质量不同，活塞及膜片均相同，实验弹丸终速存在明显差异。高压气室初始压力是指高压气室内填充气体初压，3组实验均使用氮气。泵管初始压力是指泵管内气体初压，3组实验均使用氦气。实验2和3高压气室初始压力大致相同且明显高于实验1。弹丸质量为1.51、1.55和2.43 g，弹丸终速分别为3 846、4 155和3 571 m/s。在正交计算中，弹速区间为3.2～5.2 km/s。二级轻气炮弹丸发射速度范围为2～7 km/s，但常用速度范围仍是3～5 km/s，与本文中分析的弹速区间一致。

表  2  14 mm高压气体驱动二级轻气炮的实验参数

Table  2.  Experimental parameters of the 14-mm-caliber two-stage light-gas gun

实验	高压气室初始压力/MPa	泵管初始压力/MPa	活塞质量/g	弹丸质量/g	膜片厚度/mm	弹丸终速/（m·s−1）
1	15.7	0.041	733	1.51	0.4	3 846
2	22.1	0.045	733	1.55	0.4	4 155
3	21.7	0.047	733	2.43	0.4	3 571

下载: 导出CSV 

| 显示表格

活塞模型及其尺寸如图4所示，由聚碳酸酯和金属铝组合制作完成。其质量主要集中在金属制成的尾部，形变主要发生在聚碳酸酯制成的头部。该活塞在入锥时，能保持较快速度，更好地压缩前端气体。活塞头部与锥段角度相近，可以尽量减少活塞变形及摩擦做功消耗能量^[19]。

图  4  活塞模型

Figure  4.  The piston model

下载: 全尺寸图片 幻灯片

弹丸结构实际由弹托和其前端居中的球形弹丸组成，因不是本文的分析重点，仅将它们视作整体建模。弹丸模型及其尺寸如图5所示。弹丸材料为聚碳酸酯。实验1和实验2弹丸质量大致相同且明显低于实验3。

图  5  弹丸模型

Figure  5.  The projectile model

下载: 全尺寸图片 幻灯片

2.2   材料参数

对该二级轻气炮各部件赋值材料参数。炮体采用4340不锈钢，活塞由金属铝和聚碳酸酯组合制作而成。其中，4340不锈钢和铝使用Johnson-Cook本构模型^[20-21]:

式中：$\varepsilon$为等效塑性应变，${\dot{\varepsilon }}^{\mathrm{*}}=\dot{\varepsilon }/{\dot{\varepsilon }}_{0}$，$\dot{\varepsilon }$为等效塑性应变率，${\dot{\varepsilon }}_{0}$为参考应变率，${\dot{\varepsilon }}_{0}=1.0$ s^-1，A为材料在准静态下的屈服强度，B和n为应变硬化影响因子，C为应变率敏感指数，m为温度软化系数。4340不锈钢和铝的模型参数如表3所示。因本文中关注的重点在二级气炮分段计算方法，不详细讨论活塞入锥变形过程，仅考虑活塞入锥流场变化，所以聚碳酸酯使用线弹性材料模型，密度为1200 kg/m³，泊松比为0.3902，杨氏模量为2.5 GPa^[22]。

表  3  Johnson-Cook模型材料参数^[20-21]

Table  3.  Material parameters for the Johnson-Cook model^[20-21]

材料	密度/（kg·m−3）	泊松比	杨氏模量/GPa	初始屈服应力/MPa	硬化常数/MPa	应变硬化指数	温度系数	熔化温度/K
4340不锈钢	7900	0.33	200	280	802.5	0.622	1	1673
铝	2790	0.33	73.1	369	684	0.34	1.7	1000

下载: 导出CSV 

| 显示表格

实验中，高压气体氮气、轻质气体氦气均采用理想气体状态方程：

式中：p为气体压强，V为气体体积，T为温度，n为物质的量，R为摩尔气体常数。气体参数如表4所示。3组工况气体压力均不一样，根据文献[23-24]中氮气和氦气的参数，结合状态方程计算可得：实验1中氮气密度为180.86 kg/m³，氦气密度为0.069 kg/m³；实验2中氮气密度为244.03 kg/m³，氦气密度为0.075 kg/m³；实验3中氮气密度为240.34 kg/m³，氦气密度为0.079 kg/m³。

表  4  气体材料参数

Table  4.  Material parameters for gases

气体	R/（J·（mol·K）−1）	环境压力/Pa	比定容热容/（J·（kg·K）−1）
氦气	2077	0	3116
氮气	298.3	0	743

下载: 导出CSV 

| 显示表格

需要指出，在前人工作中较多忽略摩擦，导致计算弹丸终速与实验结果存在一定偏差^[25]。为优化气炮计算方法，得到较准确的弹丸终速，本文中将考虑材料摩擦。根据经验，铝与4340不锈钢摩擦因数在0.5～0.9之间，聚碳酸酯与4340不锈钢摩擦因数在0.2～0.4之间，具体数值与加工工艺相关。该精度摩擦因数不足以得到准确弹丸终速，结果误差较大。同时，由于采用分级计算方法，膜片破膜压力需要确定。根据膜片的材料和尺寸，预计膜片的破膜压力为15 MPa。而实验手段标定复杂，且难以测得动态加载下膜片的破膜压力，因此拟通过数值方法求解上述3个参数。

3.   关键参数的确定

根据前文分析可知，需确认活塞与泵管间的摩擦因数（一级摩擦因数μ₁）、弹丸与发射管间的摩擦因数（二级摩擦因数μ₂）和膜片的破膜压力（p），以获得准确弹丸终速。因此，以μ₁、μ₂和p为因素设计正交试验，用正交试验结果拟合参数方程。通过参数方程联立求解，得出3个参数具体数值。所以，将从正交试验设计方案、气炮参数拟合以及正交试验结果分析三方面，详细阐述气炮关键参数的确定方法。

3.1   正交试验

二级轻气炮发射过程受多因素耦合作用，线性参数方程不能较好拟合弹丸终速。因此，假设弹丸终速v与一级摩擦因数μ₁、二级摩擦因数μ₂和破膜压力p之间呈平方关系：

式(3)中共有10个参数，3个变量，需采用三因素四水平正交表确定所有参数的值。同时，考虑到水平梯度，一级摩擦因数μ₁取值为0.3、0.5、0.7和0.9；二级摩擦因数μ₂取值为0.1、0.2、0.3和0.4；膜片破膜压力p取值为10、20、30和40 MPa。设计L₁₆(4³)正交表如表5所示。

表  5  正交试验的因素及水平

Table  5.  Factors and levels of orthogonal tests

水平	μ1	μ2	p/MPa
1	0.3	0.1	10
2	0.5	0.2	20
3	0.7	0.3	30
4	0.9	0.4	40

下载: 导出CSV 

| 显示表格

为防止混淆实验参数与正交试验设计工况参数，需明确几个定义：实验1～3分别表示二级轻气炮的3个实验工况；3个实验工况对应3组正交试验，将这3组正交试验定义为试验1～3；每组正交试验，有16个计算工况；确定各因素数值后的验证工况，分别定义为实验1验证工况、实验2验证工况和实验3验证工况。

根据表5设计参数组合及计算工况，如表6所示。可以看出，随着一级、二级摩擦因数和破膜压力的变化，弹丸终速v变化明显。这进一步佐证确定摩擦因数和破膜压力的必要性。

表  6  正交工况及弹丸终速

Table  6.  The orthogonal cases as well as the final velocities of projectiles

计算工况	μ1	μ2	p/MPa	v/（m·s−1）				计算工况	μ1	μ2	p/MPa	v/（m·s−1）
				试验1	试验2	试验3						试验1	试验2	试验3
1	0.3	0.1	10	4517	4835	4183		9	0.7	0.1	30	3645	4428	4597
2	0.3	0.2	20	4751	4959	4243		10	0.7	0.2	40	4230	4279	4622
3	0.3	0.3	30	4797	5194	4142		11	0.7	0.3	10	3568	4036	3396
4	0.3	0.4	40	4782	5039	4180		12	0.7	0.4	20	3469	4019	3264
5	0.5	0.1	20	4453	4396	4184		13	0.9	0.1	40	4581	4755	3842
6	0.5	0.2	10	4081	4290	3611		14	0.9	0.2	30	4412	4550	4030
7	0.5	0.3	40	3997	4375	3887		15	0.9	0.3	20	4232	4451	3786
8	0.5	0.4	30	4197	4029	3720		16	0.9	0.4	10	3840	3843	3742

下载: 导出CSV 

| 显示表格

3.2   参数拟合

根据表6，采用最小二乘法，编写程序，拟合求解式(3)中参数a₁～a₁₀。对应3组实验工况得到3个拟合方程：

上述3个参数方程的相关系数分别为r₁=0.91, r₂=0.98, r₃=0.87，说明方程拟合效果好。同时，拟合参数方程的目的是计算μ₁、μ₂和p的具体数值。因此，将表2中3个实验工况的弹丸终速v₁=3 46 m/s，v₂=4 155 m/s和v₃=3 571 m/s代入式(4)～(6)，联立求解，得到满足μ₁∈[0.3,0.9], μ₂∈[0.1,0.4], p∈[10,40] MPa的结果：μ₁=0.82, μ₂=0.30, p=11.73 MPa。

进一步，需验证3个参数取值的准确性。将这3个参数代入实验1～3验证工况，计算得到实验1～3验证工况下的弹丸终速分别为4 088、4 367和3 485 m/s。结合表2可以看出，在二级轻气炮发射区间3.2～5.2 km/s内，每组验证工况的弹丸终速与实验的弹丸终速吻合度高。同时也说明，将弹丸速度与摩擦因数和破膜压力的关系描述为平方关系的假设是可行的。

3.3   结果分析

通过正交试验及参数拟合，确定了14 mm口径高压气体驱动二级轻气炮的一级摩擦因数、二级摩擦因数和破膜压力。而正交试验本身可应用于分析各因素的主次关系，分析结果如表7所示。

表  7  正交试验结果分析

Table  7.  Analysis of orthogonal test results

v/（m·s−1）	试验 1				试验 2				试验 3
	μ1	μ2	p		μ1	μ2	p		μ1	μ2	p
k1	4711.75	4299.00	4001.50		5006.75	4603.50	4251.00		4187.00	4201.50	3688.50
k2	4182.00	4368.50	4226.25		4272.50	4519.50	4456.25		3850.50	4126.50	3869.25
k3	3728.00	4148.50	4262.75		4190.50	4514.00	4550.25		3969.75	3802.75	4122.25
k4	4266.25	4072.00	4397.50		4399.75	4232.50	4384.00		3805.50	3682.00	4132.75
R	983.75	296.50	396.00		816.25	371.00	299.25		381.50	519.50	444.25

下载: 导出CSV 

| 显示表格

表7中：k_l (l=1,2,3,4)为使用对应因素水平l的结果均值；极差R为k_l (l=1,2,3,4)中最大值与最小值的差值，该差值的大小反映各因素对结果的影响程度。由表7可知，3个因素的极差远大于零，说明这3个因素对弹丸终速的影响不可忽略。一级摩擦因数将影响活塞运动，进而影响锥段内气体的压缩过程，使膜片破裂后作用于弹丸尾部的压力发生变化，最终影响弹丸速度。破膜压力阈值则会直接影响作用于弹丸尾部的初始压力，导致弹丸速度受到影响。二级摩擦因数直接影响弹丸的加速过程，导致弹丸速度受到影响。

工况1~2中，R_μ1>R_p且R_μ1>R_μ2，即一级摩擦因数为主要影响因素；而工况3中，R_μ2>R_p>R_μ1，即二级摩擦因数为主要影响因素。观察实验参数可以发现，实验3的弹丸质量明显大于实验1～2的，这导致在工况3中二级摩擦因数为主要影响因素。

正交试验结果表明，在轻气炮数值计算中，摩擦和破膜压力不能忽略。减小发射过程中摩擦消耗和增大膜片破膜压力有助于获得更高的弹丸终速。因此，提高加工工艺，保持炮体清洁，优化发射参数组合，能在一定程度上提高弹丸速度。而通过实验优化发射参数，将会带来大量人力物力消耗、气炮损耗等问题。采用本文中提出的计算方法，可以在得到气炮几何参数后，通过数值手段优化发射参数，以较低的成本获得更高的弹丸终速。

4.   数值化气炮分析

上一节，确定了活塞与泵管间的摩擦因数、弹丸与发射管间的摩擦因数和膜片破膜压力。将这3个参数代入实验验证工况，计算结果与实验结果拟合度高，说明该计算方法较好地复现气炮发射过程。本节将以实验1验证工况结果为例，通过数值化气炮分析，给出活塞弹丸速度曲线、压力曲线、关键时刻气炮内压力云图，并着重分析活塞入锥过程。

4.1   一级模型

一级模型如图3(a)所示，实验1发射参数为：气瓶气压，15.70 MPa；泵管气压，0.041 MPa；活塞质量，733.00 g；弹丸质量，1.51 g。实验弹丸终速为3 846 m/s，实验验证工况弹丸终速为4088 m/s。

针对高压气体驱动二级轻气炮发射初期，图6给出了0～7.50 ms内活塞前后的压力云图。图例显示压力范围为0～16 MPa，压力云图时间间隔为1.50 ms。该阶段高压气体快速作用于活塞尾部，活塞速度迅速上升。0～3.00 ms高压气室内维持较高初始气压，随着活塞运动，活塞尾部受到的压力持续下降。4.50 ms后，活塞尾部压力趋于稳定，活塞在泵管内加速运动。

图  6  0～7.50 ms活塞前后压力云图

Figure  6.  Pressure nephograms around the piston at 0−7.50 ms

下载: 全尺寸图片 幻灯片

活塞入锥阶段35.30～37.70 ms压力云图如图7所示。由于活塞运动前后压力差极大，采用统一图例不能较好反映压力变化过程，因此分3组图例。如图7(a)所示，35.30～35.90 ms图例显示压力范围为0～12 MPa，36.20～36.80 ms图例显示压力范围为0～35 MPa，37.10～37.70 ms图例显示压力范围为0～150 MPa，压力云图时间间隔均为0.30 ms。

图  7  泵管和锥段35.30～37.70 ms压力云图

Figure  7.  Pressure nephograms in the pump tube and the tapered section at 35.30−37.70 ms

下载: 全尺寸图片 幻灯片

图7(b)给出35.60～35.90 ms活塞头部及尾部局部压力，由图可知：前端反射压缩波在活塞头部与管壁缝隙处汇聚，其压力远高于活塞尾部压力，是导致该阶段活塞减速运动的主要原因。图7(c)给出相同时间内（35.60～35.90 ms）锥段局部压力云图，根据正交试验得到实验1膜片破膜压力为11.73 MPa，对应图7(c)中35.90 ms时刻，显然这时活塞还未入锥。图7(d)给出活塞入锥后（对应于37.40～37.70 ms时间内）锥段的局部压力云图，随着活塞运动，锥段内压力显著陡增。当活塞停止运动时，整体压力高于220 MPa。需指出的是，一旦膜片破裂，后续情形将不再真实。对应于本工况，即35.90 ms后的一级模型计算结果，如图7(d)所示，仅有参考意义，锥段处后续的真实压力云图及活塞运动将依赖于二级模型的计算。

活塞速度和位移随时间的变化如图8所示。该曲线可分为两部分：第1部分为活塞真实运动状态，对应的速度、位移曲线分别采用蓝色、红色实线表示；第2部分为一级模型边界条件下（膜片未破裂），活塞的后续速度、位移曲线分别采用蓝色、红色虚线表示。活塞从开始运动到入锥减速为0，共历时37.70 ms，运动位移11.89 m。0～33.30 ms活塞处于加速运动阶段；33.30～35.90 ms活塞速度小幅度波动，呈减速运动趋势；35.90～37.70 ms活塞处于减速运动阶段，历时极短，速度陡降为0。

图  8  活塞的速度/位移-时间曲线

Figure  8.  Velocity- and displacement-time curves of the piston

下载: 全尺寸图片 幻灯片

膜片处的压力变化如图9所示。该曲线同样分为2个部分：第1部分为活塞真实运动状态，对应的压力曲线采用蓝色实线表示；第2部分为一级模型边界条件下（膜片未破裂），活塞的后续压力变化曲线，采用蓝色虚线表示。可以观测到，在10.00 ms前，膜片处压力变化为零，活塞前端压缩波还未传到。随后压力出现小幅振荡，产生一系列压缩波，在膜片、锥段和活塞间来回反射。30.00 ms后，压力明显升高。当压力升至11.73 MPa时，膜片破裂，对应时刻t_p=35.90 ms。

图  9  膜片处的压力-时间曲线

Figure  9.  Pressure-time curve at the diaphragm

下载: 全尺寸图片 幻灯片

同时结合活塞速度曲线发现，活塞还未入锥，就已开始减速运动。破膜时刻35.90 ms，活塞速度为402.23 m/s，此时活塞接近锥段，但还未入锥。在该工况气炮实验中，也发现相同结果，活塞还未入锥，膜片已经破裂。说明该气炮发射能力还有提升可能，通过更换其他材质或厚度的膜片，提高膜片破膜压力，可使弹丸获得更高的弹尾压力。

如前所述，一级模型将膜片处理为刚体，计算中不考虑膜片破裂过程。根据膜片破膜压力，在一级模型中确立膜片破裂时刻。实验1膜片破膜压力为11.73 MPa，对应破膜时刻35.90 ms。因此，35.90 ms之后一级模型计算结果，并非实验1工况下活塞真实运动过程，活塞真实运动将在二级模型中反映。当更换其他膜片时，只需根据破膜压力，确定破膜时刻，将该时刻下活塞运动状态、气体参数赋值于二级模型，即可完成二级模型计算。

4.2   二级模型

二级模型初始状态为膜片完全破裂，弹丸开始加速运动。由于分析对象是弹丸，因此将初始时刻重新定义为0 ms，对应一级模型膜片破裂时刻35.90 ms。

二级模型0～1.95 ms压力云图如图10所示，红圈标识为弹丸位置。其中0～0.30 ms图例显示压力范围为0～12 MPa，压力云图时刻间隔为0.1 ms；0.40～0.70 ms图例显示压力范围为0～20 MPa，压力云图时间间隔为0.1 ms；0.80～1.40 ms图例显示压力范围为0～80 MPa，压力云图时间间隔为0.2 ms；1.60～1.95 ms图例显示压力范围为0～150 MPa，压力云图时间间隔为0.1 ms，最后一张图截取活塞速度降为0时的压力云图，间隔为0.15 ms。当压力大于图例最大值时仍显示为红色。

图  10  二级模型0～1.95 ms的压力云图

Figure  10.  Pressure nephograms of the second-stage model at 0−1.95 ms

下载: 全尺寸图片 幻灯片

从图10可以观测到，膜片破裂后，锥段泄压，弹尾压力迅速上升，在发射管内加速运动。泄压初期0～0.30 ms，弹丸运动不明显，压缩波在活塞与弹丸间来回反射，锥段压力急剧上升。同时，受锥段（变截面）和弹丸形状（内凹形）的影响，在压力云图（0.10、0.60 ms时尤为明显）中可观测到波阵面凹凸变化，向中轴线汇聚。0.80～1.95 ms，伴随活塞入锥，弹丸速度快速上升，活塞速度急速下降，锥段内压力激增，在1.70 ms（37.60 ms）达到最大值，为150.30 MPa。活塞在1.95 ms停止运动，此时锥段内压力降至约60 MPa。

由4.1节可知，需通过二级模型反映膜片破裂后活塞的真实运动情况。因此，将对比分析膜片破裂后一级模型和二级模型的活塞运动状态，其时间-速度/位移曲线如图11所示。图中蓝色实线为二级模型的活塞速度，蓝色虚线为一级模型的活塞速度，红色实线为二级模型的活塞位移，红色虚线为一级模型的活塞位移。可以看出，二级模型的活塞在1.95 ms（37.85 ms）停止运动，减速耗时比一级模型的活塞（37.70 ms）多0.15 ms，且二级模型活塞位移11.94 m比一级模型活塞位移11.89 m多0.05 m，活塞入锥更充分。

图  11  活塞的速度和位移时间里程曲线

Figure  11.  Velocity- and displacement-time curves of the piston

下载: 全尺寸图片 幻灯片

导致上述变化的原因如下：一级模型将膜片处理为刚体，未考虑膜片破裂，导致在锥段形成密闭气室；随着活塞入锥，内部压力激增，整体压力达到220 MPa，远高于二级模型锥段内的最高压力，导致活塞速度快速下降。因此，一级模型描述的活塞减速过程并非真实状态，通过二级模型可观测到活塞真实减速过程。需注意，本文计算未考虑活塞塑性变形和热损耗，活塞入锥时未发生挤压变形。活塞速度降为0后，锥段内压力远高于活塞尾部压力，活塞有反向运动趋势，但不明显。

活塞停止后，弹丸仍在发射管内加速运动，压缩波在活塞与弹尾间来回反射。以2.20～3.20 ms弹丸尾部压力云图为例，分析活塞入锥后弹丸运动状况和弹尾压力变化历程。2.20～3.20 ms以弹丸为中心建立随体坐标系，压力云图如图12所示。其中2.20～2.35 ms图例显示压力范围为10～30 MPa，压力云图时刻间隔为0.05 ms；2.40～2.55 ms图例显示压力范围为10～20 MPa，压力云图时间间隔为0.05 ms；2.60～3.20 ms图例显示压力范围为0～12 MPa，压力云图时间间隔为0.20 ms。当压力超过图例显示上限时仍显示为红色。可观测到：随着弹丸加速运动，发射管内波阵面存在明显分层。这表明，发射管内压缩波来回反射，导致弹丸尾部压力呈阶跃式下降变化。弹丸离开发射管时，弹丸尾部压力仅约6 MPa。

图  12  弹丸尾部2.20～3.20 ms的压力云图

Figure  12.  Pressure nephograms at the tail of the projectile at 2.20－3.20 ms

下载: 全尺寸图片 幻灯片

弹丸尾部压力和速度随时间的变化曲线如图13所示。从弹尾压力-时间曲线可观测到，尾部压力出现5次较明显的波峰，对应时刻分别为0.50、0.85、1.30、1.55和2.25 ms。波峰压力分别为43.19、47.67 、71.09、82.65和48.19 MPa。其中弹丸尾部压力在1.55 ms即第4次波峰达到最大值82.65 MPa，随后急剧下降。从弹丸速度-时间曲线可观测到，弹丸运动可分为3个阶段：第1阶段弹丸加速初期，时间为0～1.50 ms；第2阶段弹丸加速中段，时间为1.50～2.40 ms，弹速迅速提高，明显快于第1阶段；第3阶段弹丸加速末段，时间为2.40～3.20 ms，弹速缓慢提高。最后，弹丸在3.20 ms（39.10 ms）时离开发射管，速度为4 088.29 m/s。

图  13  弹丸时间-速度/弹尾压力曲线

Figure  13.  Pressure-time curve at the tail of the projectile and velocity-time curve of the projectile

下载: 全尺寸图片 幻灯片

由上述分析可知，压缩波在活塞与弹丸尾部之间来回反射，导致发射管内压力变化极大，流场分布复杂，难以分析压缩波的叠加过程。在发射管设置观测点，记录压力峰值到达的时刻和压力变化历史，结合弹丸运动情况，可对发射管内的流场作进一步分析。在发射管内沿模型中轴线近似等分均布11个观测点，观测点分别与膜片的轴向距离S如表8所示，膜片处为观测点0。监测各观测点的压力到达时刻及变化历史，同时记录弹丸的运动过程。将每个观测点首次出现压力阶跃式变化的时间t记录于表8。图14为各观测点和弹丸弹尾压力变化历史。根据相邻观测点间距ΔS及压力阶跃式变化时间间隔Δt，可求得气体流场到达各观测点时的波阵面平均速度c。

表  8  观测点的参数及波阵面平均速度

Table  8.  Parameters for observation points and mean wavefront velocities at observation points

观测点	S/m	ΔS/m	t/ms	Δt/ms	c/(m·s−1)		观测点	S/m	ΔS/m	t/ms	Δt/ms	c/(m·s−1)
0	0		0				6	3.260	0.650	2.35	0.20	3250
1	0.018	0.018	0.05	0.05	360		7	3.900	0.640	2.50	0.15	4267
2	0.670	0.652	1.40	1.35	483		8	4.550	0.650	2.65	0.15	4333
3	1.310	0.640	1.75	0.35	1829		9	5.200	0.650	2.80	0.15	4333
4	1.960	0.650	1.95	0.20	3250		10	5.850	0.650	3.00	0.20	3250
5	2.610	0.650	2.15	0.20	3250		11	6.490	0.640	3.15	0.15	4267

下载: 导出CSV 

| 显示表格

图  14  等分均布观测点和弹丸弹尾的压力变化对比

Figure  14.  Comparison of pressure changes between the evenly-distributed observation points and the projectile tail

下载: 全尺寸图片 幻灯片

结合图10和表8综合分析可知：0～0.018 ms泄压初期，锥段内压缩波迅速进入发射管，作用于弹丸尾部，波阵面平均速度为360 m/s。由于观测点0与观测点1时间间隔极短，计算平均速度可反映出锥段内波阵面速度，由图8可知破膜时活塞速度为410.80 m/s，与该值接近。0.018～1.95 ms波阵面平均速度及弹丸速度变化明显，该阶段压缩波在发射管内来回反射，导致管内及弹尾压力激增，从图14也可观测到，观测点1及弹丸尾部压力波动式陡增。1.95～2.35 ms波阵面平均速度均为3 250 m/s，对应弹丸速度由2 833 m/s升高到3 648 m/s。需注意，1.95 ms后活塞停止运动，观测点1处压力降至约32 MPa，明显低于发射管内压力。2.50～3.15 ms波阵面平均速度最高达4 333 m/s，临近出膛时为4 267 m/s，与弹丸速度4 080 m/s接近。

各观测点波阵面压力变化最终趋于一致，弹丸离开发射管时，管内压力在4.00～7.00 MPa之间，略低于弹丸尾部压力。这是因为，观测点压力及弹丸尾部压力均为静压，而压缩波传到弹丸尾部时会发生反射，部分动压转为静压，导致在弹尾测得的静压高于观测点处的静压。弹丸在加速初期（0～1.50 ms）和加速末期（2.40～3.20 ms）的速度变化明显慢于加速中期（1.50～2.40 ms）的，延长加速中期时长，使弹丸有效加速时间增长，能更高效地加速弹丸。

5.   讨　论

基于CEL算法，以高压气体驱动二级轻气炮为例，提出一种适用于多级轻气炮的分级数值模拟方法。计算中合理简化、解耦气炮模型，设计正交试验，确定气炮重要参数数值，以此建立数字化二级轻气炮，获得与实验高度吻合的模拟结果。选取实验1验证工况计算结果，详细描述二级轻气炮发射过程，给出泵管、锥段和发射管的压力云图，直观呈现活塞和弹丸在气炮内的运动过程以及活塞和弹丸前后流场的变化，着重分析活塞入锥过程压力变化。

实验1验证工况二级轻气炮发射总时长为39.10 ms，膜片于35.90 ms时破裂，弹丸开始运动，其运动时长为3.20 ms。弹丸运动历时短，仅占发射总时长的8.18%。整个发射过程，弹丸尾部压力显著变化，过载随时间分布很不均匀。实验1验证工况弹丸终速（4 088 m/s）对于相应的实验弹丸终速（3 846 m/s），计算误差为6.29%，导致该计算误差的原因可能是：(1)二级模型省去部分泵管，将活塞尾部动态压力加载简化为稳态压力加载，此时泵管省去部分气体压力变化虽小，但仍与赋予活塞尾部的稳态压力有一定差异；(2)根据弹丸与发射管间的摩擦因数定义二级模型摩擦，导致活塞与泵管间的摩擦因数小于实际值；(3)计算工况中假定泵管、发射管内完全清洁，但这在实验中难以保证。

图7反映出活塞减速过程分为2个阶段：第1阶段是由膜片处反射的压缩波在头锥处汇聚引起，此时减速不明显；第2阶段是因为锥段内压力上升，远高于活塞尾部压力，导致活塞速度迅速下降。在第2阶段，锥段压力进一步激增，更加高效加速弹丸。通过图7和图9可以发现，活塞还未入锥膜片就已破裂，当更换破膜压力更高的膜片时，可以使弹丸获得更高弹尾压力。因此，选取合适的发射参数，可以降低气炮损耗，同时最大限度地发挥气炮的发射潜能。仅通过发射实验来优化，会带来巨大人力物力消耗。而采用本文中提出的分级计算方法，建立数字化高压气体驱动二级轻气炮，来判断实验发射参数的可行性。结合数值计算和实验， 可高效优化发射参数，大量减少实验次数。

膜片破裂时刻，活塞动能为60 949.26 J，弹丸离开发射管时，动能为12 619.16 J，能量利用率为20.70%。实验2～3中，膜片破裂时刻，活塞动能分别为74 744.95和76 409.20 J，弹丸动能分别为14 779.17和14 758.99 J，能量利用率分别为19.77%和19.32%。可以看出，3组实验活塞动能只有约20%转化为弹丸动能，其他转化为压缩气体的内能以及活塞与管壁、弹丸与管壁摩擦做功耗散等。

综上所述，采用本文中提出的分级计算方法和气炮关键参数确认方法，可建立数字化高压气体驱动二级轻气炮。通过数值方法，认识气炮发射的物理过程，优化气炮发射参数，达到提高弹丸动能、降低膛内壁温度和优化弹丸尾部压力的目的。

需要指出，本文的重点是论述气炮解耦计算方法及气炮关键参数确认方法。采用该方法，更换活塞材料模型，可用于研究活塞入锥变形过程。而活塞入锥过程流场变化极其复杂，采用解耦计算思路，构建部分泵管及锥段模型，加密输出时间步及减小欧拉域单元尺寸，可详细分析活塞入锥过程流场变化。由于锥段类似于变截面激波管过渡段，该计算方法也可用于激波管流场研究。欧拉场的计算结果包括气体温度，可得到气炮内的温度场。以实验1验证工况为例，锥段内气体最高温度为3 956.8 K，如图15所示；36.90～37.70 ms锥段的瞬态温度场如图16所示。同时，本文中采用CEL算法，计算中考虑了气体与管壁间的相互作用；当更换活塞和气炮管壁的材料模型，考虑活塞运动过程材料熔化、变形，结合气炮瞬态温度场，可以得到管壁各阶段的准确温度，用以研究气炮的烧蚀问题。

图  15  0～39.10 ms锥段气体温度曲线

Figure  15.  Temperature-time curve of the gas in the tapered section at 0−39.10 ms

下载: 全尺寸图片 幻灯片

图  16  36.90～37.70 ms锥段瞬态温度场

Figure  16.  Transient temperature fields in the tapered section at 36.90−37.70 ms

下载: 全尺寸图片 幻灯片

6.   结论与展望

采用CEL算法，对14 mm口径高压气体驱动二级轻气炮展开三维计算。计算中，以膜片破裂与否为分级依据，将气炮完整模型解耦为一级模型和二级模型。一级模型包括高压气室、活塞、泵管和锥段；二级模型包括部分泵管、锥段和发射管。设计正交试验，确定活塞与泵管间的摩擦因数为0.82，弹丸与发射管间的摩擦因数为0.30，膜片破膜压力为11.73 MPa。将这3个参数代入实验验证工况，计算得到的弹丸终速与实验结果吻合度高。

同时，通过呈现实验1验证工况的计算结果，展现了该方法的优势：在充分考虑摩擦、流场的条件下，开展14 mm口径高压气体驱动二级轻气炮发射过程的数值计算，可以得到活塞和弹丸的运动速度曲线、气体压力云图以及活塞入锥的完整过程。

值得一提，本文中的弹速区间为3～5 km/s。对于3 km/s以下工况，计算方法完全适用；对于更高弹速工况，该方法具有可推广性。需要注意，更高弹速意味更高的气室初压，将导致整个发射过程气体压力峰值升高，可能会影响计算稳定性。调小时间缩放因子或调整网格尺寸比例，可避免算例崩溃。在该方法框架下，更换活塞和气炮管壁的材料模型，结合气炮各阶段瞬态温度场及气体温度，可在后续工作中探究管壁烧蚀等问题，更好地优化气炮的发射性能。对于固体发射药驱动、爆轰驱动等其他驱动形式的二级/多级轻气炮，本文中提出的气炮简化方法、分级思想和气炮关键参数确定方法同样适用。

有以下几个创新点：

(1)整合建立数字化高压气体驱动二级轻气炮，化繁为简将其解耦为一级、二级模型。同时，对气炮发射参数的优化和气炮发射过程中流场的变化展开了详尽分析。

(2)采用CEL算法，从三维层面对二级轻气炮计算，更真实地模拟气炮的发射过程。

(3)结合正交试验，推导确定了二级轻气炮的关键参数，包括活塞与泵管摩擦因数、弹丸与发射管摩擦因数和膜片破膜压力。