固體炸藥爆轟的一種考慮熱學非平衡的反應流動模型?

李詩堯 于明

1)(中國工程物理研究院研究生院,北京 100088)2)(北京應用物理與計算數學研究所,北京 100094)(2017年11月21日收到;2018年8月20日收到修改稿)

基于固體炸藥爆轟過程中化學反應混合區內的固相反應物與氣相生成物處于力學平衡狀態及熱學非平衡狀態的事實,提出一種考慮熱學非平衡效應的反應流動模型來描述固體炸藥的爆轟流動現象.該爆轟流動模型的主要特點是,在反應混合物Euler方程和固相反應物質量守恒方程的基礎上,通過附加一套關于固相反應物的組分物理量的流動控制方程來表達固相反應物與氣相生成物之間的熱學非平衡效應.根據反應混合區內固相反應物與氣相生成物這兩種化學組分保持各自內能守恒的混合規則,并借助它們具有壓力相等的性質以及滿足體積分數總和為1的條件,推導獲得的附加方程有:固相反應物的內能演化方程、體積分數演化方程及反應混合物的壓力演化方程.這樣,建立的爆轟模型包括:反應混合物的質量守恒方程、動量守恒方程、總能量守恒方程、壓力演化方程,以及固相反應物的質量守恒方程、內能演化方程、體積分數演化方程.對所獲得的爆轟模型方程組采用一個時空二階精度的有限體積法進行數值求解,典型爆轟問題算例結果表明本文提出的固體炸藥爆轟模型是合理的.

1 引 言

固體炸藥的爆轟是劇烈的化學反應與復雜的激波流體力學高度耦合的非線性動力學現象,其宏觀理論模型從粗略到精密不斷深入發展著.早在19世紀末及20世紀初,英國工程師Chapman和法國力學家Jouguet各自獨立地提出一個簡化的理論模型,認為爆轟過程中化學反應是在一個無限薄的間斷面上瞬時完成的,化學反應過程無須加以考慮,其作用如同一個外加的能源反映到流體力學方程中來,流體力學的質量、動量和能量守恒方程通過跨過爆轟波陣面的間斷面建立,這就是Chapman-Jouguet(CJ)模型[1,2].后來在20世紀40年代至50年代,俄國物理學家Zeldovich、美國數學家von Neumann以及德國力學家D?ring分別建立了考慮化學反應過程的理論模型,認為爆轟化學反應具有有限反應速率,爆轟波由惰性的前導激波以及緊隨其后的具有一定寬度的化學反應區構成,化學反應區內物質看作是由固相反應物與氣相生成物這兩種化學組分組成的反應混合物,并且假設固相反應物與氣相生成物處于局域熱力學平衡狀態,即固相反應物與氣相生成物具有相同的速度、壓力以及溫度,這就是Zeldovich-von Neumann-D?ring(ZND)模型[3,4].直到目前,CJ模型和ZND模型仍是固體炸藥工程中應用最為廣泛的兩個爆轟理論模型[5,6].

CJ模型的缺點是沒有考慮化學反應過程;ZND模型盡管考慮了化學反應過程,但它使用的固體炸藥爆轟化學反應混合區內的固相反應物與氣相生成物處于局域熱力學平衡狀態的假設被認為是不合理的.眾所周知,ZND模型本身來源于氣體爆轟動力學學科,在氣體爆轟過程中反應物與生成物均為氣體物質,由于氣體分子的快速而密集的碰撞運動使得氣體爆轟化學反應混合區內的反應物與生成物能夠在極其短暫的時間內達到熱力學平衡狀態,因此氣體爆轟過程中假設反應物與生成物處于局域熱力學平衡狀態是合理的[7,8].然而,固體炸藥爆轟現象具有與氣體爆轟現象截然不同的動力學圖像.固體炸藥爆轟過程中,未反應固體炸藥在激波沖擊作用下起爆時,首先是炸藥晶體形成高溫高壓的“熱點”進而分解為細觀尺度的顆粒,然后炸藥顆粒表面快速燃燒使得固相反應物轉變成為氣相生成物[9?11],這樣爆轟化學反應混合區包含的固相反應物與氣相生成物這兩種化學組分物質可以看作以接觸間斷形式進行相互作用,于是固相反應物與氣相生成物之間能夠達到力學平衡狀態,即化學反應過程中固-氣兩相組分物質具有相同的速度和壓力,同時又由于固相反應物與氣相生成物的材料物性和微觀特性差異極大,氣體分子與固體顆粒之間往往需要經過較多次以及較長時間的碰撞作用才能達到熱學平衡狀態,即化學反應過程中固-氣兩相物質會具有各自不同的溫度或內能,因此,固體炸藥爆轟過程中固相反應物與氣相生成物應當處于力學平衡狀態及熱學非平衡狀態.這一狀況已經得到熱力學特征時間尺度估算[12,13]、細觀沖擊起爆數值模擬[14?16]以及反應過程實驗觀測[17?19]的證實.

