铝粉尘云团爆轰温压效应的数值模拟

昝文涛， 洪滔， 董贺飞

(1.北京理工大学 机电学院， 北京 100081； 2.北京应用物理与计算数学研究所， 北京 100094)

0 引言

温压武器作为一种新型武器近些年来得到了广泛的关注，其利用高温、高压可以造成大范围杀伤效果，铝粉尘由于其质量轻、含能高的特点被广泛用于温压武器固体药剂中。而近些年来铝粉尘在工业领域造成的爆炸事故也越来越多，造成了巨大的人员伤亡及经济损失，因此研究铝粉尘云团的爆轰温压效应有着重要的意义。

裴明敬等[1]通过实验研究了含铝温压燃料的能量释放效率及爆炸冲击波的形成及扩散过程，发现铝粉对冲击波能量贡献很大，能量释放率较高；郑波等[2]通过高速分析系统观察了温压炸药的抛洒过程，分析了温压炸药抛洒规律；李秀丽等[3]通过外场实验研究了云爆剂的爆炸冲击波参数，发现装药量30 kg的实验弹爆炸火球直径达到17.4 m，并存在两个正压作用区；洪滔等[4-6]研究了铝颗粒的点火机制，分析了铝颗粒的点火温度及悬浮铝粉尘爆轰波参数，并模拟了爆轰波管中的铝粉尘爆轰过程。

由于温压武器毁伤过程包含抛洒、起爆等复杂过程，目前对于温压燃料杀伤效果的研究主要以实验为主，对于其毁伤效应的数值模拟基本以商业软件模拟含铝温压炸药爆炸为主。

为此，本文中采用时空守恒元解元(CE/SE)算法求解两相流方程组，编制程序模拟了半径3 m的均匀铝粉尘云团在空气中爆炸后的冲击波传播过程及产生的高温火球特性，分析了铝粉尘云团爆轰的温压效应。

1 模型

模拟过程采用了两相流模型，假设颗粒为球形，初始直径都相同，单个颗粒的温度都是均匀的，忽略粒子间的作用，粒子化学反应产生的能量假定都被气体吸收，气体的组分都是均匀混合的，忽略粒子与气体间的辐射作用，固相和气相都满足各自守恒方程[7-9]。

气相守恒方程:

(1)

(2)

(3)

(4)

式中：下角标g、a分别代表气体和铝颗粒;ρ为密度；u为速度；e为内能；p为压力；φ为气体和固体颗粒体积分数，φg+φa=1；qa为单位质量的铝粒子反应能；Q为气体与粒子间的热传导，

(5)

λ为气体导热系数，Nu为Nusselt数。

固相守恒方程:

(6)

(7)

(8)

(9)

(10)

式中：I为单位体积内颗粒的质量变化率；n为单位体积内粒子数；F为气体对粒子的压力，

(11)

(12)

Cd为拖曳系数，

(13)

Re为雷诺数。

气体组分方程:

(14)

式中：j为气体组分O、N、Al；y为组分浓度。

气体状态方程：

(15)

式中：w为分子量；R为气体常数；T为气体温度；k为气体中组分数量。

对于铝颗粒，采用如(16)式和(17)式所示的燃烧模型[11]，由于铝颗粒表面被Al2O3包裹，铝颗粒温度达到Al的熔点时铝颗粒膨胀发生破裂从而点火。

(16)

(17)

式中：r为粒子半径；Ti为点火温度；d0为粒子初始直径；Φ为气体中氧气的摩尔份额；m=1.75.

本文中所有物理量均采用国际单位制，模拟采用了CE/SE算法[12-14]求解欧拉方程组，对于两相流方程组的源项，使用4阶龙格- 库塔方法求解。

2 数值模拟结果

2.1 程序验证

对于理想气体，气体常数1.4，分别计算了网格数量为400、300、200、100和32的算例，计算结果如图1所示。

图1 t′=0.2时Sod激波管压力Fig.1 Pressures of Sod shock tube for t′= 0.2

从图1中可以看出，网格数量200时可以很好地模拟Sod问题，与理论解符合较好，网格数量继续减小到32时也能模拟Sod激波管问题，但是与理论解的相比，波阵面宽度变宽了。

模拟了铝粉尘在管道内的爆轰波传播过程，铝粉尘密度为0.304 kg/m3，颗粒半径为1.7 μm，管道直径为15.2 cm，从左端起爆，左端闭合，右端开口，上下为固壁，模拟得到的爆轰波速度D=1 630 m/s，详见文献[7]， Tulis等[15]由实验中得到的铝粉尘爆速为1 650 m/s.

