鋁粉塵爆轟的反應(yīng)模型*

岳軍政，洪 滔，吳先前，黃晨光

（1. 中國科學院力學研究所，北京 100190；2. 北京應(yīng)用物理與計算數(shù)學研究所，北京 100094）

鋁粉作為一種常見的工業(yè)原料，在其生產(chǎn)過程中容易與空氣混合發(fā)生兩相爆轟造成巨大危害[1]，同時因其含能較高，且價格便宜，也常作為燃料加入燃料空氣炸彈(fuel air explosive, FAE)中形成大規(guī)模粉塵爆轟以制造破壞。針對鋁粉的工業(yè)生產(chǎn)安全及其國防應(yīng)用，對懸浮鋁粉塵的兩相爆轟已開展了大量研究[2-3]。由于實驗觀測難度大以及實驗重復性較低，數(shù)值模擬成為研究粉塵爆轟的重要手段。

懸浮鋁粉塵兩相爆轟是非理想爆轟過程，常用經(jīng)典的ZND (Zel’dovich-von Neumann-D?ring)爆轟模型描述其爆轟波結(jié)構(gòu)，其數(shù)值模擬的難點在于鋁粉的燃燒模型。在目前的研究中，鋁粉燃燒模型主要有擴散模型和動力學模型兩種。在擴散燃燒模型中，化學動力學反應(yīng)速率很快，燃燒速率主要受氧化氣體的擴散速率控制[4]；而在動力學燃燒模型中，鋁顆粒表面的氧化氣體擴散速度較快，主要由化學動力學反應(yīng)速率控制鋁顆粒的燃燒速度[5-6]。Glorian 等[7]考慮顆粒大小以及燃燒環(huán)境，對表面反應(yīng)動力學和氣相機制進行了理論分析。在現(xiàn)有的數(shù)值模擬工作中，Veyssiere 等[2]提出了一個兩步反應(yīng)模型，通過擬合動力學參數(shù)得到了與實驗結(jié)果一致的計算結(jié)果。Zhang 等[8]提出了同時考慮擴散機制和動力學機制的混合模型，并且討論了兩個燃燒機制的適用條件。此外，Briand 等[9]改進了兩步反應(yīng)模型，并且發(fā)展了一個混合模型。然而，在這些模型中有一個指前參數(shù)Zhyb需要擬合，而目前對這個參數(shù)并沒有統(tǒng)一的認識[10]。

本文中，擬通過考慮鋁粉燃燒產(chǎn)物氧化鋁(Al2O3)在高溫下的分解，針對鋁粉的擴散燃燒模型，增加燃燒產(chǎn)物的逆反應(yīng)吸熱機制。此外，將改進的鋁粉反應(yīng)模型嵌入到并行的粉塵爆轟數(shù)值模擬程序中，分別對鋁粉/空氣混合物以及鋁粉/氧氣混合物的兩相爆轟進行三維數(shù)值模擬，將計算得到的穩(wěn)定爆轟波速度與實驗結(jié)果和文獻值進行對比，獲得爆轟流場的物理量分布。

1 計算模型

1.1 模型假設(shè)

由于鋁粉塵氣-固兩相爆轟是復雜的物理化學過程，因此對計算模型作以下假設(shè)[11]：（1）鋁粉不可壓縮，且均勻分布在氣體中；（2）鋁粉為直徑相等的球形顆粒，并且片狀顆粒等效為球形；（3）單個鋁顆粒內(nèi)部沒有溫度梯度和壓力，重力、顆粒之間的作用以及熱輻射作用可忽略；（4）鋁粉燃燒所釋放的能量均被氣體吸收。

1.2 控制方程

數(shù)值模擬中采用兩相流模型[12]描述懸浮鋁粉塵爆轟的發(fā)展過程。模型中將氣、固兩相均看作連續(xù)介質(zhì)，氣體和固體滿足各自的守恒方程，并且考慮氣體和固體顆粒之間的質(zhì)量、動量和能量的交換。描述鋁粉塵爆轟的三維非定常兩相流模型方程如下。

氣相質(zhì)量、動量和能量守恒方程分別為：

固相質(zhì)量、動量和能量守恒方程分別為：

式中：下標1 和2 分別指氣相和固相鋁顆粒；ρ 為密度；φ為體積分數(shù)，且φ1+φ2=1；e為比內(nèi)能；u和v分別為氣相和鋁顆粒的速度矢量；p為氣體壓力；I為燃燒引起的單位體積內(nèi)鋁粉的質(zhì)量變化率；F為氣體對固體顆粒的作用力矢量；Q為兩相間的對流熱傳導速率；qAl為單位質(zhì)量鋁粉燃燒釋放的熱量。