值得指出的是,文獻[13]在估算固體炸藥爆轟過程的熱力學特征時間尺度的基礎上推導出一種考慮熱學非平衡的爆轟反應流動模型,本文借鑒了這種熱學非平衡思想,并且在如下三方面做了改進:1)文獻[13]組分體積分數通過物質等熵混合假設得出,本文通過物質內能演化得出,更符合物理實際情況;2)文獻[13]對反應混合物定義了混合狀態方程來封閉反應流動控制方程組,本文將混合物壓力直接作為一個自變量而無須定義混合物狀態方程,這樣的處理方式有更好的適應性;3)正是定義了混合物狀態方程使得進入混合區的激波往往產生附加的體積壓縮,為克服這一問題文獻[13]在數值求解過程中進行了壓力松弛和重新初始化處理,本文沒有定義混合物狀態方程因而無須進行壓力松弛和重新初始化處理,從而減小了誤差.

由此,基于固體炸藥爆轟過程中化學反應混合區內的固相反應物與氣相生成物處于力學平衡狀態及熱學非平衡狀態的事實,本文提出一種考慮熱學非平衡效應的反應流動模型來描述固體炸藥的爆轟流動現象.建立的爆轟模型包括:化學反應混合物的質量守恒方程、動量守恒方程、總能量守恒方程、壓力演化方程,以及固相反應物的質量守恒方程、內能演化方程、體積分數演化方程.典型爆轟問題算例結果表明本文提出的固體炸藥爆轟模型是合理的.

2 固體炸藥爆轟模型方程的推導

首先給出物理量的定義.考慮一個控制體包含有固相反應物與氣相生成物兩種組分物質,固相反應物的物理量標記下標“s”,氣相反應物的物理量標記下標“g”,則有組分密度ρs和ρg、組分內能es和eg、組分壓力ps和pg、組分速度us和ug、組分體積分數fs和fg、組分溫度Ts和Tg、組分比容τs和τg,且有fs+fg=1.對該控制體內的混合物密度ρ、混合物內能e、混合物比容τ有如下定義:

如果控制體內固相反應物與氣相生成物處于力學平衡狀態,則對該控制體的混合物壓力p及速度u有如下關系式:

不考慮各類耗散因素及外力做功情況,則在力學平衡狀態下固體炸藥爆轟反應流動方程組可用如下反應混合物的Euler方程與固相反應物的質量守恒方程來表達[20,21]:

(3a)式為反應混合物的質量守恒方程,(3b)式為反應混合物的動量守恒方程,(3c)式為反應混合物的總能量守恒方程,(3d)式為固相反應物的質量守恒方程.

對于(3a)—(3d)式,如果再補充反應混合區內的固相反應物與氣相生成物按照等溫等壓(熱學平衡及力學平衡)的規則進行混合,即有以下各式:

則(3a)—(3d)式滿足偏微分方程組封閉條件因而能夠得到定解,其中考慮了固體炸藥發生化學反應過程中單位質量的固相反應物轉變成為氣相生成物時釋放了數量為q的熱能.(3a)—(3d)式連同(4)式一起就構成了前述的ZND模型.

顯然,在熱學非平衡狀態下,固相反應物與氣相生成物之間的等溫條件已經不再滿足,必須尋找另外的混合規則以使得(3a)—(3d)式能夠滿足偏微分方程組封閉條件從而得到定解.根據爆轟反應流動的物理機理,可以認為反應混合區內的固相反應物與氣相生成物按照各自化學組分的內能保持守恒這一規則進行混合,以此為基礎可以給出與固相反應物相關的組分物理量的流動控制方程.為此,由(3a),(3b)和(3c)式可以把反應混合物的總能量守恒方程轉化為反應混合物的物質導數形式的內能方程:

