昝文濤, 洪 滔, 董賀飛

(1. 北京理工大學(xué)機電學(xué)院, 北京 100081; 2. 北京應(yīng)用物理與計算數(shù)學(xué)研究所, 北京 100094)

1 引 言

粉塵爆炸是工業(yè)生產(chǎn)中常見的事故,2009～2013年我國共發(fā)生金屬粉塵引發(fā)的爆炸共有20起,其中大部分為鋁鎂粉塵爆炸引起[1]。近年軍事領(lǐng)域中溫壓武器發(fā)展迅速,其中鋁粉等固體燃料由于其質(zhì)輕含能高的優(yōu)點得到越來越廣泛的應(yīng)用,因此研究鋁粉塵爆轟在工業(yè)安全和軍事領(lǐng)域均具有重要意義。

Pawel Kosinski[2]研究了容器中產(chǎn)生的爆轟波通過管道傳播進入另一充滿粉塵的空間內(nèi)的數(shù)學(xué)模型。洪滔[3-4]研究了鋁顆粒激波點火機制,模擬了爆轟波管中的鋁粉塵爆轟。劉慶明[5]實驗研究了管道中鋁粉塵在弱起爆條件下的燃燒轉(zhuǎn)爆轟(DDT)過程。韋偉[6]基于CE/SE方法模擬二維管道內(nèi)鋁粉塵爆轟過程。滕宏輝[7]研究了真實比熱模型下鋁粉塵的兩相爆轟波。李鑫等[8]實驗研究了不同類型微/納米鋁粉的點火燃燒特性,對比了不同粒徑以及不同介質(zhì)包覆后的鋁粉點火延遲時間。

由于條件限制,上述研究大多以爆轟波管實驗為主,數(shù)值模擬也大多集中在簡單的二維平面管道方面,而關(guān)于鋁粉塵爆轟波在較復(fù)雜的幾何空間內(nèi)傳播研究較少,無論是工業(yè)生產(chǎn)車間,還是巷道等防御工事,其中存在由通風(fēng)管道或巷道連接的幾何空間,在其中的一處發(fā)生爆轟后容易通過管道等發(fā)生連鎖爆炸,造成大范圍毀傷破壞,對通過管道連接的空間內(nèi)鋁粉塵的爆轟研究在粉塵爆轟傳播及毀傷效應(yīng)及防護工程中有著重要的實際意義。

為此,本研究采用CE/SE(時-空守恒元解元)算法[9-11]求解兩相流模型方程組,對空氣中懸浮鋁粉塵在管道連接空間內(nèi)的爆轟進行了數(shù)值模擬,對鋁粉塵爆轟波的傳播過程及爆轟波的繞射等物理過程進行了研究。

2 兩相流方程及反應(yīng)模型

模擬采用了二維兩相流模型[12],假設(shè)鋁顆粒為球形,初始半徑都相同,單個顆粒溫度均勻,氣體中各組分都是均勻混合的,忽略了粒子間的相互作用,不考慮顆粒與氣體間的輻射作用,忽略了粒子和固壁間的熱傳導(dǎo)。

氣相方程:

(1)

(2)

(3)

IqAl-Fxu2-Fyv2

(4)

固相方程:

(5)

(6)

(7)

(8)

(9)

組分方程:

(10)

氣體狀態(tài)方程:

(11)

式中,角標(biāo)1、2分別代表氣體和Al。變量ρ為密度,kg·m-3;u為橫向速度,m·s-1;v為縱向速度,m·s-1;p為壓力,Pa;e為內(nèi)能,J;φ(φ1+φ2=1)為體積分?jǐn)?shù),n為單位體積內(nèi)鋁顆粒粒子數(shù);T為溫度,K;y為組分濃度,w為分子量。

源項中I為單位體積內(nèi)鋁顆粒的質(zhì)量變化率,qAl為單位質(zhì)量的鋁顆粒的反應(yīng)能,J。τ為氣體中各組分質(zhì)量生成率[13],當(dāng)溫度達(dá)到或超過Al2O3沸點時,Al2O3會發(fā)生分解[14],溫度將保持在Al2O3沸點。

對于鋁顆粒,采用以下反應(yīng)模型[15]:

(12)

(13)

式中，Tign為鋁粉塵點火溫度[2];r2為粒子半徑,m;d0為粒子初始直徑,m;ψ為氣體中氧氣的摩爾份數(shù);Fx為氣體對粒子在x方向的拖曳力,N;Fy為氣體對粒子在y方向的拖曳力,N。

(14)

(15)

拖曳系數(shù):

