基于气固两相反应流的温压炸药能量释放规律数值模拟及实验验证

李 根，卢芳云，李翔宇，李志斌

(国防科技大学 文理学院，湖南 长沙 410073)

引 言

温压炸药(TBX)以其在密闭空间的高效毁伤能力广泛应用于打击藏匿于山洞、地道中的目标。温压炸药一般由高能炸药、高热值金属颗粒、氧化剂、黏结剂构成，爆轰产物包含大量可燃组分，能够与环境中的氧气持续发生反应，在密闭空间形成高温、高压、缺氧环境，对目标造成“软杀伤”效果。

根据温压炸药的组成，其能量释放过程包括3个阶段[1]：首先，高能炸药组分快速分解发生无氧爆轰反应，此过程通常在几微秒内完成，金属颗粒在此阶段几乎不参与反应；接着，爆轰产物中的氧化性气体与金属颗粒进行无氧燃烧反应，此阶段不需要外界氧气参与，反应持续数百微秒；最后，前两个阶段生成的可燃组分与周围空气发生有氧燃烧反应，此阶段可持续数毫秒。

国内外学者对温压炸药的爆炸问题进行了大量实验研究，发现其能量释放过程受到气体环境[2-4]、颗粒含量[5]、空间体积[6]等因素的影响。实验虽然可以通过测量压强、温度等流场参数侧面反映温压炸药在不同条件下的能量释放差异，但是对温压炸药能量释放机理研究不足。数值模拟是另一种研究温压炸药能量释放过程的重要方法。国内通常采用JWL-Miller模型对温压炸药爆炸问题进行数值模拟[7-9]，通过调整模型参数能够得到与实验比较相符的计算结果，但是该方法并不能真实地反映流场中气体和颗粒的温度、速度的差异以及由此产生的对燃烧速率的影响。国外以Khasainov、Kuhl等[10-11]为代表建立了温压炸药爆炸流场的两相流模型，能够比较详细地给出流场压强、温度、产物浓度等信息，但是由于对固相计算采用欧拉方法，不能跟踪颗粒运动以及与冲击波的相互作用过程。

本研究基于气固两相反应流模型，对温压炸药能量释放规律开展数值研究，建立了有限差分-物质点耦合算法，其中气相和固相分别采用能够捕捉波系结构和追踪颗粒运动的有限差分法和物质点法，通过与实验结果对比验证算法的可靠性和准确性。在此基础上，研究气体环境、铝粉含量、空间体积对温压炸药能量释放的影响规律，以期加深对温压炸药密闭空间爆炸及后燃效应的认识。

1 温压炸药爆炸流场两相反应流模型

1.1 两相流控制方程

温压炸药爆炸后流场物质包括气体爆轰产物、固体颗粒、环境气体。根据两相流理论，将气体爆轰产物和环境气体作为气相，固体颗粒作为固相，两相具有不同的速度和温度场，并存在质量、动量和能量交换。根据连续介质力学原理以及对离散颗粒介质的“拟流体”假设，可建立气固两相的控制方程[10-11]:

气相质量方程：

(1)

气相动量方程：

(2)

气相能量方程：

(3)

气相组分质量方程：

(4)

固相质量方程：

(5)

固相动量方程：

(6)

固相能量方程：

(7)

(8)

对于铝颗粒，K=150s/cm2。由于液态铝的热膨胀系数大于固态氧化铝，在温度升高过程中液态铝的膨胀使固态氧化膜破裂，液态铝与空气提前接触，因此铝颗粒的点火温度介于铝熔点和氧化铝熔点之间，取932K[12]。

此外，气固相互作用还包括作用力fs和传热速率qs，具体表达形式可参考文献[11]。

1.2 化学反应模型

1.2.1 爆轰反应

在无氧爆轰阶段，高能炸药组分快速分解释放能量，反应在几微秒之内完成，相比于无氧和有氧燃烧反应时间很短，因此在数值计算时可采用瞬时爆轰模型。无氧爆轰能量释放可通过爆炸反应方程式和盖斯定律确定，本研究采用Brinkly-Wilsion(B-W)原则确定炸药爆轰产物组成[13]。对于RDX，其爆炸反应方程式为：

