基于五阶WENO格式的燃气在药床中流动过程二维两相流研究*

程申申，王 浩，薛 绍，陶如意

（南京理工大学能源与动力工程学院，江苏 南京 210094）

膛内点传火过程是一个复杂的高温高压传热传质过程，包括点火管火药燃烧、点火管燃气流出以及膛内药床的传热燃烧过程，对火炮的药床燃烧以及后续推进过程有重要的影响。在实际实验中，现有测试手段很难得到实际状态下这一阶段的实验现象。因此，开展药床燃气流动的内弹道模拟对整个火炮发射过程的研究具有重要意义[1-3]。

近年来针对点火管和膛内火药床燃烧的研究成果有很多。王浩等[4]通过点传火特性实验，对比了不同装填条件下的传火管多点压力曲线及火焰传播速率，分析了底火焰在传火管中的传输特性，是早期对点传火进行的实验研究。王珊珊等[5-6]对点火管和中心管建立了一维两相流模型进行内弹道两相流研究，分析了不同结构点火管的点传火性能。郭梦婷等[7]对中心炸管定容阶段进行了二维两相流研究，探究了中心管式抛撒机构的内弹道流场特性。但是以上研究都着眼于点火管内火药燃烧和传播，对于燃气在膛内流动及火药燃烧的研究较少。Miura 等[8]使用带破裂膜片的燃烧装置进行实验，并采用两相流体动力学代码进行模拟计算，用于探索膛内负压差产生的机理,并得出结论,压差是由固体推进剂运动导致的，压差波动取决于底火喷口的流量和喷口的位置，但实验和模拟过程未对火焰发展及流场规律进行详细描述和观测。Jaramaz 等[9]对整个火炮膛内流动过程进行实验和一维两相流研究，但是未考虑径向效应。程诚等对某制导炮弹二维两相流内弹道性能进行研究，并采用高阶黎曼近似模型求解两相流方程，详细研究了带推进系统的整个内弹道过程中膛内流场特性及其发展[10-11]，但主要研究火药床燃烧及弹丸推进过程，并未仔细对初始阶段点传火过程进行流场分析。因此，开展中心传火管装药点火过程的理论与试验研究，对进一步分析火炮发射起始阶段膛内的点火过程燃气流动规律，降低起始压力波的产生，保证火炮射击的安全性至关重要。

本文中，拟设计低压可视化实验装置，观察点火药燃气在药床中的流动过程。考虑到可视化装置的承压能力，采用假药粒装填结构，并结合二维两相流数值模拟来预测这一阶段的实际过程。使用五阶精度的WENO格式，相比于经典的Mac Cormack 格式与TVD格式，其具有精度高、间断分辨清晰，而且效率高的优点，是十分有效的算法[12-14]。以期研究结果对进一步开展火炮起始阶段燃气流动及压力变化规律研究具有应用价值。

1 数理模型

1.1 物理模型

对膛内初始阶段点传火过程建立的计算模型如图1所示。膛内填充由尼龙材料仿真实药粒大小制成的假药粒，中心为点传火管，内装小粒黑火药，左端用端盖密封，装有底火以点燃点火管内火药，右端用铝合金膜片密封。

图1 装药结构示意图Fig.1 Schematic of the charge structure

根据上述模型建立真实药粒运动燃烧和假药粒运动的两相流模型，为了简化内弹道数值模拟的复杂性，需对实际物理过程进行适当假设，基本假设如下：

（1）采用双流体模型假设，将发射药颗粒群作为具有连续介质特性的拟流体来处理;

（2）忽略湍流影响和气固两相间黏性；

（3）火药颗粒服从几何燃烧和指数燃烧定律[2]；

（4）气相状态方程采用Nobel-Abel状态方程[2]；

（5）假设底火内点火药完全以燃气形式加入点火管内，点火管采用一维两相流计算，破膜后源项在膛内当地释放；

（6）假药粒不产生燃气，尺寸与真实药粒相同。

1.2 膛内轴对称二维两相流方程组

对于内弹道膛内气-固两相流过程，由于膛内沿身管轴线方向呈对称形式，因此在内弹道膛内气-固两相流二维计算时，往往采用轴对称模型，守恒形式可以表示如下：

式中：U为时间守恒矢量，F、G分别为径向和轴向的矢通量，H为轴对称变换矢量，S为源项。具体可表示如下：