(16)

Re為雷諾數(shù),Q為氣體與粒子間的熱傳導(dǎo):

(17)

式中，λ=0.1 J·(m·s·K)-1為氣體導(dǎo)熱系數(shù)[16],Nu為Nusselt數(shù)[12],Nu=2+0.459Re0.55Pr0.33，Pr是普朗特數(shù)。

通過CE/SE方法求解歐拉方程組,采用四階龍格庫塔方法求解方程組的源項,編制程序模擬爆轟波在二維空間的發(fā)展過程。

3 數(shù)值模擬結(jié)果

3.1 程序驗證

為檢驗程序合理性,模擬了激波管問題,氣體密度、壓力和速度等進行了無量綱化,初始條件:t=0時刻,

對于理想氣體,絕熱指數(shù)1.4,網(wǎng)格數(shù)量200,計算結(jié)果如圖1所示。

圖1t=0.2時刻激波管的壓力分布

Fig.1 The pressure distribution of shock tube att= 0.2

由圖1可看出,網(wǎng)格數(shù)量200時的數(shù)值解與理論解符合較好,可以很好地模擬激波管問題。

對鋁粉塵與空氣當(dāng)量比為1時在管道內(nèi)的爆轟波傳播過程進行了模擬驗證,鋁粉塵濃度為0.304 kg·m-3,鋁顆粒半徑為1.7 μm,爆轟波管直徑為15.2 cm。數(shù)值模擬得到的爆轟波參數(shù)為[17]:D=1.63 km·s-1,ρCJ=2.43 kg·m-3,uCJ=673 m·s-1,pCJ=2.04 MPa,T=3800 K,p=3.31 MPa,與Tulis[18]等由實驗中得到的鋁粉塵的爆速為1650 m·s-1的結(jié)果符合良好。

3.2 計算條件

計算了管道連接空間的二維模擬模型,區(qū)域尺寸如圖2所示,左側(cè)區(qū)域為3 m×3 m空間,除管道出口外均為固壁; 中間管道直徑0.2 m,長2 m; 右側(cè)空間3 m×1 m,上下為固壁,右端開口。模擬網(wǎng)格尺寸為5 m×5 m,時間步長Δt=0.2×min(dxdy)/(5×(u1+c)),c為聲速。鋁粉塵的密度為0.304 kg·m-3,鋁顆粒半徑為2.0 μm。起爆區(qū)域位于左上角,起爆條件為φ1=1,ρ1=2.2 kg·m-3,u=1400 m·s-1,v=-1400 m·s-1,T1=3200 K,該起爆條件是為了符合實驗中起爆炸藥的質(zhì)量與能量而確定的。A-Q為數(shù)據(jù)采集點。

圖2 數(shù)值模擬區(qū)域尺寸模型

Fig.2 Area size model used for numerical simulation

3.3 數(shù)值模擬結(jié)果

3.3.1 左側(cè)空間中的爆轟波數(shù)值模擬結(jié)果

圖3是數(shù)值模擬得到的不同時刻的爆轟波流場壓力圖。由圖3可見,2.36 ms時刻(圖3a),由于起爆區(qū)域較小位于左上角,爆轟波陣面接近球形,此時爆轟波波陣面同時到達(dá)右壁面和下側(cè)壁面發(fā)生反射,在固壁處形成的高壓區(qū)壓力達(dá)到6.5 MPa; 2.53 ms時刻(圖3b),此時爆轟波反射區(qū)更明顯,爆轟波與空間右側(cè)固壁產(chǎn)生的1號反射波向左傳播,爆轟波與空間下固壁產(chǎn)生的2號反射波向上傳播,而未發(fā)生反射區(qū)域爆轟波壓力為2.6 MPa,速度為1414 m·s-1,由于傳播距離不夠長,此時未達(dá)到穩(wěn)定爆轟,爆轟波速度及壓力均低于表1中的結(jié)果; 3.16 ms時刻(圖3c),此時爆轟波到達(dá)右下角處,同時爆轟波與右固壁和下固壁產(chǎn)生的1號和2號反射波在此區(qū)域交匯形成局部高壓區(qū),壓力可達(dá)18 MPa; 3.90 ms時刻(圖3d),此時左側(cè)空間中主要存在四類區(qū)域,Ⅰ區(qū)為左上角處區(qū)域,此次空間為爆轟波后流場; Ⅱ區(qū)為右上角和左下角區(qū)域,這兩部分區(qū)域為反射波后流場,右上角區(qū)域是1號反射波后形成,左下角區(qū)域是2號反射波后形成的; Ⅲ區(qū)為對角線中間區(qū)域的局部高壓區(qū),右上角和左下角區(qū)域為1號和2號反射波此時相遇形成; Ⅳ區(qū)為位于右下角高壓區(qū),是由1號和2號反射波與固壁作用后形成的多次反射波疊加形成的復(fù)雜流場。