將(1)式代入到(5)式,有

(6)式也可化為

(7)式兩個括號中式子分別表示固相反應物和氣相生成物的內能變化,再考慮到發生化學反應時固相反應物轉變成為氣相生成物會釋放熱量,以及由于固相反應物與氣相生成物之間存在溫度差而導致的熱量擴散,可設單位質量的固相反應物向氣相生成物進行了總數為q的熱量交換,則可以將(7)式分裂為固相反應物和氣相生成物的內能演化方程:

(8a)和(8b)式說明,盡管固相反應物和氣相生成物進行了熱量交換,但它們的總能量是保持不變的.

針對固相反應物流動過程,(8a)式借助混合物質量守恒方程可以化成如下形式:

(9)式中出現了固相反應物的體積分數變化這一變量.為了求出固相反應物體積分數的控制方程,需借助化學反應混合區內固相反應物與氣相生成物具有壓力相等的性質以及滿足體積分數總和為1的條件.

為簡化推導過程可不失一般性地認為,固相反應物具有ps(ρs,es)形式的狀態方程,氣相生成物具有pg(ρg,eg)形式的狀態方程,則在固相反應物與氣相生成物壓力相等即

的條件下,可以得到如下物質導數形式的固相反應物的壓力控制方程:

由(3d)式可得

由(3d)和(9)式可得

將(11)和(12)式代入(10)式,可得

(13)式也可以轉化成物質導數形式的固相反應物的體積分數控制方程:

采用與(14)式相同的思路,可以得到物質導數形式的氣相生成物的體積分數控制方程:

如果定義化學組分等熵聲速:

則(14)和(15)式變成

在固相反應物與氣相生成物滿足體積分數總和為1的條件下,即有fs+fg=1,則由(16)和(17)式可以獲得反應混合物的壓力演化方程:

將(18)式代入(16)式,可以獲得固相反應物的體積分數演化方程:

將(19)式代入到(9)式中,可以獲得固相反應物的內能守恒方程:

綜上推導過程,本文所獲得的考慮熱學非平衡效應的固體炸藥爆轟模型的反應流動控制方程組由(3a)—(3d)式、(18)—(20)式組成,具體包括:反應混合物的質量守恒方程、動量守恒方程、總能量守恒方程、壓力演化方程,以及固相反應物的質量守恒方程、內能演化方程、體積分數演化方程.

3 爆轟模型方程的數值方法與典型爆轟問題算例

3.1 爆轟模型方程的數值方法

為了便于數值求解,獲得的固體炸藥爆轟反應流動控制方程組(3a)—(3d),(18)—(20)式在二維直角坐標系下可以寫成如下矢量-矩陣形式的偏微分方程組:

式中狀態變量q=[ρ,ρu,ρv,ρE,ρsfs,ρsesfs,fs,p]T;矩陣A(q)與矩陣B(q)是具有類似表達形式的Jacobi矩陣,其中

源項s(q)=[0,0,0,0,?ρr,?ρr(q?p?),?ρr?,?ρrφ]T.這里,u,v為速度分量;

Jacobi矩陣A(q)的8個特征值和特征向量如下.

右特征向量為

相應的左特征向量為

Strang分裂方法[22]被用來求解這個帶源項的非守恒型雙曲律方程組(21),求解的方程組分裂為雙曲律方程組求解步與常微分方程組求解步.在常微分方程組求解步中,根據固體炸藥爆轟通常具有快反應過程和慢反應過程的特點[9?11],可以將其化學反應率分解成快反應率和慢反應率兩部分,即有r=rF+rS,相應地(21)式的源項也分解成與快反應率和慢反應率相關的兩部分,由此(21)式可以分裂為

其中,

在雙曲律方程組求解步(22a)中,一個具有波傳播特征的二階精度Godunov型格式被采用[23],該格式能夠很好地處理雙曲律方程組中的非守恒型對流項:

這里A±=RΛ±L,L和R是Jacobi矩陣A的左右特征向量構成的矩陣,Λ±是A的非負與非正特征值構成的對角矩陣,且?qi,j=qi,j?qi?1,j,?α=φ(θ)L?q,其中φ(θ)是限制函數;?Gi,j為橫向通量修正.B±有與之類似的性質.