2.2 计算条件

本文数值模拟了半径为3.0 m的铝粉尘云团均匀分布在空气中。铝粉尘密度为0.304 kg/m3，与分布范围内空气中氧气当量比为1，颗粒半径为2 μm，铝密度为2 700.000 kg/m3. 空气初始密度为1.210 kg/m3，气体压力为0.10 MPa，温度为298 K. 起爆在中心位置起爆，起爆区为2×2个网格，起爆条件φg=1，ρg=2.200 kg/m3，ug=1 400 m/s，vg=1 400 m/s，Tg=3 200 K. 由于模型中心对称，为减小计算量，采用计算模型1/4区域，左侧与下侧采用全反射对称边界，上侧与右侧采用透射边界。如图2所示，左下角区域为半径3.0 m铝粉尘云团，计算模型尺寸为25 m×25 m，网格数量800×800. 计算中采用的参数[16]：λ=0.1 J/(m·s·K)，qa=ηQa，Qa=31.5 kJ/g，η=80%.

图2 模拟区域尺寸模型Fig.2 Geometry model of simulation area

2.3 数值模拟结果

2.3.1 不同时刻爆炸流场物理量分布

图3为t=2.13 ms时爆炸流场物理量分布。图3(a)为爆炸流场压力分布，此时冲击波到达3.0 m位置，到达铝粉尘云团边缘，由于起爆区域为正方形，因此向上和向右传播的冲击波阵面接近平面波，最大压力为2.10 MPa左右，中间部分冲击波阵面为圆弧形，压力为1.75 MPa，波后流场压力逐渐下降。图3(b)为爆炸流场密度分布，向上和向右传播的冲击波阵面密度较高，约为2.900 kg/m3，中间部分冲击波波阵面密度为2.600 kg/m3，波后流场密度逐渐下降，中心位置密度较低，为1.000 kg/m3. 图3(c)为流场温度分布，可以看出流场中冲击波阵面处温度迅速上升，波后流场中大部分区域温度达到3 800 K. 图3(d)是流场中铝粉尘体积分数分布，此时铝粉尘由于激波作用堆积在波阵面附近，冲击波阵面上侧与右侧处体积分数约为5.2×10-5，与初始时刻浓度比为46.2%，右上部体积分数约为2.0×10-5，与初始时刻浓度比为17.8%，中心处铝粉尘基本反应完毕剩余较少。

图4为t=3.49 ms时爆轰流场物理量分布，此时冲击波到达4.7 m位置。图4(a)是冲击波流场压力分布，此时向上和向右传播的冲击波阵面压力下降到0.90 MPa，中间部分冲击波阵面压力为0.75 MPa，波后压力缓慢下降，而后小幅上升达到稳定，中心区域压力稳定在0.75 MPa. 图4(b)为爆轰流场密度分布，向上和向右传播的冲击波阵面处密度为3.800 kg/m3，中间部分冲击波阵面处密度为3.200 kg/m3，波后流场密度下降至0.800 kg/m3，中心区域密度约1.300 kg/m3. 图4(c)为流场温度分布，冲击波附近温度逐渐上升到3 800 K，对比图4(b)发现此部分区域同时为低密度区域，将此部分高温低密度区域定义为火球区域，将2 500 K处设为火球边界，可以看出此时火球与冲击波开始发生分离。图4(d)是流场中铝粉尘体积分数分布，极少量铝粉尘在冲击波作用下运动到4.0 m处，可以看出此时铝粉尘基本反应完毕，只有在冲击波后上侧与右侧部分区域存在少量铝粉尘，浓度与初始时刻比例约0.14%. 此时铝粉尘速度为480 m/s，冲击波速度为900 m/s，由于冲击波波速较铝粉尘颗粒速度快，因此此时冲击波到达4.7 m位置而铝粉尘颗粒只到达4.0 m处。当冲击波离开云团时，冲击波阵面与粉尘云团分离，反应区的能量无法使冲击波保持增长，冲击波压力降下降。

图3 t=2.13 ms时刻的流场压力、密度、温度和铝粉尘体积分数Fig.3 Pressure, density and temperature of flow field, and volume fraction of aluminum particles for t=2.13 ms

图4 t=3.49 ms时刻的流场压力、密度、温度和铝粉尘体积分数Fig.4 Pressure, density and temperature of flow field, and volume fraction of aluminum particles for t=3.49 ms