(9)

无氧爆轰阶段生成的气体产物和释放能量作为初始条件用于之后流场演化的计算。

1.2.2 燃烧反应

CaHbAlcOd=
Cn1+n2+n3H2(n4+n5)Aln6+2n7On2+2n3+n5+3n7+2n8

(10)

原则上，燃烧产物的组成可根据当前单元的压强和温度,依据最小自由能原理[14]通过反复迭代计算确定。然而在流场计算时单元数量众多，为提高计算效率，采用含金属炸药产物组分的经验确定方法[13]，即首先将金属元素完全氧化，剩余元素再按B-W原则处理，因此燃烧产物生成顺序为Al2O3、H2O、CO、CO2。于是，根据质量守恒可以得到这些组分在燃烧反应后的分子数密度m1，m2，m3…m8。在已知反应前后分子数密度及组分生成焓ΔHi的前提下，根据盖斯定律可以得到该单元发生燃烧反应释放的热量ΔQb：

ΔQb=-∑(mi-ni)ΔHi

(11)

(12)

式中：Mi为组分相对分子质量。

1.3 状态方程

由于流体运动及燃烧反应的存在，流场各处组分不同且时刻发生变化。单一的高能炸药爆轰产物的JWL状态方程并不能反映组分改变带来的影响。为此，Donahue等[15]在分析了JWL状态方程各项的物理意义后提出一种描述燃烧产物的JWL状态方程：

可以看出，与爆轰产物JWL方程相比，燃烧产物JWL方程低压项参数ω由组分的质量分数Fi、定压比热容Cp,i、定容比热容Cv,i共同决定，这反映了低压条件下组改变的影响。除此之外，中高压项参数A、B、R1、R2与爆轰产物JWL方程参数保持一致。RDX爆轰—燃烧产物高压项参数为：A=573.2GPa，B=14.6GPa，R1=4.6，R2=1.4，爆热Ed=5.18MJ。

2 数值模拟方法及实验验证

2.1 有限差分-物质点耦合算法

温压炸药爆炸后流场中存在稀疏波、冲击波、接触间断等复杂波系结构，颗粒在其中运动的同时发生后燃反应释放能量，这些能量反过来影响波的传播。因此，流场波系结构的捕捉以及颗粒运动的跟踪对于分析认识流场演化规律及后燃反应过程非常重要。本文采用有限差分-物质点耦合算法(FDM-MPM)对两相流控制方程(1)求解，其中气相有限差分计算采用具有良好间断捕捉特性的Wave Propagation格式[16]，固相采用能够跟踪质点运动且不受网格畸变影响的物质点法[17]，气固相互作用直接在网格节点计算，最终形成计算程序。下面分别进行介绍。

2.1.1 有限差分计算

气相控制方程的有限差分计算采用均匀的笛卡尔网格，不考虑源项作用时其一阶离散格式为：

(14)

式(14)中：Qg为气相守恒变量；Fg、Gg、Hg为气相数值通量。右端通量差分计算采用wave propagation方法。该方法属于Godunov格式的一种，通过求解单元界面的局部黎曼问题对通量差分进行分裂，x方向的计算格式为：

(15)

式(15)中：λx为黎曼分解波波速；Wx为波间断；“+”号代表右行波，“-”代表左行波。同理，y和z方向的通量差分分裂计算分别为：

(16)

(17)

高阶空间离散可通过MUSCL、WENO等方法对界面变量进行重构实现，时间离散可采用具有TVD性质的龙格-库塔方法。

2.1.2 物质点计算

(18)

其中：η代表相对质量(以当前时刻为参考)；Cs为颗粒材料比热容；Ts为固相温度。

物质点法采用欧拉-拉格朗日双重描述，首先将固相介质划分为有限个拉格朗日质点，质点变量包括质量mp，位置xp，速度vp，温度Tsp。在计算时采用欧拉背景网格，通过插值函数NI,J,K将质点变量映射到背景网格节点，建立固相常微分方程如下：