在常微分方程組求解步(22b)中,一個具有顯-隱離散特性的二階精度Runge-Kutta格式被采用[24],顯式離散被用來處理非剛性源項,隱式離散被用來處理剛性源項,該格式能夠很好地處理爆轟化學反應過程中出現的快慢反應差異:

3.2 典型爆轟問題算例

對固體炸藥PBX9404的爆轟流動情況進行數值考察.該炸藥的主要特征參數如下(本節均采用cm-g-μs單位制):ρ0=1.842,pVN=0.563,uVN=0.347,ρVN=3.041,DCJ=0.880,pCJ=0.370,uCJ=0.229,ρCJ=2.488,q=0.102.炸藥爆轟化學反應率選用Lee-Tarver模型[25]:r=I(ρs/ρ?1?a)mλx+Gλy(1?λ)zpn,其中λ=ρsfs/ρ為炸藥爆轟固相反應物的質量分數,這里I,G,a,m,n,x,y,z均為與炸藥相關的常數;固相反應物和氣相生成物選擇應用廣泛的Jones-Wilkins-Lee形式狀態方程[26]:

這里Ak,Bk,R1k,R2k,ωk均為與炸藥相關的常數.

3.2.1 平面一維爆轟問題

考慮長度為4.0的一維炸藥爆轟起爆與爆轟傳播問題,左端以CJ條件起爆,在如下幾種網格分點條件下:5 cells/mm,10 cells/mm,30 cells/mm,50 cells/mm,100 cells/mm,獲得爆轟波傳播過程中幾個典型時刻反應混合物的壓力、密度、速度分布以及固相反應物的質量分數分布,如圖1—圖3所示(文中僅展示5 cells/mm,30 cells/mm,100 cells/mm三種情況),圖中每條曲線對應的時刻為t=0.4,0.8,1.2,1.6,2.0,2.4,2.8,3.2,3.6,4.0.從計算結果可以看出:

1)隨著網格變細,計算獲得的壓力von Neumann尖峰值變大,爆轟傳播速度變大,達到定常狀態的時間變短,化學反應區變窄;這五種網格分點條件下的壓力尖峰值分別約為0.461,0.503,0.543,0.555,0.561;爆轟傳播速度分別約為0.810,0.840,0.865,0.870,0.875;定常時間分別約為1.2,0.8,0.5,0.3,0.2;化學反應區寬度分別約為0.0175,0.0115,0.0091,0.0087,0.0085;爆轟壓力尖峰值、傳播速度、定常時間、化學反應區寬度隨網格尺度的變化如圖4所示;

圖1 5 cells/mm條件下爆轟波傳播過程物理量變化趨勢 (a)壓力;(b)速度;(c)密度;(d)固相反應物質量分數Fig.1.Physical variables of detonation wave under 5 cells/mm:(a)Pressure;(b)velocity;(c)density;(d)mass fraction.

圖2 30 cells/mm條件下爆轟波傳播過程物理量變化趨勢 (a)壓力;(b)速度;(c)密度;(d)固相反應物質量分數Fig.2.Physical variables of detonation wave under 30 cells/mm:(a)Pressure;(b)velocity;(c)density;(d)mass fraction.

圖3 100 cells/mm條件下爆轟波傳播過程物理量變化趨勢 (a)壓力;(b)速度;(c)密度;(d)固相反應物質量分數Fig.3.Physical variables of detonation wave under 100 cells/mm:(a)Pressure;(b)velocity;(c)density;(d)mass fraction.

圖4 不同網格尺度下爆轟波傳播過程物理量變化 (a)尖峰壓力;(b)傳播速度;(c)定常時間;(d)反應區寬度Fig.4.Physical variables of detonation wave under different mesh scales:(a)Peak pressure;(b)propagation velocity;(c)stationary time;(d)width of reaction area.

圖5 柱面散心爆轟波爆轟波傳播過程物理量變化趨勢 (a)壓力;(b)徑向速度;(c)密度;(d)固相反應物質量分數Fig.5.Physical variables of cylindrically divergent detonation wave:(a)Pressure;(b)radial velocity;(c)density;(d)mass fraction.