图5为t=10.99 ms时刻流场的压力、密度和温度分布图，可以看出此时冲击波到达10.0 m处。图5(a)为流场压力分布图，冲击波阵面压力衰减到0.36 MPa，波后区域压力逐渐下降，中心区域压力约0.17 MPa. 图5(b)为流场密度分布图，冲击波阵面处密度约为2.700 kg/m3，冲击波后密度逐渐下降，中心区域密度约为0.400 kg/m3，此区域半径约为6.5 m，对比图5(c)可以看出此部分区域为火球区域。图5(c)为流场温度分布图，冲击波后流场温度上升至800 K左右，在火球处温度迅速上升，火球中心区域温度在3 800K左右，火球与冲击波产生明显分离。图5(d)为流场中O2体积分数φO分布，可以看出火球内部中心区域中O2体积分数只有0.02左右，火球内部靠外区域中O2体积分数稍高为0.03，火球外O2体积分数为0.23. 在铝粉尘爆炸过程中会产生气态铝，而在火球扩散过程中气态铝与O2发生反应释放能量，维持火球中心的高温。

图5 t=10.99 ms时刻的流场压力、密度、温度和氧气体积分数Fig.5 Pressure, density and temperature of flow field,and volume fraction of O2 for t=10.99 ms

图6 t=29.32 ms时刻的流场压力、密度、温度和氧气体积分数Fig.6 Pressure, density and temperature of flow field, and volume fraction of O2 for t=29.32 ms

图6为t=29.32 ms时刻流场物理量分布。图6(a)为流场压力分布图，此时冲击波传播至20.0 m处，冲击波阵面压力衰减至0.23 MPa，冲击波后流场压力逐渐下降，中心区域压力为0.10 MPa. 图6(b)为流场密度分布图，冲击波阵面密度下降至2.100 kg/m3，波后流场密度逐渐下降，到达中心火球区域时密度约为0.200 kg/m3. 图6(c)为流场温度分布图，冲击波后流场中温度约800 K，此时火球传播至8.0 m位置，火球附近流场温度由800 K迅速上升至3 000 K，由于火球进一步扩散，火球内部温度下降到3 500 K以上。图6(d)为流场中O2体积分数分布，流场中O2含量在火球处急剧下降，从火球外0.23下降至0.02左右，可以看出火球区域为一个低压、低密度、低氧、高温区域。

2.3.2 流场中不同位置的物理量随时间的变化曲线

图7 数据采集点及不同时刻流场近场压力和铝粉尘体积分数曲线Fig.7 Data collection points, and pressure in near field and volume fraction of aluminum distribution lines at different times

选取对角线上的点作为研究对象，图7(a)为数据采集点示意图，图7(b)和图7(c)为不同时刻近场压力曲线以及铝粉尘体积分数曲线。起爆后，初始冲击波向外传播，铝粉尘被点火发生反应，冲击波压力在3.0 m处粉尘云团边界位置时上升达到最大值1.75 MPa，随冲击波继续向外传播，云团内铝粉尘被气体加速向外扩散并消耗减少，波后压力逐渐下降。随时间增加，如图7(c)所示,铝粉尘云团可加速扩散至3.5 m处，极少量粉尘颗粒可以到达4.0 m处反应完毕。

图8 不同时刻的流场压力、密度和温度曲线Fig.8 Pressure, density and temperature distributions at different times

图8(a)为不同时刻远场的流场压力曲线，可以看出在小于4.0 m时，铝粉尘反应释放能量，流场压力较高，4.0 m以后冲击波压力迅速下降。到达28.0 m位置时，冲击波阵面压力下降至0.19 MPa. 图9(a)为冲击波和火球随时间传播距离s曲线，图9(b)为冲击波和火球随时间传播速度曲线，从图中可以看出，冲击波传播速度随时间急速下降，45 ms时刻冲击波传播速度下降至500 m/s，此时冲击波传播至28.0 m处。根据冲击波对暴露人员的损伤程度[17]，冲击波超压达到0.07 MPa时就会对人体产生严重杀伤，达到0.1 MPa时会造成50%致死。根据图8(a)，冲击波在24.5 m处时冲击波阵面压力下降至0.20 MPa，冲击波在28.0 m处时波阵面压力为0.19 MPa，因此可以认为半径24.5 m范围内致人死亡，24.5～28.0 m范围内对人体造成严重损伤。