提供辅助方程来封闭基本方程。膛内点传火过程源项由点火管提供，点火管为一维线源，从破膜后点火燃气流入膛内过程视为源项当地释放。源项与膛内燃烧过程进行耦合计算，具体一维点火管控制方程及流出流量公式可参考文献[5]，其他如状态方程、燃烧规律、相间阻力、相间传热、颗粒间应力、颗粒表面温度等辅助方程等参见文献[2]。

2 数值计算方法

基于有限差分法，对上述模型在时间项上采用三阶TVD型龙格库塔法进行时间推进，空间项采用五阶WENO格式进行空间离散。

2.1 五阶WENO格式

WENO格式基本思路是将基架点分为多个组，每个组独立计算某一点的导数逼近，得到多个差分，由于流场中间断不可能处处存在，所以根据每个模板的光滑程度，设定权重，如果光滑度越高，则权重越大，存在间断则权重趋近零。

对于向量形式的守恒方程：

采用全局Lax-Friedrichs分裂对空间项F、G进行流通矢量分裂：

利用分裂后的矢通量构造矢通量空间导数项的差分逼近式，以F项为例：

采用Jiang 等[14]的五阶精度WENO格式构造矢通量导数逼近式的通量形式为（以正特征值为例）：

hj+1/2为正通量，由3个通量加权平均而得：

加权系数为（k=1,2,3）：

式中：C1=1/10，C2=6/10，C3=3/10，ε0=10−6，β=2，Ik为光滑度量因子。Ik的表达式为：

当特征值为负时，利用对称性，将上式中所有下标“j+m”换成“j−m”即可。对G项的离散同理。

2.2 时间方向离散

空间导数离散化后的半离散方程为：

式中：Q n为离散后空间导数项，R n为源项。使用三阶TVD龙格库塔方法进行离散：

∆t为时间步长，可以根据下式确定：

式中：C为柯朗数，dx为计算网格尺度；λ*为全局最大特征值，由上述算法可计算得出。

2.3 初始条件和边界条件

膛内装填条件采用如图1所示的装填结构，并结合实验装填数据和环境参数，确定模拟过程中的初始条件。膛底和出口处边界条件采用镜面反射方法，计算至膜片破裂。由于WENO格式采用多个基架点进行计算，在边界外设置3个虚拟点可将边界点转化为内点参与计算，还可以保持边界点的计算精度。中心轴线边界满足轴对称关系式。

3 实验验证与数值结果分析

3.1 假药粒实验

为了验证所建立模型和数值方法的可靠性，对假药粒情况进行实验验证。采用了针对105 mm 口径火炮药室缩比的可视化实验装置，如图2所示。装置一端端盖内侧连接点火管，外侧可以装底火，用以激发点火点燃点火管内的点火药；另一端端盖为中通结构，用以锁紧闭压膜片。装置中间装有耐高压玻璃管，管内装有定量的假药粒，通过装置中部预留的矩形窗口可以有效观测装置内部的实验现象。装置上方设置3 个传感器接头，分别对应药室的膛底、中部和弹底（膜片）处，可以通过压力传感器测得这几处的压力。

整体实验装置如图2所示，可视化窗口为400 mm×120 mm，高速摄像系统可通过此窗口拍摄实验过程中火焰发展传播过程。传感器距膛底分别为40、247.5和455 mm，可通过数据采集系统获取各位置压力。整体实验装填点火药为2#小粒黑火药，质量为30 g；假药粒为PA 尼龙，质量为3 kg；膜片为铝合金，厚度为0.5 mm。

图2 可视化实验装置Fig.2 Visual experimental system

3.2 数值计算与实验数据比对

以上述实验装置进行可视化模拟实验的内弹道过程二维两相流动为研究对象，其具体结构可参见图1中的简化物理模型，计算参量设置为：点火药质量，30.0 g；点火药密度，1650 kg/m3；点火药燃速指数，0.98；点火药燃速系数，3.68 mm·s−1·MPa−0.98；火药力，300 J/g；点火药平均直径，1 mm；假药粒密度，1100 kg/m3；假药粒质量，3.0 kg；药室长度，495 mm；药室直径，110 mm；点火管长度，25 mm；点火管直径，19 mm；气体绝热指数，1.3；网格长度，5 mm；气体余容，1 cm3/g。