顆粒數(shù)守恒方程為：

式中：n為單位體積內(nèi)的鋁顆粒數(shù)。

氣體組分方程為：

式中：yj為氣體組分j的質(zhì)量分數(shù)，ωj為單位體積內(nèi)氣體組分j的質(zhì)量變化率。

氣體狀態(tài)方程為：

式中：T為溫度，R為普適氣體常數(shù)，wj為氣體組分j的摩爾質(zhì)量。

兩相比內(nèi)能分別為:

式中：cV1為 氣體比定容熱容，隨氣體溫度T變化[13]；cV2為鋁的比熱容，e2的方程中已考慮相變。

氣體對鋁顆粒的作用力F表達式為[12]：

式中：r為鋁顆粒半徑，Cd為拖拽系數(shù)。拖拽系數(shù)的表達式為：

式中：Re為雷諾數(shù)，且Re= 2ρ1|u?v|r/μ，μ為氣體黏性系數(shù)。

兩相間的對流熱傳導速率為[12]：

式中：λ1為氣體導熱系數(shù)；Nu為Nusselt 數(shù)，且Nu=2+0.459Re0.55Pr0.33；Pr為普朗特數(shù)，且Pr=μcp/λ1，cp為氣體比定壓熱容。

式（13）中的因子β 是針對片狀鋁顆粒而言，計算中將片狀鋁顆粒等效為相同質(zhì)量的球形顆粒，β 的表達式為：

式中：Sf為片狀鋁顆粒的實際表面積，Ss為等效球形顆粒的表面積，Tm為鋁的熔點溫度。之所以引入β 是由于片狀顆粒的表面積比等效球形顆粒的表面積大，熱傳導時吸熱較多；當鋁熔化后，鋁液滴呈球形，β=1。而對于球形鋁顆粒，式（13）中不用考慮β，或者β=1。

1.3 鋁粉反應(yīng)模型

鋁粉的反應(yīng)模型是本文的創(chuàng)新點。兩步反應(yīng)模型應(yīng)用較廣泛，即將爆轟中鋁粉的反應(yīng)過程分為點火前的誘導階段和點火后的燃燒階段。針對誘導階段，Levitas 等[14]考慮氧化鋁殼的保護作用，通過理論分析認為鋁顆粒達到其熔點溫度后即可點火；洪滔等[15]通過分析鋁顆粒的激波點火機制，也認為鋁顆粒在爆轟動力學流場中的點火溫度為其熔點溫度。因此，本文中也采用該點火判據(jù)，當鋁粉達到其熔點溫度變?yōu)橐簯B(tài)后開始燃燒反應(yīng)。對于燃燒階段，鋁粉化學反應(yīng)簡化為 4 Al+3O2→2Al2O3，反應(yīng)速率I采用擴散燃燒模型[16]：

式中：Lv為氧化鋁的蒸發(fā)潛熱。對于鋁粉/空氣混合物爆轟，假設(shè)氣體中存在4 種組分，分別為O2、N2、Al(g)和Al2O3(l)；對于鋁粉/氧氣混合物爆轟，不存在N2組分。

計算模型的主要參數(shù)為：鋁密度ρ2=2 700 kg/m3，鋁熔化溫度Tm=931 K，鋁比熱容cV2=905 J/(kg·K)，鋁的燃燒熱qAl=31.5 kJ/g，氧化鋁蒸發(fā)潛熱Lv=4.77 kJ/g。

上述控制方程(1)～(16)為含有源項的歐拉方程，為了較好地對爆轟流場中的強間斷問題進行數(shù)值模擬，本文中采用三維的時-空守恒元解元(space-time conservation element and solution element method,CE/SE)差分格式求解不含源項的歐拉方程，然后用四階Runge-Kutta 法求解其源項部分，開發(fā)了懸浮鋁粉塵兩相爆轟的三維數(shù)值模擬程序，并且基于消息傳遞接口(message passing interface, MPI)技術(shù)，實現(xiàn)了程序的并行化設(shè)計。

2 數(shù)值方法驗證

為了驗證開發(fā)的數(shù)值模擬程序的可靠性，分別采用1 和2 mm 兩種網(wǎng)格尺寸對激波管問題進行了數(shù)值模擬，并與理論結(jié)果[20]進行了對比。在一根尺寸為1.0 m×0.1 m×0.1 m 的激波管中充滿理想氣體，距左端0.5 m 處有一隔膜將氣體分為兩部分，t=0 時刻兩部分氣體的狀態(tài)為：