(19)

(20)

(21)

式中：VI,J,K为网格节点体积，节点固相变量为：

(22)

对常微分方程(19)～(21)求解采用显式计算，计算格式为：

(23)

(24)

(25)

质点变量的更新通过将更新的节点变量映射回质点来实现：

质量更新计算：

(26)

位置更新计算：

(27)

速度更新计算：

(28)

温度更新计算：

(29)

至此，固相物质点的质量、位置、速度、温度计算全部给出。

2.1.3 算法耦合及计算流程

考虑到气固两相占据相同区域，且有限差分法和物质点法计算都需要使用欧拉网格，因此可以采用一种新的网格系统，如图1所示。该网格在计算过程中作为有限差分离散网格，始终储存气相信息。同时，又作为物质点法的背景网格，临时储存固相信息，在更新变量映射回质点变量之后清零。

图1 有限差分-物质点耦合算法网格Fig.1 Grid system of FDM-MPM

(30)

计算气固两相反应流的有限差分-物质点耦合算法在Visual Studio 2019平台实现，采用C/C++语言编写，算法计算流程如图2所示。

图2 有限差分-物质点耦合算法计算流程图Fig.2 Computational flow chart of FDM-MPM

2.2 实验验证

温压炸药装药构型包括活性外壳型和复合相爆炸混合型[18-19]。活性外壳型采用中心高能炸药和外围活性颗粒的装药形式，而复合相爆炸混合型采用高能炸药和活性颗粒混合混合形式。由于单质炸药爆轰过程不受颗粒的影响，因此活性外壳构型中的颗粒质量分数可以达到60%以上，从而有效提升温压炸药的温压效应。本小结采用的验证实验和后续温压炸药能量释放规律研究讨论的都是活性外壳装药构型。

对Kuhl等[11]进行的温压炸药密闭容器爆炸试验(见图3)进行数值模拟，以验证有限差分-物质点耦合算法的有效性和准确性。

图3 温压炸药装药密闭容器爆炸试验图Fig.3 Explosion test of a thermobaric explosive charge in a closed container

试验容器为圆柱体，直径20cm，长21cm，监测点位于容器顶部距离中心5cm(r=5cm，z=10.5cm)处。温压炸药由高能炸药太安(PETN)和铝(Al)粉装药，为圆柱体，位于容器中心位置。其中球形PENT位于装药中心，质量0.5g，密度1.0g/cm3，外围Al粉质量1.0g，填充密度0.6g/cm3。考虑到容器和装药形状满足轴对称条件且爆炸位置关于z=0平面对称，为提高计算效率，采用1/4模型，如图4所示。网格尺寸为0.0625cm。模型中心PETN和空气采用网格离散，外围Al粉采用物质点离散。PETN采用瞬时爆轰模型，初始爆轰产物采用JWL状态方程，高压项状态方程参数和气体组分的物性参数与原文献[11,20]一致。

图4 PETN/Al爆炸问题数值离散模型Fig.4 Numerical model of PETN/Al explosion

监测点处超压的实验与数值模拟结果对比如图5(a)所示。可以看出，两者的冲击波到达时间和一次峰值比较一致，0.4ms内压力曲线吻合良好。数值模拟结果能够反映出冲击波早期的多峰现象，当压力变化平稳时数值模拟结果与实验结果较为符合。

图5 监测点处超压及冲量的实验和数值模拟对比Fig.5 Comparison of experimental and numerical simulation of overpressure and impulse history at the gauge point

监测点冲量反映了压力变化的平均特征，实验和数值模拟的冲量对比如图5(b)所示。可以看出，冲量曲线近似为一次函数，这说明容器内压力经过一段时间后趋于稳定，其斜率可代表准静态压力。数值模拟得到的冲量曲线与实验结果基本一致，验证了其对温压炸药在密闭空间爆炸后空间准静态压力预测的准确性。

