基于CE/SE方法的铝粉尘爆轰二维两相数值计算

韦 伟，翁春生

（1.南京理工大学 瞬态物理国家重点实验室，南京210094；2.江苏海事职业技术学院，南京211170）

铝粉因其特殊的工业与军事应用背景，其燃烧与爆炸性能一直受到密切关注.研究铝粉燃烧转爆轰（DDT）的性能，对提高武器的威力及性能具有重要意义.然而该过程中，燃烧波的形成、加速以及转变为爆轰波的机理十分复杂.国外的STRAUSS[1]、TULIS[2]等，国内的刘晓利[3]、李小东等[4]采用实验的方式研究铝粉在氧气或空气中爆轰的现象.同时在数值计算方面，WOLANSKI P[5]、VEYSSIERE B[6]、FEDOROV A V[7]等用两相流体力学模型研究铝粉尘爆轰波的定常问题.而在非定常问题上，洪涛[8]用二阶精度的 MacCormack差分格式加上FCT技术在求解一维铝粉尘爆轰波方面有了一定的突破.本文采用CE/SE方法对铝粉/空气两相爆轰进行了二维数值模拟，不仅考虑了爆轰波的非定常性，还考虑了气相的组份变化，计算精度高、捕获爆轰波能力强.

1 管内铝粉尘爆轰理论模型

管内铝粉尘爆轰过程异常复杂，为简化问题，模型做以下假定：①所研究的流场是二维轴对称非定常流场；②铝粉尘颗粒在气体中均匀分布，并将其视为具有连续介质特性的拟流体；③铝粉尘颗粒均为初始直径相同的球形，并且颗粒内部温度均匀分布；④忽略铝粉尘颗粒间的相互作用力以及与壁面的作用；⑤燃烧产物氧化铝粉末视作为气相产物，但不考虑其对压力的影响；⑥不考虑反应的中间产物，放能与压力无关，化学反应释放的能量全部被气体吸收[8～11].模型中考虑了气相、固体颗粒相的质量、动量和能量的交换，以及气体内部O2、N2、Al2O3等的组分浓度变化.相关质量守恒方程、组分方程、动量守恒方程、能量方程为

式中，下标“g”表示气相，下标“s”表示固相；φ、ρ、v、u分别为体积分数、密度、径向速度和轴向速度，φg＋φs＝1；Id为单位体积内铝粉颗粒因燃烧引起的质量变化率；wO和wA分别为O2和Al2O3的组分浓度，它们与N2的组分浓度一起满足归一化条件；Fdx、Fdy分别为气相与固体颗粒间轴向、径向的相互作用力；p为气相压力；Eg为气相能量；Es为固体颗粒相能量；Qd为气相与固体颗粒相间的对流传热量；Qc为固体颗粒相燃烧释放的热量.

2 理论模型的求解方法

本文应用二维守恒元－求解元方法求解，其思想、求解方程和加权函数的取得见参考文献[12].在此仅列出离散量的表达式[12]：

本文推导了满足A＝∂F/∂U，B＝∂G/∂U的10×10的雅克比系数矩阵表达式，亦可将Ft和Gt用雅克比矩阵表示为U、Ux和Uy的函数，即[9]：

具体雅克比系数矩阵在此不赘述.

在源项的处理上应用四阶龙格－库塔方法求解，时间步长为[10]

式中，ΔtCE为二维CE/SE方法的计算时间步长，而N一般取为5～20.

3 计算条件及结果分析

由于爆轰管是轴对称的，计算区域取爆轰管的一半.爆轰管长为1 000mm，半径为30mm.初始时刻爆轰管内按化学当量比均匀充满常温、常压的铝粉颗粒/空气混合物，颗粒直径取1μm，两相速度均为0.

通过计算得到不同时刻爆轰管内中心轴线处的压力、速度和温度沿轴向的分布曲线图，如图1所示.由图可知，铝粉被点燃后，随着能量的释放，燃烧波压力值不断上升，速度不断加快，最终形成稳定的爆轰波.通过计算，稳定爆轰波的传播速度为1 285m/s.对比图1（a）和图1（c）可知，爆轰管内火焰阵面（以温度表示）和压力阵面耦合较好，说明在爆轰管内铝粉和空气能够形成自持传播的爆轰波.图1（a）还显示，爆轰波阵面一般只有3～4个网格，很陡.可见，在捕获如爆轰波这类强间断波方面，CE/SE方法是切实、有效的.

图2为爆轰管中心轴线0.5m处的气相温度和颗粒相半径随时间的变化图.由图可见，爆轰波到达后首先对未燃混合物进行加热，t1＝347.54μs时刻，气相温度开始上升.随后t2＝349.19μs时刻，铝粉开始燃烧，颗粒半径开始减小.铝粉燃烧释放的能量持续使爆轰波的强度增强，温度在t3＝365.36μs达到峰值2 936K，此时的铝粉半径为0.15μm.在爆轰波扫过后的一段时间内铝粉燃烧依然持续，直到t4＝426.67μs时燃烧完全.