計算開始時刻隔膜突然打開，左右兩部分氣體開始混合。圖1 分別顯示了某一時刻數(shù)值計算得到的密度和壓力分布與理論值的對比。由圖1 可以看出，采用這兩種網(wǎng)格計算得到的結(jié)果與理論結(jié)果均吻合較好，表明采用1～2 mm 的網(wǎng)格均可以達到精度要求。計算結(jié)果表明，采用本文中開發(fā)的程序可以較好地捕捉模擬激波問題，從而驗證了該程序的可靠性。

圖1 數(shù)值模擬得到的密度和壓力分布與理論計算結(jié)果[20]的比較Fig. 1 Distributions of density and pressure by simulation and theory[20]

3 計算結(jié)果與分析

為了驗證改進的鋁粉反應(yīng)模型，分別對鋁粉/空氣混合物爆轟以及鋁粉/氧氣混合物爆轟進行數(shù)值模擬。

3.1 鋁粉/空氣兩相爆轟

首先對Tulis 等[21]開展的爆轟管內(nèi)片狀鋁粉/空氣混合物爆轟實驗進行數(shù)值模擬。實驗所用的爆轟管長度為5.487 m，內(nèi)徑為152 mm。鋁粉的質(zhì)量濃度σ2為0.3 kg/m3，電子顯微鏡顯示該片狀鋁粉的比表面積為3～4 m2/g，本文計算中將該鋁粉等效為直徑為3.4 μm 的球形顆粒，相應(yīng)地，β=5。

考慮到爆轟管的軸對稱特性，建立1/4 計算模型。圖2(a)顯示了計算區(qū)域的橫截面(yOz)，尺寸為76 mm×76 mm，在y和z兩個方向上的網(wǎng)格數(shù)量均為61，即網(wǎng)格尺寸為1.246 mm。在笛卡爾坐標系中對爆轟管的圓周作近似處理，由圖2 可以看出，爆轟管壁由水平和垂直的單元組成。計算域x方向長度為6 m，網(wǎng)格數(shù)量為4 800，即網(wǎng)格尺寸為1.25 mm。為了提高計算效率，管壁外的節(jié)點單元（見圖2(a)中灰色部分）不參與計算。基于MPI，采用120 個進程并行計算，每個進程計算50 mm 長的計算域，并且記錄爆轟管軸上每隔1 m 位置處的壓力變化。

圖2 計算域橫截面和局部區(qū)域初始壓力云圖Fig. 2 Transverse cross section of the computational domain and initial pressure contour of partial domain

爆轟管的左端和圓周均采用固壁條件，右端6 m 處開口，并且在面xOy和xOz上的邊界采用周期邊界方法。圖2(b)為初始時刻局部區(qū)域(x=0～300 mm)壓力云圖，其中左邊31 mm 是點火段。本文參考文獻[22]，點火條件模擬炸藥起爆后的高溫高速爆轟產(chǎn)物，屬于強點火條件。由于實驗中采用2.8 g RDX的點火源，并且模擬采用1/4 模型，所以計算中采用高溫高壓的0.7 g N2作為點火源，保證點火源質(zhì)量和能量均為實驗點火源的1/4。點火區(qū)域的初始參數(shù)為：φ1=1, ρ1=5 kg/m3,T1=1 500 K,ux=3 235 m/s,uy=uz=0。右邊藍色區(qū)域為靜止的懸浮鋁粉塵和空氣兩相混合區(qū)域，初始參數(shù)和實驗值相同：ρ1=1.21 kg/m3, σ2=0.3 kg/m3,T1=T2=298 K, |u|=0,|v|=0。

圖3 顯示了數(shù)值計算得到的爆轟管軸上1～6 m 各位置處壓力隨時間的變化，鋁粉在前導激波的作用下點燃并燃燒釋放能量，支持爆轟波壓力逐漸升高，并使其在5～6 m 處發(fā)展為穩(wěn)定爆轟波，這和實驗觀測現(xiàn)象一致。

圖3 爆轟管軸線不同位置處壓力歷史Fig. 3 Pressure histories at different locations on the tube axis

圖4 顯示了計算得到的t=3.83 ms 時刻爆轟波陣面附近氣體壓力、氣體溫度、氣體密度、氣體軸向速度、鋁粉體積分數(shù)以及鋁粉軸向速度的分布。由圖4 可以看出，此時爆轟波陣面到達5.84 m位置，在前導激波的作用下，氣體被迅速壓縮至高溫高壓狀態(tài)，氣體和鋁粉塵均被加速，鋁粉塵的聚集導致局部濃度升高。當前導激波過后，鋁粉并不會立刻燃燒釋放能量，而是在其點燃前有個誘導點火區(qū)域，因而前導激波后氣體的溫度相對后面鋁粉燃燒區(qū)域的氣體溫度較低，這一點也可以從氣體溫度云圖中看出。由于管壁的摩擦作用，可以看到兩相速度均存在徑向的梯度分布。

