基于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）