铝粉尘爆轰的反应模型*

岳军政，洪 滔，吴先前，黄晨光

（1. 中国科学院力学研究所，北京 100190；2. 北京应用物理与计算数学研究所，北京 100094）

铝粉作为一种常见的工业原料，在其生产过程中容易与空气混合发生两相爆轰造成巨大危害[1]，同时因其含能较高，且价格便宜，也常作为燃料加入燃料空气炸弹(fuel air explosive, FAE)中形成大规模粉尘爆轰以制造破坏。针对铝粉的工业生产安全及其国防应用，对悬浮铝粉尘的两相爆轰已开展了大量研究[2-3]。由于实验观测难度大以及实验重复性较低，数值模拟成为研究粉尘爆轰的重要手段。

悬浮铝粉尘两相爆轰是非理想爆轰过程，常用经典的ZND (Zel’dovich-von Neumann-Döring)爆轰模型描述其爆轰波结构，其数值模拟的难点在于铝粉的燃烧模型。在目前的研究中，铝粉燃烧模型主要有扩散模型和动力学模型两种。在扩散燃烧模型中，化学动力学反应速率很快，燃烧速率主要受氧化气体的扩散速率控制[4]；而在动力学燃烧模型中，铝颗粒表面的氧化气体扩散速度较快，主要由化学动力学反应速率控制铝颗粒的燃烧速度[5-6]。Glorian 等[7]考虑颗粒大小以及燃烧环境，对表面反应动力学和气相机制进行了理论分析。在现有的数值模拟工作中，Veyssiere 等[2]提出了一个两步反应模型，通过拟合动力学参数得到了与实验结果一致的计算结果。Zhang 等[8]提出了同时考虑扩散机制和动力学机制的混合模型，并且讨论了两个燃烧机制的适用条件。此外，Briand 等[9]改进了两步反应模型，并且发展了一个混合模型。然而，在这些模型中有一个指前参数Zhyb需要拟合，而目前对这个参数并没有统一的认识[10]。

本文中，拟通过考虑铝粉燃烧产物氧化铝(Al2O3)在高温下的分解，针对铝粉的扩散燃烧模型，增加燃烧产物的逆反应吸热机制。此外，将改进的铝粉反应模型嵌入到并行的粉尘爆轰数值模拟程序中，分别对铝粉/空气混合物以及铝粉/氧气混合物的两相爆轰进行三维数值模拟，将计算得到的稳定爆轰波速度与实验结果和文献值进行对比，获得爆轰流场的物理量分布。

1 计算模型

1.1 模型假设

由于铝粉尘气-固两相爆轰是复杂的物理化学过程，因此对计算模型作以下假设[11]：（1）铝粉不可压缩，且均匀分布在气体中；（2）铝粉为直径相等的球形颗粒，并且片状颗粒等效为球形；（3）单个铝颗粒内部没有温度梯度和压力，重力、颗粒之间的作用以及热辐射作用可忽略；（4）铝粉燃烧所释放的能量均被气体吸收。

1.2 控制方程

数值模拟中采用两相流模型[12]描述悬浮铝粉尘爆轰的发展过程。模型中将气、固两相均看作连续介质，气体和固体满足各自的守恒方程，并且考虑气体和固体颗粒之间的质量、动量和能量的交换。描述铝粉尘爆轰的三维非定常两相流模型方程如下。

气相质量、动量和能量守恒方程分别为：

固相质量、动量和能量守恒方程分别为：

式中：下标1 和2 分别指气相和固相铝颗粒；ρ 为密度；φ为体积分数，且φ1+φ2=1；e为比内能；u和v分别为气相和铝颗粒的速度矢量；p为气体压力；I为燃烧引起的单位体积内铝粉的质量变化率；F为气体对固体颗粒的作用力矢量；Q为两相间的对流热传导速率；qAl为单位质量铝粉燃烧释放的热量。

颗粒数守恒方程为：

式中：n为单位体积内的铝颗粒数。

气体组分方程为：

式中：yj为气体组分j的质量分数，ωj为单位体积内气体组分j的质量变化率。