在爆轰过程中气相组份发生相应的变化.初始时刻气相组份中O2占21%，N2及其它不参与反应的气体占79%.爆轰过程的化学反应产物Al2O3视作气相部分，初始时刻其组份含量为0.图3为Al2O3在爆轰管0.798 m处含量的变化图.在t＝506.7μs时刻之前，由于爆轰波尚未传到此处，气相只含有O2和N2等；在t＝506.7μs时刻之后，气相产物Al2O3的含量逐渐增加.图4为t＝10.55μs时刻，即爆轰初始时刻，管内铝粉颗粒含量变化图.由图可见，在距离封闭端很短的范围内（x＝0.005m），由于燃烧不充分，铝粉含量较高.但随着快速的化学反应，铝粉含量下降速度很快.在x＝0.01m之后，铝粉含量降为原有的14%.

图1 爆轰管内轴向压力、速度、温度随时间变化图

图2 中心轴线0.5m处的气相温度、铝粉半径随时间的变化图

图3 爆轰管出口处Al2O3百分含量变化图

图4 t＝10.55μs时刻铝粉颗粒体积分数变化图

4 结论

本文用CE/SE方法对铝粉尘/空气二维轴对称两相管内爆轰过程进行了数值模拟.计算结果表明，在捕获如爆轰波这类强间断波方面，CE/SE方法是切实、有效的.

计算结果还反映出，铝粉和空气被点燃后形成的燃烧波通过与爆轰管壁面碰撞后形成反射波和局部强爆轰点，加快铝粉燃烧速率，经过燃烧向爆轰的转变，最终形成爆轰波.爆轰管内火焰阵面和压力阵面能很好地耦合，爆轰管内铝粉和空气能够形成自持传播的爆轰波.

[1]STAUSS W A.Investigation of the detonation of aluminum powder－oxygen mixtures[J].AIAA，1968，6（9）：1 753－1 756.

[2]TULIS A J.Detonation tube studies of aluminum particles dispersed in air[C].19th International Symposium on Combustion.USA：Standing Committee on Combustion Symposia，1982：633－655.

[3]刘晓利，李鸿志，郭建国，等.铝粉－空气混合物燃烧转爆轰过程的实验研究[J].爆炸与冲击，1995，15（3）：218－228.LIU Xiao－li，LI Hong－zhi，Guo Jian－guo，et al.An experimental investigation of deflagration to detonation transition in aluminum dust－air mixture[J].Explosion and Shock Waves，1995，15（3）：218－228.（in Chinese）

[4]李小东，刘庆明，白春华，等.铝粉－空气混合物的燃烧转爆轰过程[J].火炸药学报，2009，32（6）：58－61 LI Xiao－dong，LIU Qing－ming，BAI Chun－hua，et al.Deflagration to detonation transition process in aluminum dust－air mixture[J].Journal of Explosives ＆ Propellants，2009，32（6）：58－61.（in Chinese）

[5]WOLANSKI P，LEE D，SICHEL M.The structure of dust detonations[C].Dynamics of Shock Waves，Explosions and Detonations：AIAA Progress in Astronautics and Aeronautics.New York：AIAA，1984：214－263.

[6]VEYSSIERE B，KHASAINOV B A.Steady plane double－front detonations in gaseous detonable mixtures containing a suspension of aluminum particles[J].Combustion and Flame，1991，85：241－253.

[7]FEDOROV A V.Non－equilibrium model of steady detonations in aluminum particles－oxygen suspensions[J].Shock Waves，1999，9（5）：313－318.

[8]洪涛，秦承森.爆轰波管中铝粉尘爆轰的数值模拟[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）

[9]马丹花，翁春生.二维守恒元和求解元方法在两相爆轰流场计算中的应用[J].燃烧科学与技术，2010，16（1）：85－91.MA Dan－hua，WENG Chun－sheng.Application of two－dimensional CE/SE method to calculation of two－phase detonation flow field[J].Journal of Combustion Science and Technology，2010，16（1）：85－91.（in Chinese）

[10]白桥栋，翁春生.二维CE/SE方法在内弹道湍流两相流动中的应用[J].兵工学报，2009，30（6）：682－687.BAI Qiao－dong，WENG Chun－sheng.Application of two dimensional space－time conservation element and solution element method（CE/SE）to turbulent two－phase flow in interior ballistics[J].Acta Armamentarii，2009，30（6）：682－687.（in Chinese）

[11]彭振，翁春生.等离子体点火对燃烧转爆轰影响的数值计算[J].燃烧科学与技术，2011，17（1）：84－89.PENG Zhen，WENG Chun－sheng.Numerical calculation of effect of plasma ignition on DDT[J].Journal of Combustion Science and Technology，2011，17（1）：84－89.（in Chinese）

[12]翁春生，王浩.计算内弹道学[M].北京：国防工业出版社，2006.WENG Chun－sheng，WANG Hao.Computational Interior Ballistics[M].Beijing：National Defense Industry Press，2006.（in Chinese）