图9 不同时刻冲击波与火球的传播距离和速度Fig.9 Propagation distances and velocities of shock wave and fireball at different times

图8(b)和图8(c)为不同时刻流场密度与温度曲线。起爆后，冲击波阵面处密度随冲击波传播逐渐上升，而后开始下降；冲击波波后流场密度逐渐下降，到达火球区域时密度产生急速下降，可以看出火球区域密度逐渐稳定在0.150 kg/m3左右。从不同时刻温度曲线可以看出随时间增大，火球扩散速度降低，从图8(b)可以看出，33 ms以后火球传播速度基本接近为0 m/s，到45 ms时，火球半径为10.0 m，火球中心温度在3 500 K以上，在半径10.0 m范围内以高温及超压杀伤为主。

图10 t=43.87 ms时刻沿对角线处流场参数分布Fig.10 Parameters of flow field at the points of diagonal for t=43.87 ms

图10为43.87 ms时刻沿对角线区域的流场参数曲线。图10(a)为密度分布曲线，此时冲击波到达27.0 m位置，冲击波阵面上密度为1.900 kg/m3，波后流场密度下降，在火球区域密度急速下降至0.120 kg/m3. 图10(b)为压力分布曲线，冲击波阵面处压力降为0.19 MPa，波后流场压力下降。在冲击波后压缩波，在实验中也观察到了此现象[3]。在对图11中流场压力分布随时间变化曲线分析后得出，在冲击波与火球发生分离时，波后压力下降，而火球内由于铝粉尘及气态铝反应放热温度较高，压力较高，因此冲击波后压力下降后会产生小幅跃升，随冲击波传播距离增大，火球内部铝反应完毕，火球内中心区域压力下降，因此在火球边缘处形成一个压缩波，而火球传播速度随时间增大而急剧下降，造成压缩波与火球发生脱离，在冲击波后形成二次波，在火球边缘附近压力有小幅跃升而后下降。图10(c)为流场气体速度分布曲线，冲击波阵面处气体速度165 m/s，波后流场中气体速度下降，在压缩波处气体速度为70 m/s，在火球边缘处气体速度由10 m/s跃升至30 m/s，而后气体速度下降。图10(d)为流场中温度曲线，冲击波阵面后气体温度在350 K左右，到达火球处时气体温度迅速上升，火球中心温度达到3 500 K以上。

图11 不同时刻的局部流场压力曲线Fig.11 Pressure distributions of local flow field at different times

3 结论

本文数值模拟了当量比为1、密度为0.304 kg/m3、颗粒半径为2 μm的铝粉尘在空气中均匀分布形成的半径3.0 m的云团，在被中心起爆后对周围环境的温压效应，得出以下结论：

1)初始冲击波在粉尘云团中传播并引起铝粉尘点火发生化学反应形成爆轰波，波阵面压力在冲击波到达云团边界3.0 m处时达到最大值1.75 MPa. 随冲击波传播的距离增加，冲击波压力下降, 这是由于冲击波与云团分离，粉尘中的能量无法支持冲击波造成的。

2) 铝粉尘反应形成的火球随时间可向外扩散到达10.0 m位置，火球内部为密度0.120 kg/m3、温度3 500 K以上的低密度、低氧区域。

3) 铝云团爆轰在半径10.0 m区域内对周边毁伤以高温及冲击波超压为主，在半径24.5 m区域内冲击波超压可致人死亡，24.5～28.0 m范围内可对人体造成严重杀伤。

本文对空气中悬浮铝粉尘云团的爆轰进行了数值模拟，分析了自由场中铝粉尘云团爆轰对周边环境的杀伤规律，对模拟温压弹的杀伤及铝粉尘爆炸安全事故对周边环境的杀伤及防护提高了认识。

)

[1] 裴明敬, 毛根旺, 胡华权, 等, 含铝温压燃料性能研究[J]. 含能材料, 2007, 15(5):442-446.

PEI Ming-jing, MAO Gen-wang, HU Hua-quan, et al. Characteristic of the thermobaric explosive contained aluminum powders[J]. Chinese Journal of Energetic Materials, 2007, 15(5):442-446.(in Chinese)

[2] 郑波, 陈力, 丁雁生, 等, 温压炸药爆炸抛洒的运动规律[J]. 爆轰与冲击, 2008, 28(5): 433-437.

ZHENG Bo, CHEN Li, DING Yan-sheng, et al. Dispersal process of explosion production of thermobaric explosive[J]. Explosion and Shock Waves, 2008, 28(5):433-437.(in Chinese)