气体状态方程为：

式中：T为温度，R为普适气体常数，wj为气体组分j的摩尔质量。

两相比内能分别为:

式中：cV1为 气体比定容热容，随气体温度T变化[13]；cV2为铝的比热容，e2的方程中已考虑相变。

气体对铝颗粒的作用力F表达式为[12]：

式中：r为铝颗粒半径，Cd为拖拽系数。拖拽系数的表达式为：

式中：Re为雷诺数，且Re= 2ρ1|u−v|r/µ，µ为气体黏性系数。

两相间的对流热传导速率为[12]：

式中：λ1为气体导热系数；Nu为Nusselt 数，且Nu=2+0.459Re0.55Pr0.33；Pr为普朗特数，且Pr=µcp/λ1，cp为气体比定压热容。

式（13）中的因子β 是针对片状铝颗粒而言，计算中将片状铝颗粒等效为相同质量的球形颗粒，β 的表达式为：

式中：Sf为片状铝颗粒的实际表面积，Ss为等效球形颗粒的表面积，Tm为铝的熔点温度。之所以引入β 是由于片状颗粒的表面积比等效球形颗粒的表面积大，热传导时吸热较多；当铝熔化后，铝液滴呈球形，β=1。而对于球形铝颗粒，式（13）中不用考虑β，或者β=1。

1.3 铝粉反应模型

铝粉的反应模型是本文的创新点。两步反应模型应用较广泛，即将爆轰中铝粉的反应过程分为点火前的诱导阶段和点火后的燃烧阶段。针对诱导阶段，Levitas 等[14]考虑氧化铝壳的保护作用，通过理论分析认为铝颗粒达到其熔点温度后即可点火；洪滔等[15]通过分析铝颗粒的激波点火机制，也认为铝颗粒在爆轰动力学流场中的点火温度为其熔点温度。因此，本文中也采用该点火判据，当铝粉达到其熔点温度变为液态后开始燃烧反应。对于燃烧阶段，铝粉化学反应简化为 4 Al+3O2→2Al2O3，反应速率I采用扩散燃烧模型[16]：

式中：Lv为氧化铝的蒸发潜热。对于铝粉/空气混合物爆轰，假设气体中存在4 种组分，分别为O2、N2、Al(g)和Al2O3(l)；对于铝粉/氧气混合物爆轰，不存在N2组分。

计算模型的主要参数为：铝密度ρ2=2 700 kg/m3，铝熔化温度Tm=931 K，铝比热容cV2=905 J/(kg·K)，铝的燃烧热qAl=31.5 kJ/g，氧化铝蒸发潜热Lv=4.77 kJ/g。

上述控制方程(1)～(16)为含有源项的欧拉方程，为了较好地对爆轰流场中的强间断问题进行数值模拟，本文中采用三维的时-空守恒元解元(space-time conservation element and solution element method,CE/SE)差分格式求解不含源项的欧拉方程，然后用四阶Runge-Kutta 法求解其源项部分，开发了悬浮铝粉尘两相爆轰的三维数值模拟程序，并且基于消息传递接口(message passing interface, MPI)技术，实现了程序的并行化设计。

2 数值方法验证

为了验证开发的数值模拟程序的可靠性，分别采用1 和2 mm 两种网格尺寸对激波管问题进行了数值模拟，并与理论结果[20]进行了对比。在一根尺寸为1.0 m×0.1 m×0.1 m 的激波管中充满理想气体，距左端0.5 m 处有一隔膜将气体分为两部分，t=0 时刻两部分气体的状态为：

计算开始时刻隔膜突然打开，左右两部分气体开始混合。图1 分别显示了某一时刻数值计算得到的密度和压力分布与理论值的对比。由图1 可以看出，采用这两种网格计算得到的结果与理论结果均吻合较好，表明采用1～2 mm 的网格均可以达到精度要求。计算结果表明，采用本文中开发的程序可以较好地捕捉模拟激波问题，从而验证了该程序的可靠性。