a. 2.36 ms

b. 2.53 ms

c. 3.16 ms

d. 3.90 ms

圖3 數(shù)值模擬得到的不同時刻的模擬區(qū)鋁粉塵爆轟波流場的壓力分布

Fig.3 Pressure distribution of the aluminum dust detonation wave flow field in simulation area at different time obained by numerical simulation

圖4為3.90 ms時刻爆轟波流場的溫度分布圖。圖4中顯示在左側(cè)空間內(nèi)波后流場為3400 K以上的高溫。

圖4t=3.90 ms時刻鋁粉塵爆轟波流場的溫度分布

Fig.4 Temperature distribution of the aluminum dust detonation wave flow field at 3.90 ms

圖5分別為空間中不同數(shù)據(jù)采集點處的壓力隨時間變化曲線。圖5a中六個點(A,B,C,D,E,F)均勻分布在左側(cè)空間的左上角與右下角斜對角線上,從圖5中可以看出在鋁粉塵起爆后生成的爆轟波壓力不斷增長,到達(dá)E(2.5,0.5)處時壓力達(dá)到2.25 MPa,當(dāng)爆轟波到達(dá)F(3.0,0.0)處時,爆轟波與1號、2號反射波相交形成三波點,此時壓力躍升到18 MPa,為空間區(qū)域壓力最高點。而后由1號反射波與下固壁作用產(chǎn)生的1-1號反射波和由2號反射波與右固壁作用產(chǎn)生的2-1號反射波同時到達(dá)E(2.5,0.5)點,此時雙波疊加壓力為4.9 MPa,隨著反射波逐漸傳播,壓力逐漸下降。

圖5b為距離左側(cè)空間上固壁處0.5 m處的一排點(A,G,H,I,J)的壓力隨時間變化曲線,在爆轟向右側(cè)固壁傳播,爆轟波波陣面壓力逐漸增大,到達(dá)J(2.5,2.5)處時壓力為2.3 MPa; 在爆轟波與壁面作用發(fā)生反射后會生成反射波,反射波經(jīng)過J(2.5,2.5)時壓力為2.1 MPa,而后反射波強度隨傳播距離增大逐漸降低; 由于1號和2號反射波同時到達(dá)A(0.5,2.5)點,因此壓力躍升至2.1 MPa,而后可以看出由于2號反射波經(jīng)過,G-J點壓力也會出現(xiàn)躍升并且逐漸降低。從圖中可以看出,A點由于反射波的作用,壓力還有一次明顯的躍升過程。

a. pressure at collecting points on the diagonal line at the left side of space

b. pressure at horizontal collection points at the left side of space

圖5 不同采集點壓力隨時間變化

Fig.5 Change in pressure with time at different collecting points

3.3.2 連接管道中爆轟波傳播的數(shù)值模擬結(jié)果

圖6為2.62 ms時刻管道入口附近的流場壓力圖和密度圖。右側(cè)的波為進入管道內(nèi)部傳播的初始爆轟波。此時爆轟波沿管道左上角繞射進入管道內(nèi)部,在管道左上角區(qū)域形成低壓低密度區(qū)域,在管道左下角區(qū)域由于爆轟波與固壁作用形成高壓高密度區(qū)域。3.90 ms時刻爆轟波到達(dá)管道出口處,爆轟波波后流場溫度從圖4中可以看出穩(wěn)定在3400 K。

圖7為沿管道中心線上均勻分布的5個點(K,L,M,N,O)和右側(cè)空間處2個點(P,Q)壓力隨時間變化,爆轟波到達(dá)管道入口時壓力為2.6 MPa,在爆轟波繞射進入管道后壓力略下降,而后爆轟波繼續(xù)增長,在到達(dá)管道出口處時爆轟波速度為1571 m·s-1,壓力達(dá)到2.85 MPa。

對比模擬了0.2 m寬管道中鋁粉塵的爆轟波傳播過程。表1為爆轟波達(dá)到穩(wěn)定時爆轟波參數(shù),模擬得到網(wǎng)格尺寸5 mm時爆速為1591 m·s-1,爆轟波陣面壓力3.01 MPa,CJ點氣體密度為2.07 kg·m-3,CJ點壓力為1.73 MPa,CJ點溫度是3794 K,CJ點氣體速度為510 m·s-1。