需要说明的是，对于活性外壳装药构型的温压炸药，中心高能炸药爆炸为理想爆轰，用于爆炸流场计算的物性参数与装药尺寸无关，因此可以适用于大尺寸装药的温压炸药。这为后续研究温压炸药能量输出规律奠定了基础。

3 温压炸药能量释放规律

3.1 计算条件

数值计算采用轴对称模型，如图6所示。以目前国内外常用的RDX/Al型温压炸药为对象，装药总质量100g，采用分离构型，中心为RDX，密度为1.6g/cm3；外围为Al粉，填充密度0.6g/cm3。爆炸容器为长径比1∶1的圆柱体，侧壁中间位置布有监测点。RDX爆轰产物状态方程中高压项参数见上述1.3节，铝粉物性参数为：颗粒密度ρm=2.7g/cm3，点火温度Tign=935K，颗粒直径ds=4.3μm，燃烧热Q=30MJ/kg。

图6 RDX/Al型温压炸药爆炸计算模型Fig.6 Computational model of RDX explosion

3.2 气体环境对能量释放的影响

设爆炸容器尺寸为Ф80cm×80cm，对铝粉质量分数40%的温压炸药在空气和氮气环境中的爆炸问题进行数值模拟，气体环境参数见表1。

表1 气体环境参数

监测点处超压和冲量时程曲线如图7所示，其中冲量曲线可以用线性函数来拟合，斜率代表容器内准静态压力。

图7 不同气体环境监测点压力冲量曲线Fig.7 Pressure and impulse history at the gauge in different gas environment

由图7可以看出，相比于氮气，温压炸药在空气中的爆炸冲击波峰值压力和空间准静态压力都大于氮气，且准静态压力增幅(增幅80%)高于峰值压力(增幅20%)。

图8为气体环境对温压炸药铝粉燃烧的影响。可以看出，在空气中铝粉会与其中的氧气发生燃烧反应，在温压炸药爆炸后3ms内反应基本结束，而在氮气中铝粉质量几乎不发生变化。

图8 不同气体环境容器内铝粉剩余质量Fig.8 Residual mass of Al powder in a closed container in different gas environment

图9展示了气体环境对温压炸药后燃能量释放的影响。可以看出在空气中温压炸药后燃能量释放非常明显，而在氮气中几乎不发生后燃反应。

图9 不同气体环境温压炸药后燃能量释放Fig.9 Afterburning energy of TBX in different gas environment

为进一步研究气体环境对温压炸药爆炸特性的影响，表2给出了不同空间体积的峰值压力和准静态压力。可以发现，峰值压力增幅20%以上，且随着空间体积的增大而增大。准静态压力增幅都在80%，呈现先增大后减小趋势。

表2 不同尺寸容器的峰值压力和准静态压力增幅

从表2还可以看出，准静态压力增幅明显高于峰值压力。这是因为铝粉后燃反应速率较慢，对峰值压力的提升体现为减缓冲击波在传播过程中的衰减，且空间体积越大(冲击波传播距离越远)，提升越明显，但增幅作用有限；而空间准静态压力与温压炸药释放总能量相关，与能量释放速率无关，因此增幅更加显著。以上结果表明，在密闭空间中温压炸药后燃效应主要表现为空间准静态压力的提升。

3.3 铝粉质量分数对能量释放的影响

对不同铝粉质量分数(20%、40%、60%)的温压炸药在不同体积容器(V1=Ф80cm×80cm、V2=Ф120cm×120cm)内的爆炸问题进行数值模拟，得到监测点处压力和冲量曲线如图10所示，峰值压力和准静态压力见表3。可以发现，峰值压力随铝粉含量的变化规律与空间体积相关：在V1容器中，由于容器体积较小、壁面监测点距离装药较近，高能炸药爆轰反应释放的能量决定冲击波强度，峰值压力随铝粉含量增加而减小；在V2容器中，容器体积较大、壁面监测点距离装药较远，此时铝粉燃烧释放能量减缓了冲击波在传播过程中的衰减，峰值压力随铝粉含量增加而增加。以上规律说明，铝粉含量影响温压炸药爆炸冲击波的衰减，铝粉含量越高，冲击波衰减得越慢。而空间准静态压力只与温压炸药能量释放总量有关，因此随铝粉含量的增加而增加。