为了简化计算，忽略点火管的二维效应，即在传火孔破膜后将点火管作为线性点火源加入膛内进行两相耦合计算，弹底与膛底为垂直于轴向的固壁平面，采用反射壁面条件。通过压力传感器对膛内各位置采集的压力与数值模拟获取压力进行比对，由于模拟计算只到破膜后截止，因此只取破膜前的压力进行对比分析。

图3所示为中间测点压力随时间变化的实验结果与模拟结果的对比曲线。可以看出，模拟结果与实验结果符合较好，全局平均误差为5.35%，最高峰值处（即2.5 ms后区域）误差仅为1.42%，说明所建立的模型和采用的数值方法合理有效。

图3 计算压力曲线与实验对比图Fig.3 Comparison of pressure-time curves between numerical simulation and experiment

整个膛内点传火过程从点火管火药点燃开始，高温高压燃气经由点火管小孔流入膛内，点燃主装药。燃气在点火管燃烧破膜以及在膛内流动需要一定的时间，而压力传感器测点在管壁处，所以初始阶段压力基本没有变化，从点火开始到1 ms时间段内，均保持在一个大气压左右；在随后的1～2 ms时间段内，压力快速上升，这是由于假药粒装填几乎占满整个药室，药室透气性很差，药室自由空间很小，膛内压力上升也较快；压力会在到达峰值时，会产生波动现象，这种波动现象主要由于假药粒装填密度过高，装药床透气性差，且假药粒不会燃烧导致燃气流动性变差，局部压力会产生振荡。

3.3 膛内两相流动过程分析

通过高速摄像拍摄可视化实验装置内火药燃气流动过程，与计算得到的压力场进行对比，可以清晰反映火药燃气在膛内的发展状况。由于可视化装置只能观察到膛内部分区域火焰传播过程，而火焰反映了高温燃气的流动过程，进而反映膛内的压力场的变化过程，因此对计算得到的压力云图进行相对应位置截取对比，如图4所示。

图4 数值模拟压力云图与高速摄影图像的对比Fig.4 Comparison between pressure clouds by numerical simulation and images by high-speed photography

整个点传火可分为3 个阶段；初始破孔阶段，在压力云图中可以看出随着点火管依次破孔，小孔破孔位置出现不同程度压力梯度，而在高速摄影中，由于点火管被假药粒包围，很难直接观测到初始阶段的破孔现象；随着点火管高温燃气不断向膛内喷射，高温燃气在膛内向膛壁和药室两端传播，在压力云图中可以观察到其压力梯度沿径向和轴向都有明显的变化，相对应的高速摄影图像中，火焰传播也呈现与压力云图同样的特征；到达2.5 ms时，点火管对应的膛内区域燃气流动在径向上基本一致，即沿径向的二维效应基本消失，而在远离点火管的区域仍存在一定的径向效应，从点火管开孔处向膛底和弹底形成明显的轴向压力梯度。

为有效呈现膛内压力随时间的发展变化，在径向上取r=20 mm 处观察压力沿轴向的变化，如图5所示。图5(a)展示点火破膜初期阶段（0.5～1.2 ms）膛内压力沿轴向的变化，这一阶段膛内径向效应明显，压力从破孔区域不断上升，且随着点火管依次破孔，压力出现了相对应的峰值，并逐渐向周围发展。随着点火管燃气不断向膛内喷射，点火管对应的膛内区域压力基本趋于稳定，如图5(b)所示：2.5 ms后中心位置压力约18 MPa，径向压力梯度基本消失，此时在轴向压力梯度作用下，燃气向药室两端传播，膛底和弹底的压力不断升高，沿轴向压力变化区域也不断增大。

图5 r=20 mm 处沿轴向方向的压力分布Fig.5 Pressure distributions along the axis at r=20 mm

图6和图7描述了r=20 mm 处膛内气相的径向速度和轴向速度沿药室轴向的变化规律。从图6可以看出，点火管破孔初始阶段，燃气不断从点火管内流向膛内，此时点火管附近区域燃气的径向速度要远大于其轴向速度，沿径向传播也比较快。初始破孔阶段，由于膛内压力远低于管内压力，气体喷射出来速度极大（由于压力差导致由点火管流入膛内气相速度很大），随着燃气不断流入，径向压力差变小，径向的燃气速度也逐渐降低。

图6 0.5～1.0 ms 不同时刻r=20 mm 处沿轴向方向的气相速度分布Fig.6 Gas phase velocity distributions along theaxis at r=20 mm and different times from 0.5 msto 1.0 ms