图1 数值模拟得到的密度和压力分布与理论计算结果[20]的比较Fig. 1 Distributions of density and pressure by simulation and theory[20]

3 计算结果与分析

为了验证改进的铝粉反应模型，分别对铝粉/空气混合物爆轰以及铝粉/氧气混合物爆轰进行数值模拟。

3.1 铝粉/空气两相爆轰

首先对Tulis 等[21]开展的爆轰管内片状铝粉/空气混合物爆轰实验进行数值模拟。实验所用的爆轰管长度为5.487 m，内径为152 mm。铝粉的质量浓度σ2为0.3 kg/m3，电子显微镜显示该片状铝粉的比表面积为3～4 m2/g，本文计算中将该铝粉等效为直径为3.4 µm 的球形颗粒，相应地，β=5。

考虑到爆轰管的轴对称特性，建立1/4 计算模型。图2(a)显示了计算区域的横截面(yOz)，尺寸为76 mm×76 mm，在y和z两个方向上的网格数量均为61，即网格尺寸为1.246 mm。在笛卡尔坐标系中对爆轰管的圆周作近似处理，由图2 可以看出，爆轰管壁由水平和垂直的单元组成。计算域x方向长度为6 m，网格数量为4 800，即网格尺寸为1.25 mm。为了提高计算效率，管壁外的节点单元（见图2(a)中灰色部分）不参与计算。基于MPI，采用120 个进程并行计算，每个进程计算50 mm 长的计算域，并且记录爆轰管轴上每隔1 m 位置处的压力变化。

图2 计算域横截面和局部区域初始压力云图Fig. 2 Transverse cross section of the computational domain and initial pressure contour of partial domain

爆轰管的左端和圆周均采用固壁条件，右端6 m 处开口，并且在面xOy和xOz上的边界采用周期边界方法。图2(b)为初始时刻局部区域(x=0～300 mm)压力云图，其中左边31 mm 是点火段。本文参考文献[22]，点火条件模拟炸药起爆后的高温高速爆轰产物，属于强点火条件。由于实验中采用2.8 g RDX的点火源，并且模拟采用1/4 模型，所以计算中采用高温高压的0.7 g N2作为点火源，保证点火源质量和能量均为实验点火源的1/4。点火区域的初始参数为：φ1=1, ρ1=5 kg/m3,T1=1 500 K,ux=3 235 m/s,uy=uz=0。右边蓝色区域为静止的悬浮铝粉尘和空气两相混合区域，初始参数和实验值相同：ρ1=1.21 kg/m3, σ2=0.3 kg/m3,T1=T2=298 K, |u|=0,|v|=0。

图3 显示了数值计算得到的爆轰管轴上1～6 m 各位置处压力随时间的变化，铝粉在前导激波的作用下点燃并燃烧释放能量，支持爆轰波压力逐渐升高，并使其在5～6 m 处发展为稳定爆轰波，这和实验观测现象一致。

图3 爆轰管轴线不同位置处压力历史Fig. 3 Pressure histories at different locations on the tube axis

图4 显示了计算得到的t=3.83 ms 时刻爆轰波阵面附近气体压力、气体温度、气体密度、气体轴向速度、铝粉体积分数以及铝粉轴向速度的分布。由图4 可以看出，此时爆轰波阵面到达5.84 m位置，在前导激波的作用下，气体被迅速压缩至高温高压状态，气体和铝粉尘均被加速，铝粉尘的聚集导致局部浓度升高。当前导激波过后，铝粉并不会立刻燃烧释放能量，而是在其点燃前有个诱导点火区域，因而前导激波后气体的温度相对后面铝粉燃烧区域的气体温度较低，这一点也可以从气体温度云图中看出。由于管壁的摩擦作用，可以看到两相速度均存在径向的梯度分布。

图4 3.83 ms 时刻铝粉/空气爆轰波阵面附近气体压力、气体温度、气体密度、气体轴向速度、铝粉体积分数以及铝粉轴向速度的分布Fig. 4 Distributions of gas pressure, gas temperature, gas density, gas velocity in the x direction, volume fraction of Al particles, and Al particles velocity in the x direction around the detonation wave for Al/air mixtures at t=3.83 ms