2)30 cells/mm分點條件下von Neumann壓力值和爆轟傳播速度已接近理論值.

3.2.2 柱面一維散心爆轟問題

考慮半徑為4.0的圓柱狀炸藥爆轟起爆與爆轟傳播問題,中心線以CJ條件起爆,遵守上節得出的30 cells/mm網格分點條件.計算獲得爆轟波傳播過程中幾個典型時刻反應混合物的壓力、密度、徑向速度分布以及固相反應物的質量分數分布,如圖5所示,圖中每條曲線對應的時刻為t=0.05,0.1,0.2,0.4,0.8,1.2,1.6,2.0,2.4,2.8,3.2,3.6,4.0.從計算結果可以看出:

1)炸藥爆轟波傳播過程中爆速和爆壓均隨波面半徑的增加而增加,但壓力傳播3.0μs后,速度和密度傳播2.0μs后增加不顯著;

2)散心爆轟傳播過程中最大壓力值約為0.527,低于平面爆轟最大壓力值;反應區寬度約為0.013,比平面爆轟反應區寬50%;這是由于散心爆轟波在傳播過程中的面積效應使得化學反應變慢,因而傳播壓力更低,化學反應區更寬.

3.2.3 柱面散心爆轟的反射問題

考察無限厚炸藥表面上有相距2l長度的兩線狀起爆器A和B同時起爆,如圖6所示,產生兩列相同的柱面散心爆轟波a,這兩列散心爆轟波傳播到對稱面Ox處產生碰撞爆轟波b,接著繼續向前傳播形成正規反射波c,之后進一步向前傳播形成爆轟波Mach反射d.由于在相互碰撞作用點處壓力相對增高,所以在對稱面附近會引起傳播速度增加,初始落后的Mach桿會逐漸追趕上散心入射波,最后兩散心爆轟波陣面趨于一致而形成平面爆轟波.

圖6 柱面散心爆轟相互作用Fig.6.Interaction of two cylindrically divergent detonation waves.

圖7 柱面散心爆轟波反射壓力圖Fig.7.Pressure contours of two cylindrically divergent detonation waves.

采用文獻[27]中的算例,兩起爆點距離為5.08 cm,計算域取下半平面并設四邊為固壁條件,仍采用30 cells/mm網格分點條件,得到不同時刻爆轟流場的壓力等值線(如圖7(a)—(f)所示).其中圖7(a)為炸藥起爆形成擬定常散心爆轟波,圖7(b)為兩散心爆轟波在對稱面處發生對碰,圖7(c)為兩爆轟波在碰撞角度較小條件下發生的正規反射,圖7(d)為爆轟波的正規反射向Mach反射過渡狀態,圖7(e)為兩爆轟波在碰撞角度較大條件下發生的Mach反射,圖7(f)顯示整個爆轟波陣面趨于一致而形成平面爆轟波.該數值模擬結果與文獻[27]結果符合.

當兩爆轟波碰撞作用的夾角大于100?時,出現Mach反射現象,在爆轟波相互作用點附近形成三波點A(圖8),其中AB為進行入射的爆轟波陣面,AC為Mach桿,AD為反射的激波.Mach桿后的流場壓力比入射爆轟波后的流場壓力高約25%—30%,這種高壓狀態驅動爆轟波使得Mach桿運動速度比散心爆轟波運動速度快約11%—13%,因此Mach桿會逐漸趕上散心爆轟波.

圖8 柱面散心爆轟波Mach反射壓力圖Fig.8.Mach reflected pressure for two cylindrically divergent detonation waves.

4 結 論

本文提出一種考慮熱學非平衡效應的反應流動模型來描述固體炸藥的爆轟流動現象,該爆轟模型通過附加固相反應物的內能演化方程、體積組分演化方程以及化學反應混合物的壓力演化方程來表達固相反應物與氣相生成物之間的熱學非平衡效應.從固體炸藥的平面一維爆轟、柱面一維散心爆轟以及柱面散心爆轟波反射等幾個典型算例來看,該模型能夠正確預測固體炸藥爆轟的起爆及傳播過程中的關鍵特征,如von Neumann壓力、爆轟波傳播速度、化學反應區寬度、爆轟波反射等,數值模擬結果證明本文提出的固體炸藥爆轟的理論模型是合理的.