图10 在不同容器中铝粉含量对监测点压力和冲量的影响Fig.10 Influence of Al content on pressure and impulse history at the gauge position in different vessels

表3 不同铝粉含量对应的峰值压力和准静态压力

3.4 空间体积对能量释放的影响

对铝粉质量分数40%的温压炸药在不同体积容器(V1=Ф80cm×80cm、V2=Ф100cm×100cm、V3=Ф120cm×120cm、V4=Ф320cm×320cm)的爆炸进行数值模拟，得到空间体积对温压炸药后燃能量释放的影响，如图11所示。

图11 温压炸药在不同体积容器后燃能量释放Fig.11 Afterburning energy of TBX in different specific volumes

由图11可以看出，温压炸药在V1、V2、V3容器内后燃反应最终释放的能量趋于一致，而在V4容器内释放的能量较少，这说明当容器体积增大到一定程度时温压炸药后燃反应释放的能量会降低。图12给出了温压炸药中铝粉在不同体积容器内燃烧时剩余质量随时间的变化情况。可以看出铝粉反应进程在初始阶段(< 0.193ms)相同，之后V1、V2、V3容器内的铝粉分别在0.193、0.343、0.505ms时刻加快了反应速率，这说明反应进程受空间体积影响。

图12 不同体积容器内剩余铝粉质量Fig.12 Residual mass of Al powder in a closed container with different specific volumes

图13给出了这些时刻不同体积容器内的压力和温度分布，其中黑色圆点代表固相物质点。由图3可以看出，在这些时刻不同体积容器内的冲击波在到达壁面后开始反射，与铝粉颗粒发生相互作用。壁面反射波加剧了铝粉与空气的混合，促进了铝粉的燃烧，因此在0.193、0.343、0.505ms时刻，V1、V2、V3容器内的铝粉反应速率发生明显变化。而V4容器由于体积过大，冲击波的壁面反射作用不能抵消爆轰产物在膨胀过程中温度降低带来的影响，最终导致铝粉反应不完全。

图13 不同体积容器内压力和温度分布Fig.13 Pressure distribution in a closed container with different specific volumes

图14 空间体积对温压炸药铝粉反应程度的影响Fig.14 Influence of specific volume on the reaction degree of Al powder of TBX

根据Feldgun理论模型[6]，V1、V2、V3、V4容器内的氧气可以使RDX爆轰产物和Al粉完全燃烧，而计算结果表明，V4容器由于体积过大导致铝粉没有完全反应，这说明铝粉反应度(铝粉反应度=1-剩余铝粉质量/初始铝粉质量)与空间体积相关。对铝粉质量分数20%、40%、60%的温压炸药在不同体积容器的爆炸进行数值模拟，得到铝粉反应度与空间体积的关系，结果如图14所示。为使结果具有普适性，横坐标采用比空间体积(比空间体积=空间体积/炸药质量)。可以看出，铝粉反应度随比空间体积增加呈现下降趋势，当比空间体积超过100m3/kg时，铝粉反应度低于90%左右。当比空间体积一定时，铝粉反应度随着铝粉含量的增加而降低。

4 结 论

(1)有限差分-物质点耦合算法能够准确反映温压炸药在密闭空间爆炸的冲击波多峰现象，压力和冲量曲线的计算结果与实验结果符合良好。

(2)温压炸药中铝粉与空气中氧气的燃烧反应提高能量输出。相比于氮气环境，在空气中冲击波峰值压力增幅在20%以上，随空间体积的增大而增大；空间准静态压力增幅在80%以上，随空间体积的增大呈现先增大后减小趋势。

(3)铝粉含量越高，温压炸药爆炸冲击波衰减得越慢，空间准静态压力越高。

(4)密闭空间中壁面反射波有利于加快铝粉燃烧速率，提高反应程度。当比空间体积超过100m3/kg时，铝粉反应度下降到90%以下，且铝粉含量越高，其反应程度越低。