数值模拟得到的稳定爆轰波速度为1 536 m/s，与实验结果1 510 m/s 吻合较好。计算结果表明，本文改进的铝粉反应模型适用于铝粉/空气两相爆轰的数值模拟。此外，根据爆轰波的C-J (Chapman-Jouguet)理论，D=uCJ+cCJ（uCJ和cCJ分别表示气体速度和气体声速），爆轰波CJ 参数分别为：pCJ=1.47 MPa，ρCJ=1.81 kg/m3，uCJ=534 m/s，TCJ=3 899 K。根据数值模拟，CJ 面离爆轰波阵面1.58 m，并且CJ 面上剩余铝粉的体积分数为初始值的0.7%。这说明两相爆轰的反应区远大于气相爆轰的反应区。

图5 为计算得到的t=3.83 ms 时刻爆轰管轴上两相的轴向速度、气体温度以及气相组分和铝粉的质量浓度分布。模拟结果表明，前导激波过后，气体速度急剧升高，铝颗粒在气体拖拽力的作用下也提高速度。但是由于固体颗粒具有较大的惯性，所以铝颗粒的速度下降较慢，并且在离波阵面较远处铝颗粒的速度甚至高于气体速度。

图5 3.83 ms 时刻爆轰管轴上两相的轴向速度、气体温度以及氮气、氧气和铝粉的质量浓度分布Fig. 5 Longitudinal distributions of two phases velocities in the x direction, gas temperature, mass concentrations of nitrogen, oxygen and Al particles along the tube axis at t=3.83 ms

在激波压缩和铝颗粒燃烧放热的作用下，气体温度逐渐升高，但是并不会无限制的升温，而是限制在约4 000 K，这正是新模型考虑Al2O3在高温下分解反应吸热的缘故，该数值模拟结果符合Glassman 观察到的实验现象[18]。此外，气体温度到达最高点后，会随着距爆轰波阵面距离的增加而有所降低，这是因为Al2O3的沸点随着压力的降低而降低[19]，而它的沸点温度决定了分解吸热反应是否发生以限制气体的最高温度。

同样地，氮气、氧气以及铝粉尘的浓度在前导激波过后迅速升高，此后铝粉和氧气由于燃烧反应快速消耗。根据数值模拟结果，3.83 ms 时刻5.5 m 位置处剩余铝粉尘的浓度为初始值的12.5%。

3.2 新模型的改进效果

为了体现新模型的改进效果，本节中采用改进前不考虑氧化铝分解反应的原模型对3.1 节中的工况铝粉/空气两相爆轰进行数值计算，并和上节的计算结果进行对比，分析模型改进后的数值仿真效果。

图6 为采用两个模型计算得到的相同时刻爆轰波附近气体温度和气体压力的分布。由图6 可以看出，采用改进前的模型计算得到的爆轰波后的气体温度高达4 800 K，该温度明显高于Al2O3的沸点。考虑到Al2O3不能以气体形式存在，所以在铝粉的燃烧模型中加入Al2O3的高温分解反应是合理并且必要的。此外，改进前的模型由于无法限制铝粉的燃烧放热，所以采用其计算得到的爆轰速度和爆轰压力均相对较高。

图6 模型改进前后计算得到的爆轰波附近气体温度和气体压力的分布Fig. 6 Distributions of gas temperature and gas pressure around the steady detonation wave calculated by the improved and previous models

基于以上对比分析可以看出，本文改进后的模型可以反映悬浮铝粉尘爆轰过程中存在的Al2O3高温分解现象，能够较准确地模拟计算爆轰流场。

3.3 铝粉/氧气两相爆轰

为了检验改进的铝粉反应模型对不同氧化性气氛的适用性，对铝粉/氧气混合物的两相爆轰开展了数值模拟。混合物由浓度为700 g/m3、直径为2.5 µm 的球形铝粉和纯氧气组成[23]，爆轰在直径为0.12 m的管道中进行。