图7 1.5～3.2 ms不同时刻r=20 mm 处气相速度沿轴向的分布Fig.7 Gas phase velocity distributions along theaxis at r=20 mm and different times from 1.5 msto 3.2 ms

在点火1.5 ms后，大量燃气由点火管破孔处向两端流动，并逐渐堆积形成高压区，在1.5～2.0 ms内，燃气逐渐向膛底和弹底运动，气相轴向速度渐渐升高，如图7(a)所示；待燃气运动到膛底或弹底时遇到阻力后又气相轴向速度慢慢降低至反向运动，如图7(b)所示，由于高密度装填，在弹底形成了药粒壅塞，所以气相速度在接近弹底处已开始有反向的波动；随着燃气轴向运动加剧，燃气的径向速度在点火管对应区域基本为零，在2.0 ms前燃气向膛底和弹底方向传播，由于燃气压力较低，初始传播时受颗粒扰动较大，径向速度波动剧烈，如图7(c)所示；2.2 ms后，由于燃气充盈膛底和中部区域，燃气径向速度除在膛底有小范围波动外，其他区域燃气径向速度降为零，燃气的径向扰动逐渐向弹底传播，如图7(d)所示。

火药颗粒运动对初始阶段膛内点传火研究很重要，为进一步研究药床中药粒运动规律，图8～9展示了r=20 mm 膛内固相速度沿药室轴向发展变化。从图8可以看出，破孔初始阶段火药颗粒轴向和径向速度非常低，最高速度不超过0.5 m/s，说明膛内的火药颗粒基本保持不动，呈现一种多孔介质状态。从实验实际拍摄的高速影像来看，固相颗粒确实未发生大范围运动。

图8 0.5～1.2 ms不同时刻r=20 mm 处固相速度沿轴向的分布Fig.8 Solid phase velocity distributionsalong the axis at r=20 mm and different times from 0.5 ms to 1.2 ms

随着火药燃气在药床内不断流动，势必会带动火药颗粒运动。在高密度装填的药床中，火药颗粒可运动的范围十分有限，且模拟和实验采用的假药粒不会燃烧且变形较小，所以颗粒间相互作用力十分大，导致药床内的固相速度非常低。如图9(a)所示，火药颗粒径向速度随时间发展逐渐为零，从破孔位置沿轴向至弹底方向有一定的径向速度，至膛底处径向速度也逐渐趋向于零，与径向气相速度相对应，说明气相速度也是影响颗粒运动的重要原因。图9(b)中火药颗粒轴向速度也验证了这一点，由于燃气轴向压力梯度大，气相轴向速度高，所以导致火药颗粒向膛底和弹底两侧运动速度较高。在真实装药情况下，会造成膛底火药粒挤压破碎，形成小范围高压区域，继而推动药床向弹底运动。

图9 2.5～3.2 ms不同时刻r=20 mm 处固相速度沿轴向的分布Fig.9 Solid phase velocity distributionsalong the axis at r=20 mm and different times from 2.5 ms to 3.2 ms

4 结 论

以火炮膛内初期点传火性能研究为背景，建立了轴对称二维两相流数理模型，对假药粒进行相应合理化假设，在空间上采用五阶WENO格式进行空间离散，时间上采用三阶TVD型龙格库塔法进行时间推进，并进行实验与数值模拟对比研究。数值分析了点火初期膛内压力的发展及气固两相速度，具体结论如下。

（1）可视化实验装置与高精度数值方法，有效模拟了初始阶段点火燃气在药床中的传播规律，通过与可视化实验平台测试数据和现象进行对比，验证了数理模型和数值方法的有效性。

（2）点火初期膛内径向效应明显，随着传火不断进行，点火管对应膛内区域径向效应逐渐消失，燃气的轴向传播开始占据主导地位，膛内压力呈现径向一致、轴向梯度分布的特征。

（3）通过数值模拟与实验观测可以得出，在高密度装填情况下，初始阶段药床药粒基本不会发生运动，最高速度为0.5 m/s，后续固相速度的发展与气相速度有紧密的关系；气相速度在初始阶段沿径向传播较快，后续在轴向压力梯度作用下，轴向速度开始占据主导，在膛底和中部区域气相速度降低为零，径向速度扰动也逐渐向弹底传播，在到达弹底之前由于药粒壅塞现象导致阻力变大，气相速度会出现振荡。