[3] 李秀丽, 惠君明, 王伯良. 云爆剂爆炸/冲击波参数研究[J]. 含能材料, 2008, 16(4):410-414.

LI Xiu-li, HUI Jun-ming, WANG Bo-liang. Blast/shock wave parameters of single-event FAE[J]. Chinese Journal of Energetic Materials, 2008, 16(4) :410-414.(in Chinese)

[4] 洪滔, 秦承森. 铝颗粒激波点火机制初探[J]. 爆轰与冲击, 2003, 23(4):295-299.

HONG Tao, QIN Cheng-sen. Mechanism of shock wave ignition of aluminum particle[J]. Explosion and Shock Waves, 2003, 23(4):295-299.(in Chinese)

[5] 洪滔, 秦承森. 悬浮铝粉尘爆轰波参数[J]. 含能材料, 2004, 12(3):129-133.

HONG Tao, QIN Cheng-sen. Parameters of detonation in suspended aluminum dust[J]. Chinese Journal of Energetic Materials, 2004, 12(3):129-133.(in Chinese)

[6] 洪滔, 秦承森. 爆轰波管中铝粉尘爆轰的数值模拟[J]. 爆轰与冲击, 2004, 24(3):193-200.

HONG Tao, QIN Cheng-sen. Numerical simulation of dust detonation of aluminum powder in explosive tubes[J]. Explosion and Shock Waves, 2004, 24(3):193-200.(in Chinese)

[7] 昝文涛, 洪滔, 董贺飞. 基于CE/SE方法关于RDX-AL悬浮粉尘在空气中的两相爆轰的数值模拟[J]. 爆炸与冲击, 2016, 36(5):603-610.

ZAN Wen-tao, HONG Tao, DONG He-fei. Numerical simulation of two phase detonation of suspending RDX-AL dust in air with CE/SE method[J]. Explosion and Shock Waves, 2016, 36(5):603-610.(in Chinese)

[8] 昝文涛, 洪滔, 董贺飞. 带管道连接的空间中悬浮铝粉尘爆轰波传播数值模拟[J]. 含能材料, 2017, 25(6): 508-514.

ZAN Wen-tao, HONG Tao, DONG He-fei. Numerical simulation of detonation of suspending aluminum dust in air in the universal spaces connected by channel[J]. Chinese Journal of Energetic Materials, 2017, 25(6):508-514.(in Chinese)

[9] Zan W T, Hong T, Dong H F. Simulation of double-front detonation of suspended mixed RDX and aluminum dust in air[J]. Chinese Physics Letter, 2017, 34(7): 074701.

[10] Steinberg T A, Wilson D B, Benz F. The combustion phase of burning particle[J]. Combustion and Flame, 1992, 91(2):200-208.

[11] Price E W. Combustion of metalized propellants[C]∥Progress in Astronautics and Aeronautics: Fundamenals of Solid-Propellant Combustion. NY, US:AIAA, 1984:479-513.

[12] Chang S C. The method of space-time conservation element and solution element-a new approach for solving the Navier-Stokes and Euler equations[J]. Journal of Computational Physics, 1995,119(2): 295-324.

[13] Zhang D L, Wang J T, Wang G. High-order CE/SE method and applications[J]. Chinese Journal of Computational Physics, 2009, 26(2):211-220.

[14] Wang J T, Zhang D L, Liu K X. A Eulerian approach based on CE/SE method for 2D multimaterial elastic-plastic flows[J]. Chinese Journal of Computational Physics, 2007, 24(4):395-401.

[15] Tulis A J, Selman J R. Detonation tube studies of aluminum particles dispersed in air[J]. Symposium (International) on Combustion, 1982, 19(1):655-663.

[16] 洪滔, 秦承森, 林文洲. 悬浮RDX炸药和铝颗粒混合粉尘爆轰的数值模拟[J]. 爆炸与冲击, 2009, 29(5):468-473.

HONG Tao, QIN Cheng-sen, LIN Wen-zhou. Numerical simulation of detonation in suspended mixed RDX and aluminum dust[J]. Explosion and Shock Waves, 2009, 29(5):468-473.(in Chinese)

[17] GJB 5212—2004 云爆弹定型实验规程[S]. 北京: 中国人民解放军总装备部, 2004.

GJB 5212—2004 Finalizing test procedures for fuel-air-explosive ammunition[S]. Beijing: General Armament Department of PLA, 2004.(in Chinese)