数值计算中，仍采用1/4 管道模型作为计算域，其几何形状和图2 相似，不同之处是该计算模型的横截面为60 mm×60 mm，长度为10 m。在y和z方向上分为60 个网格，在x方向上有10 000 个网格，即3 个维度的网格尺寸均为1 mm。边界条件和上文铝粉/空气两相爆轰的计算模型相似。为了提高计算效率，数值模拟在250 个进程上并行计算，每个进程计算0.04 m 长的计算域。计算域左端的0.15 m 长为点火区域，铝粉/氧气混合物的参数为[23]：ρ1=1.28 kg/m3, σ2=700 g/m3,T1=T2=300 K, |u|=0, |v|=0。

图7 为计算得到的t=4.43 ms 时刻爆轰波阵面附近0.24 m 长的局部区域气体压力、气体温度、气体密度、气体轴向速度、铝粉体积分数以及铝粉轴向速度的分布。由图7 可以看出，此时爆轰波阵面到达7.27 m 位置，铝氧爆轰的前导激波后的气体温度很高，表明铝粉已开始燃烧，说明铝粉的诱导点火区域很短，紧跟着前导激波阵面即开始燃烧反应。这与上文的铝粉/空气两相爆轰有所不同，原因是该工况下的铝颗粒较小，表面质量比相对较大，所以能够更快地吸收热量达到点火温度，从而开始燃烧放热。从铝粉体积分数的云图还可看出，波阵面后约0.03 m 之外铝粉几乎燃烧完毕，由于铝颗粒较小以及纯氧环境，该爆轰反应区长度远小于上节中铝粉/空气爆轰约1.60 m 的反应区长度。

图7 4.43 ms 时刻铝粉/氧气爆轰波阵面附近气体压力、气体温度、气体密度、气体轴向速度、铝粉体积分数以及铝粉轴向速度的分布Fig. 7 Distributions of gas pressure, gas temperature, gas density, gas velocity in the x direction as well as volume fraction of Al particles, and Al particles velocity in the x direction around the detonation wave for Al/O2 mixtures at t=4.43 ms

数值模拟得到的该两相爆轰波速度为1 610 m/s，比文献值1 702 m/s[23]小5.4%，这可能是由于本文的数值模拟中考虑了管壁摩擦，而文献中没有考虑该摩擦作用。根据计算结果，pCJ=2.07 MPa, ρCJ=2.63 kg/m3,uCJ=543 m/s,TCJ=3838 K，CJ 面离爆轰波阵面0.93 m，并且CJ 面上没有剩余铝粉，和铝粉/空气两相爆轰相比，该爆轰的CJ 反应区更短。

图8 为计算得到的t=4.43 ms 时刻爆轰管轴上两相的轴向速度、气体温度以及铝粉、氧气的质量浓度分布。由图8 可以看出，前导激波过后，气体和铝颗粒相继加速，在约7.18 m 位置处铝颗粒燃烧完毕，表明铝粉反应区约0.09 m 长。由气体温度分布可以发现，在极小的区域内气体温度即上升至最高点，这一现象和上文的铝粉/空气爆轰的情况不同，这是较小的铝颗粒在更短的时间被点燃并燃烧放热的缘故。

图8 4.43 ms 时刻爆轰管轴上两相的轴向速度、气体温度以及铝粉、氧气的质量浓度分布Fig. 8 Longitudinal distributions of two phase velocities in the x direction, gas temperature,mass concentrations of oxygen and Al particles along the tube axis at t=4.43 ms

4 结 语

铝粉反应模型是悬浮铝粉尘两相爆轰数值模拟研究的难点。本文中，通过引入铝粉燃烧产物氧化铝（Al2O3）在高温下的分解反应，改进了铝粉的燃烧反应模型，并且新的模型更符合物理现象。分别对铝粉/空气混合物爆轰以及铝粉/氧气混合物爆轰进行数值模拟，计算的结果和实验结果、文献值均吻合较好，表明改进的铝粉反应模型适用于不同氧化性气氛中悬浮铝粉尘爆轰的数值模拟。