圖4 3.83 ms 時刻鋁粉/空氣爆轟波陣面附近氣體壓力、氣體溫度、氣體密度、氣體軸向速度、鋁粉體積分數(shù)以及鋁粉軸向速度的分布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

數(shù)值模擬得到的穩(wěn)定爆轟波速度為1 536 m/s，與實驗結(jié)果1 510 m/s 吻合較好。計算結(jié)果表明，本文改進的鋁粉反應(yīng)模型適用于鋁粉/空氣兩相爆轟的數(shù)值模擬。此外，根據(jù)爆轟波的C-J (Chapman-Jouguet)理論，D=uCJ+cCJ（uCJ和cCJ分別表示氣體速度和氣體聲速），爆轟波CJ 參數(shù)分別為：pCJ=1.47 MPa，ρCJ=1.81 kg/m3，uCJ=534 m/s，TCJ=3 899 K。根據(jù)數(shù)值模擬，CJ 面離爆轟波陣面1.58 m，并且CJ 面上剩余鋁粉的體積分數(shù)為初始值的0.7%。這說明兩相爆轟的反應(yīng)區(qū)遠大于氣相爆轟的反應(yīng)區(qū)。

圖5 為計算得到的t=3.83 ms 時刻爆轟管軸上兩相的軸向速度、氣體溫度以及氣相組分和鋁粉的質(zhì)量濃度分布。模擬結(jié)果表明，前導激波過后，氣體速度急劇升高，鋁顆粒在氣體拖拽力的作用下也提高速度。但是由于固體顆粒具有較大的慣性，所以鋁顆粒的速度下降較慢，并且在離波陣面較遠處鋁顆粒的速度甚至高于氣體速度。

圖5 3.83 ms 時刻爆轟管軸上兩相的軸向速度、氣體溫度以及氮氣、氧氣和鋁粉的質(zhì)量濃度分布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在高溫下分解反應(yīng)吸熱的緣故，該數(shù)值模擬結(jié)果符合Glassman 觀察到的實驗現(xiàn)象[18]。此外，氣體溫度到達最高點后，會隨著距爆轟波陣面距離的增加而有所降低，這是因為Al2O3的沸點隨著壓力的降低而降低[19]，而它的沸點溫度決定了分解吸熱反應(yīng)是否發(fā)生以限制氣體的最高溫度。

同樣地，氮氣、氧氣以及鋁粉塵的濃度在前導激波過后迅速升高，此后鋁粉和氧氣由于燃燒反應(yīng)快速消耗。根據(jù)數(shù)值模擬結(jié)果，3.83 ms 時刻5.5 m 位置處剩余鋁粉塵的濃度為初始值的12.5%。

3.2 新模型的改進效果

為了體現(xiàn)新模型的改進效果，本節(jié)中采用改進前不考慮氧化鋁分解反應(yīng)的原模型對3.1 節(jié)中的工況鋁粉/空氣兩相爆轟進行數(shù)值計算，并和上節(jié)的計算結(jié)果進行對比，分析模型改進后的數(shù)值仿真效果。

圖6 為采用兩個模型計算得到的相同時刻爆轟波附近氣體溫度和氣體壓力的分布。由圖6 可以看出，采用改進前的模型計算得到的爆轟波后的氣體溫度高達4 800 K，該溫度明顯高于Al2O3的沸點?？紤]到Al2O3不能以氣體形式存在，所以在鋁粉的燃燒模型中加入Al2O3的高溫分解反應(yīng)是合理并且必要的。此外，改進前的模型由于無法限制鋁粉的燃燒放熱，所以采用其計算得到的爆轟速度和爆轟壓力均相對較高。

圖6 模型改進前后計算得到的爆轟波附近氣體溫度和氣體壓力的分布Fig. 6 Distributions of gas temperature and gas pressure around the steady detonation wave calculated by the improved and previous models

基于以上對比分析可以看出，本文改進后的模型可以反映懸浮鋁粉塵爆轟過程中存在的Al2O3高溫分解現(xiàn)象，能夠較準確地模擬計算爆轟流場。

3.3 鋁粉/氧氣兩相爆轟