a. pressure nephogram b. density nephogram

圖6t=2.62 ms時刻管道入口處流場的壓力和密度分布

Fig.6 Pressure and density distributions of the flow field at the entrance of channel at 2.61 ms

圖7 管道中不同采集點壓力隨時間變化

Fig.7 Change in pressure with time at different collecting points in channel

表1 穩(wěn)定爆轟波參數(shù)

Table 1 Parameters for stable detonation wave

gridsize/mmwavevelocity/km·s-1pressure/MPadensityofCJpoint/kg·m-3pressureofCJpoint/MPatemperatureofCJpoint/KvelocityofCJpoint/m·s-1515702．912．041．713795494

與表1比較發(fā)現(xiàn),管道出口處爆轟波基本達(dá)到穩(wěn)定狀態(tài)。爆轟波進入右側(cè)空間中爆轟波壓力會下降。

3.3.3 右側(cè)空間中的爆轟波傳播的數(shù)值模擬結(jié)果

圖8為4.03 ms時刻管道出口附近流場壓力分布和密度分布。此時爆轟波離開管道傳播進入右側(cè)空間內(nèi),爆轟波沿管道上下兩個拐角繞射進入大空間內(nèi)傳播,在拐角處形成兩個低壓低密度區(qū)域,此時爆轟波由于繞射的效應(yīng),往上下傳播的波陣面較弱,因此波陣面后的壓力密度較低,往右側(cè)傳播的波陣面較強,波陣面后壓力密度相對較高。爆轟波陣面的壓力和密度沿管道中心軸線成對稱分布,呈蘑菇狀分布。圖9為4.37 ms時刻爆轟波流場的溫度分布,爆轟波從管道進入右側(cè)空間后,波后流場溫度為3600 K,由于爆轟波繞射,在管道出口處溫度比波后流場小,為3000 K。圖10為4.62 ms時刻右側(cè)空間內(nèi)流場壓力,此時右側(cè)空間內(nèi)爆轟波壓力為1.9 MPa,速度為1400 m·s-1,在計算區(qū)域內(nèi)未達(dá)到穩(wěn)定爆轟狀態(tài)。

a. pressure distributions b. density distributions

圖8t=4.03 ms時刻管道出口附近的壓力和密度分布

Fig.8 Pressure and density distributions near the exit of channel at 4.03 ms

圖9 4.37 ms時刻鋁粉塵爆轟波流場的溫度分布

Fig.9 Temperature distribution of the aluminum dust detonation wave flow field at 4.37 ms

圖10t=4.62 ms時右側(cè)空間內(nèi)流場壓力云圖

Fig.10 Pressure distributions of the flow field in the right side of space at 4.62 ms

4 結(jié) 論

采用兩相流模型數(shù)值模擬了懸浮空氣中的鋁粉塵濃度為0.304 kg·m-3,鋁顆粒半徑為2.0 μm時爆轟波在管道連接的兩個大空間內(nèi)起爆及爆轟波傳播過程。得到如下結(jié)論:

(1)左側(cè)空間中左上角鋁粉塵點火后形成爆轟波在空間內(nèi)傳播,當(dāng)爆轟波到達(dá)固壁處后產(chǎn)生反射波,兩個壁面產(chǎn)生的反射波相遇后相互作用在右下角區(qū)域交匯形成18 MPa的局部高壓區(qū),而后兩個反射波隨傳播會產(chǎn)生交匯形成相互作用區(qū)域,同時兩個反射波會繼續(xù)與固壁產(chǎn)生反射形成復(fù)雜流場區(qū)域。

(2)爆轟波能通過管道入口處繞射進入管道,在管道中傳播并在到達(dá)管道出口處時爆轟波速度為1571 m·s-1,壓力達(dá)到2.85 MPa,基本達(dá)到穩(wěn)定傳播狀態(tài)。

(3)在爆轟波繞射后傳播進入右側(cè)空間,在管道拐角處產(chǎn)生繞射形成對稱低壓低密度區(qū)域,爆轟波波后區(qū)域壓力密度分布呈蘑菇狀,在計算區(qū)域內(nèi)未達(dá)到穩(wěn)定爆轟。

(4)在整個計算區(qū)域內(nèi),爆轟波陣面后大部分區(qū)域的溫度為3400 K以上的高溫,僅在管道出口處由于爆轟波發(fā)生繞射溫度為3000 K。