為了檢驗改進的鋁粉反應(yīng)模型對不同氧化性氣氛的適用性，對鋁粉/氧氣混合物的兩相爆轟開展了數(shù)值模擬。混合物由濃度為700 g/m3、直徑為2.5 μm 的球形鋁粉和純氧氣組成[23]，爆轟在直徑為0.12 m的管道中進行。

數(shù)值計算中，仍采用1/4 管道模型作為計算域，其幾何形狀和圖2 相似，不同之處是該計算模型的橫截面為60 mm×60 mm，長度為10 m。在y和z方向上分為60 個網(wǎng)格，在x方向上有10 000 個網(wǎng)格，即3 個維度的網(wǎng)格尺寸均為1 mm。邊界條件和上文鋁粉/空氣兩相爆轟的計算模型相似。為了提高計算效率，數(shù)值模擬在250 個進程上并行計算，每個進程計算0.04 m 長的計算域。計算域左端的0.15 m 長為點火區(qū)域，鋁粉/氧氣混合物的參數(shù)為[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 長的局部區(qū)域氣體壓力、氣體溫度、氣體密度、氣體軸向速度、鋁粉體積分數(shù)以及鋁粉軸向速度的分布。由圖7 可以看出，此時爆轟波陣面到達7.27 m 位置，鋁氧爆轟的前導激波后的氣體溫度很高，表明鋁粉已開始燃燒，說明鋁粉的誘導點火區(qū)域很短，緊跟著前導激波陣面即開始燃燒反應(yīng)。這與上文的鋁粉/空氣兩相爆轟有所不同，原因是該工況下的鋁顆粒較小，表面質(zhì)量比相對較大，所以能夠更快地吸收熱量達到點火溫度，從而開始燃燒放熱。從鋁粉體積分數(shù)的云圖還可看出，波陣面后約0.03 m 之外鋁粉幾乎燃燒完畢，由于鋁顆粒較小以及純氧環(huán)境，該爆轟反應(yīng)區(qū)長度遠小于上節(jié)中鋁粉/空氣爆轟約1.60 m 的反應(yīng)區(qū)長度。

圖7 4.43 ms 時刻鋁粉/氧氣爆轟波陣面附近氣體壓力、氣體溫度、氣體密度、氣體軸向速度、鋁粉體積分數(shù)以及鋁粉軸向速度的分布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

數(shù)值模擬得到的該兩相爆轟波速度為1 610 m/s，比文獻值1 702 m/s[23]小5.4%，這可能是由于本文的數(shù)值模擬中考慮了管壁摩擦，而文獻中沒有考慮該摩擦作用。根據(jù)計算結(jié)果，pCJ=2.07 MPa, ρCJ=2.63 kg/m3,uCJ=543 m/s,TCJ=3838 K，CJ 面離爆轟波陣面0.93 m，并且CJ 面上沒有剩余鋁粉，和鋁粉/空氣兩相爆轟相比，該爆轟的CJ 反應(yīng)區(qū)更短。

圖8 為計算得到的t=4.43 ms 時刻爆轟管軸上兩相的軸向速度、氣體溫度以及鋁粉、氧氣的質(zhì)量濃度分布。由圖8 可以看出，前導激波過后，氣體和鋁顆粒相繼加速，在約7.18 m 位置處鋁顆粒燃燒完畢，表明鋁粉反應(yīng)區(qū)約0.09 m 長。由氣體溫度分布可以發(fā)現(xiàn)，在極小的區(qū)域內(nèi)氣體溫度即上升至最高點，這一現(xiàn)象和上文的鋁粉/空氣爆轟的情況不同，這是較小的鋁顆粒在更短的時間被點燃并燃燒放熱的緣故。

圖8 4.43 ms 時刻爆轟管軸上兩相的軸向速度、氣體溫度以及鋁粉、氧氣的質(zhì)量濃度分布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 結(jié) 語

鋁粉反應(yīng)模型是懸浮鋁粉塵兩相爆轟數(shù)值模擬研究的難點。本文中，通過引入鋁粉燃燒產(chǎn)物氧化鋁（Al2O3）在高溫下的分解反應(yīng)，改進了鋁粉的燃燒反應(yīng)模型，并且新的模型更符合物理現(xiàn)象。分別對鋁粉/空氣混合物爆轟以及鋁粉/氧氣混合物爆轟進行數(shù)值模擬，計算的結(jié)果和實驗結(jié)果、文獻值均吻合較好，表明改進的鋁粉反應(yīng)模型適用于不同氧化性氣氛中懸浮鋁粉塵爆轟的數(shù)值模擬。