上述模擬結(jié)果反映了鋁粉塵爆轟的高溫,易傳播的特點,以及爆轟波傳播過程產(chǎn)生的反射和繞射的效應(yīng),為復(fù)雜幾何形狀空間內(nèi)的鋁粉塵爆轟的傳播及毀傷效應(yīng)提供了認(rèn)識。

參考文獻(xiàn):

[1] 多英全,劉垚楠,胡馨升. 2009～2013年我國粉塵爆炸事故統(tǒng)計分析研究[J]. 中國安全生產(chǎn)科學(xué)技術(shù), 2015, 11(2): 186-190.

DUO Ying-quan, LIU Yao-nan, HU Xin-sheng. Statistical analysis on dust explosion accidents occurring in China during 2009-2013[J].JournalofSafetyScienceandTechnology, 2015, 11(2): 186-190.

[2] Pawel K, Alex C. Hoffmann. Dust explosions in connected vessels: mathematical modelling[J].PowderTechnology, 2005, 155: 108-116.

[3] 洪滔, 秦承森. 鋁顆粒激波點火機制初探[J]. 爆炸與沖擊, 2003, 23(4): 295-299.

HONG Tao, QIN Cheng-sen. Mechanism of shock wave ignition of aluminum particle[J].ExplosionandShockWaves, 2003, 23(4): 295-299.

[4] 洪滔, 秦承森. 爆轟波管中鋁粉塵爆轟的數(shù)值模擬[J]. 爆炸與沖擊, 2004, 24(3): 193-200.

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

[5] LIU Qing-ming, LI Xiao-dong, BAI Chun-hua. Deflagration to detonation transition in aluminum dust-air mixture under weak ignition condition[J].CombustionandFlame, 2009,156: 914-921.

[6] 韋偉, 翁春生. 基于CE/SE方法的鋁粉塵爆轟二維兩相數(shù)值計算[J]. 彈道學(xué)報, 2012, 24(4): 99-102.

WEI Wei, WENG Chun-sheng. Two-dimension two-phase-flow numerical simulation of aluminum-dust detonation based on CE/SE method[J].JournalofBallistics, 2012, 24(4): 99-102.

[7] 滕宏輝, 楊旸, 姜宗林. 真實比熱模型中鋁粉塵兩相爆轟波的數(shù)值研究[J]. 計算物理, 2013, 30(1): 44-52.

TENG Hong-hui, YANG Yang,JIANG Zong-lin. Realistic heat capacity effects in two phase aluminum dust detonations[J].ChineseJournalofComputationalPhysics, 2013, 30(1): 44-52.

[8] 李鑫, 趙鳳起, 郝海霞, 等. 不同類型微/納米鋁粉點火燃燒特性研究[J]. 兵工學(xué)報,2014, 35(5): 640-647.

LI Xin, ZHAO Feng-qi, HAO Hai-xia, et al. Research on ignition and combustion properties of different micro/nano-aluminum powders[J].ActaArmamentaril, 2014, 35(5): 640-647.

[9] 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].JournalofComputationalPhysics, 1995, 119(2): 295-324.

[10] Zhang D L, Wang J T, Wang G. High-order CE/SE method and applications[J].ChineseJournalofComputationalPhysics, 2009, 26(2): 211-220.

[11] Wang J T, Zhang D L, Liu K X. A eulerian approach based on CE/SE method for 2D multimaterial elastic-plastic flows[J].ChineseJournalofComputationalPhysics, 2007, 24(4): 395-401.

[12] 方丁酉. 兩相流體力學(xué)[M]. 國防科技大學(xué)出版社. 1988: 117-129.

[13] 洪滔, 秦承森, 林文洲. 懸浮RDX炸藥和鋁顆?；旌戏蹓m爆轟的數(shù)值模擬[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].ExplosionandShockWaves, 2009, 29(5): 468-473.

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

[15] Price E W. Combustion of metalized propellants[A].ProgressinAstronauticsandAeronautics:FundamenalsofSolid-PropellantCombustion, AIAA,New York, 1984, 90: 479-513.

[16] 楊世銘, 陶文銓. 傳熱學(xué)[M]. 高等教育出版社. 2006: 559.

[17] 昝文濤, 洪滔, 董賀飛. 基于CE/SE方法關(guān)于RDX-AL懸浮粉塵在空氣中的兩相爆轟的數(shù)值模擬[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].ExplosionandShockWaves, 2016, 36(5): 603-610.

[18] Tulis A J, Selman J R. Detonation tube studies of aluminum particles dispersed in air[C]∥19th International Symposium on Combustion.The Combustion Institute, 1982: 655-663.