中圖分類號：0381 國標(biāo)學(xué)科代碼：13035 文獻(xiàn)標(biāo)志碼：A

Numerical simulation for shock to detonation process of explosive nitromethane containing cavities

XIAO Min1，WANG Cheng2，YANG Tonghui3 （1. College ofScience， China University of Petroleum，Beijing102249，China; 2.State KeyLboratoryofExplosioncienceandafetyrotectioneiingInstituteofchnologyBeijingChina; 3.School ofAerospace Engineering， Tsinghua University，Beijing l0oo84， China）

Abstract：Inordertostudytheinitiationprocessofliquidexplosivenitromethanecontainingcavitiesunder shock wave loading，anEulerianmulti-materialcomputationalapproachbasedonthelevelsetmethodwasdeveloped.ThereactiveEuler equation wereadoptedasthegoverming equations，the levelsetmethod wasutilizedto track themulti-medium interface between the chemicalreactionmixture and thecavity.To improve the robustnessofcalculation method，the modified ghost fluidmethod wasapplied incomputationalcelsnearthe interface.Basedonthe modified ghostfluid method，amulti-medium problem was transformed intoasingle media problem.For these twofluid phasesonboth sidesofthe interface，the highorder weightedessentialnon-oscilatoryfinitedifferencemethodwasimplemented tocalculatethenumerical fluxesoncel boundary，making the simulationresultsreliable. However，the Jones-Wilkins-Leeequationof state difers greatly from the ideal gasequationofstate.Inadition，the massfractionofdetonation productdirectlyaffects the transformationprocess betweentheconserved variablesand the primitive variables inreaction zone，making itdificultto provideanexplicit expresion for the equationof state of explosive mixture.Inorder to solve the above problems，a ghost fluidstate prediction method based on the Harten-Lax-van Ler-contact （HLLC）approximate Riemann solver was developed.Bydealing with a complexmulti-medium Riemann problemconsidering chemicalreaction，thevariable statesof ghost fluidonboth sides of the interface can be obtained.The multi-medium calculation method was used tosimulate the interaction problems between liquid nitromethaneandthecavityundertheloadingcondition withdiferentimpactstrengths.Thenumericalresultsillustratetatthe method proposed inthepapercancapture theentirefluiddynamics processofcavitycompresson，cavitycolapse，cavity closure and cavity disappearance.

Keywords： liquid nitromethane; cavity colapse; multi-medium interface; level set method; ghost fluid

炸藥在澆鑄、壓裝過程中很難避免在內(nèi)部形成空穴。根據(jù)熱點(diǎn)理論，非均質(zhì)炸藥在沖擊波作用下，空穴受壓塌陷產(chǎn)生局部高溫高壓區(qū)域，當(dāng)熱點(diǎn)區(qū)域足夠大或者高溫持續(xù)時間足夠長時，會引起炸藥局部發(fā)生點(diǎn)火，實(shí)現(xiàn)爆轟，點(diǎn)火時間會明顯短于點(diǎn)燃純凈炸藥所需的時間。因此，考察炸藥對于空穴的敏感性具有重要的研究意義和應(yīng)用價值。

科研人員對沖擊波與不同排列下孔洞的相互作用問題進(jìn)行過一些數(shù)值模擬研究。Mader[2]模擬了液體炸藥硝基甲烷內(nèi)的空穴坍塌過程，這是一項(xiàng)早期研究成果，受到分辨率的限制。Lauer等[3]分別用Tait狀態(tài)方程和理想氣體狀態(tài)方程描述了水和空氣這2種介質(zhì)，通過處理歐拉方程組模擬了沖擊波載荷下空穴陣列在水中的變形過程，探討了水平空穴陣列的坍塌機(jī)制。Ozlem等[4]采用Godunov型計算方法，重點(diǎn)關(guān)注了沖擊波與空穴相互作用形成的高壓區(qū)域。Kapila等[]延續(xù)了Ozlem等[4]的工作，引入關(guān)于反應(yīng)率、未反應(yīng)炸藥和爆轟產(chǎn)物的平流方程，模擬了嵌入橢圓形空穴的HMX型固體炸藥在沖擊波加載下的爆轟過程。Kapahi等研究了強(qiáng)瞬態(tài)載荷下微米尺度多孔高能材料的響應(yīng)過程。Michael等通過整合多組分方程以及歐拉方程，采用MUSCL-Hancock方法討論了沖擊波誘發(fā)空穴的坍塌過程。Xiang等8采用剛性狀態(tài)方程描述雙流體，通過計算Godunov型數(shù)值通量模擬了含空穴圓柱形水柱與平面沖擊波的相互作用過程。Betney等利用拉格朗日超曲面捕捉材料之間的界面，討論了沖擊波與鑄入空穴水凝膠的相互作用。Sun 等[10]利用陡度可調(diào)諧波（steepness-adjustable harmonic，SAH）技術(shù)進(jìn)行重構(gòu)，考察了空穴列對非均勻含能材料LX-17沖擊起爆過程的影響。上述研究有助于理解空穴的坍塌機(jī)制，然而對于空穴壓縮、塌陷以及消失后爆轟波傳播全過程的數(shù)值模擬依然面臨困難和挑戰(zhàn)。可壓縮流體多組分問題常使用2種計算方法，一種是多相流方法[11-12]，另一種是求解歐拉方程組的增廣形式[13-14]，歐拉方程組的增廣形式包含質(zhì)量守恒、動量守恒、能量守恒方程以及描述混合物某種性質(zhì)的方程，比如化學(xué)組分的質(zhì)量分?jǐn)?shù)演化方程。Shyue[15-17]采用歐拉方程組作為控制方程，討論了針對剛性狀態(tài)方程、Van derWaals狀態(tài)方程和Mie-Gruneisen狀態(tài)方程的多組分模型，這是非常經(jīng)典的研究成果，并沒有考慮化學(xué)反應(yīng)。而本文中采用反應(yīng)歐拉方程組作為控制方程，計算區(qū)域被看作由兩部分組成，一部分是炸藥與爆轟產(chǎn)物的混合物，另一部分是不與爆炸混合物融合的惰性成分空氣。

對比多相流方法的混合界面模型，多介質(zhì)界面追蹤方法將界面的厚度看成零，能夠更精準(zhǔn)地捕捉到界面的變形過程，但是需要面臨不同介質(zhì)間的界面難題。常見的處理方法有任意拉格朗日歐拉（arbitraryLagrange Euler，ALE）方法[18]、流體體積（volume of fluid，VOF）方法[19-20]、界面追蹤（front tracking）方法[21]以及水平集（level set，LS）方法[2-23]。由于水平集方法能夠自然地描述界面的拓?fù)渥兓⑶夷芡ㄟ^高階方法處理水平集函數(shù)的控制方程，因此采用水平集/虛擬流體方法追蹤化學(xué)反應(yīng)混合物與空穴間的界面。傳統(tǒng)虛擬流體方法應(yīng)用起來很簡便，并且容易推廣至高維空間，但是不適用于求解強(qiáng)沖擊波或者高速射流沖擊等復(fù)雜問題。劉鐵鋼等[24]率先引人求解黎曼問題的思想，陸續(xù)發(fā)展了各種虛擬流體的變體方法，比如修正的虛擬流體方法[25]、真實(shí)的虛擬流體方法[2]以及實(shí)用的虛擬流體方法[2]等，更適用于處理含強(qiáng)激波或強(qiáng)間斷的計算難題。在修正的虛擬流體方法中，需要在界面附近的計算單元內(nèi)構(gòu)造并求解多介質(zhì)黎曼問題，給出接觸間斷的近似狀態(tài)，再根據(jù)預(yù)測值定義虛擬流體狀態(tài)。該方法已經(jīng)成功地應(yīng)用于可壓多介質(zhì)流、氣體化學(xué)反應(yīng)流、水下爆炸[28-30]等多種問題的數(shù)值模擬。然而，當(dāng)化學(xué)反應(yīng)涉及復(fù)雜形式的狀態(tài)方程時，只能運(yùn)用數(shù)值迭代方法實(shí)現(xiàn)守恒變量與原始變量的相互轉(zhuǎn)化，通過求解多介質(zhì)黎曼問題定義虛擬流體的狀態(tài)成為一種難題。本文中，運(yùn)用一種能夠處理復(fù)雜多介質(zhì)界面問題的歐拉計算方法，實(shí)現(xiàn)對炸藥內(nèi)空穴壓縮、塌陷以及消失后爆轟波傳播全過程的數(shù)值模擬。

1控制方程組

1.1 反應(yīng)歐拉方程組

考察如下形式的可壓縮反應(yīng)歐拉方程組：

其中：

式中： ρ，u，ν，p，E 和 λ 分別為密度、 x 方向的運(yùn)動速度、 y 方向的運(yùn)動速度、壓力、比總能以及爆轟產(chǎn)物的質(zhì)量分?jǐn)?shù)， R 為化學(xué)反應(yīng)率函數(shù)。比總能 E 滿足：

式中： e 為比內(nèi)能， 為熱量參數(shù)。

空氣與炸藥和爆轟產(chǎn)物的混合物被認(rèn)為是不會融合的，并且空氣不參與化學(xué)反應(yīng)。通過水平集方法追蹤化學(xué)反應(yīng)混合物與空氣的界面時，計算區(qū)域被分割為2個子區(qū)域。其中，將水平集函數(shù)大于零的區(qū)域定義為化學(xué)反應(yīng)混合物所在區(qū)域;將水平集函數(shù)小于零的區(qū)域定義為空氣所在區(qū)域。這2種流體的質(zhì)量、動量、總能量和爆轟產(chǎn)物的總質(zhì)量分?jǐn)?shù)關(guān)于反應(yīng)歐拉方程組（1）分別進(jìn)行更新，其中對于水平集函數(shù)小于零的介質(zhì)，即空氣，方程組（1）中的化學(xué)反應(yīng)率函數(shù) R 取為0，控制方程退化為無反應(yīng)歐拉方程組。

對于水平集函數(shù)大于零的介質(zhì)，即炸藥與爆轟產(chǎn)物的混合物，化學(xué)反應(yīng)率函數(shù)取為簡化的點(diǎn)火增長模型，它是依賴壓力項(xiàng)的單步反應(yīng)模型[：

R=G₁（1-λ）^c（λ+λ₀）^dp^y

式中： G₁，c，d，y 和 λ₀ 均為常數(shù)。

1.2 狀態(tài)方程

為使方程組（1）封閉，需要給出狀態(tài)方程的表達(dá)形式。硝基甲烷是一種液體炸藥，其內(nèi)部比較均勻，微介觀尺度下不存在裂紋、晶錯等結(jié)構(gòu)缺陷。因此，未反應(yīng)炸藥與爆轟產(chǎn)物可以用 Jones-Wilkins-Lee（JWL）狀態(tài)方程描述[7]：

式中： ω，ρ₀，A，B，R₁ 和 R₂ 為常數(shù)。用式（4）分別描述未反應(yīng)炸藥和爆轟產(chǎn)物時需要不同的狀態(tài)方程參數(shù)。

空氣滿足理想氣體狀態(tài)方程：

p=（γ-1）ρe

式中： γ 為比熱比。

2水平集方法和虛擬流體方法

2.1 水平集方法

通過水平集方法確定化學(xué)反應(yīng)混合物與空氣的界面，水平集函數(shù)的控制方程為對流方程[3]：

式中： ? 為水平集函數(shù)。零水平集標(biāo)志著界面所在處。

使用式（6更新下一時刻的水平集函數(shù)后，還需要根據(jù)下式：

對水平集函數(shù)進(jìn)行重新初始化，維持距離函數(shù)的性質(zhì)，式中 S 為符號函數(shù)。

2.2 虛擬流體的構(gòu)造方法

1999年，虛擬流體方法由Fedkiw等[32]首次提出，通過熵外推技術(shù)處理虛擬流體的密度。結(jié)合水平集和虛擬流體方法將多介質(zhì)問題轉(zhuǎn)化為單介質(zhì)問題，有助于提高算法的魯棒性，從而能夠運(yùn)用高精度數(shù)值方法分別處理各流體。沖擊波加載下的空穴塌陷以及爆轟波傳播過程是含強(qiáng)激波和強(qiáng)間斷的計算難題，修正的虛擬流體方法考慮了界面附近各介質(zhì)的物質(zhì)屬性，因此采用修正的虛擬流體方法處理多介質(zhì)界面。在界面附近的計算單元內(nèi)引入黎曼問題，通過求解接觸間斷兩側(cè)的物理量預(yù)測虛擬流體的狀態(tài)。

首先通過延拓方程，將當(dāng)?shù)亓黧w的變量外推到虛擬流體所在計算單元：

式中： V 為需要構(gòu)造的虛擬流體變量； τ 為人工時間步長； N 為單位法向向量，法向定義為由負(fù)水平集指向正水平集。

在二維空間，需要基于真實(shí)流體和虛擬流體的狀態(tài)沿著界面法方向構(gòu)造黎曼問題。通常，各介質(zhì)的狀態(tài)方程是確定的，求解黎曼問題并沒有特別困難。然而，由于本文中需要處理的未反應(yīng)炸藥和爆轟產(chǎn)物具有不同的狀態(tài)方程參數(shù)，不得不通過一些數(shù)值迭代方法實(shí)現(xiàn)守恒變量和原始變量的相互轉(zhuǎn)化，這使得構(gòu)造黎曼解成為一種難題。JWL狀態(tài)方程與理想氣體狀態(tài)方程形式差別很大，特別是JWL狀態(tài)方程形式較復(fù)雜，HLLC（Harten-Lax-vanLer contact）近似黎曼解法器僅通過Rankine-Hugoniot關(guān)系式就可以計算出接觸間斷區(qū)域的物理變量，因此基于HLLC黎曼解法器發(fā)展了一種涉及化學(xué)反應(yīng)率的虛擬流體變量求解方法。

接下來將介紹HLLC黎曼解法器，假設(shè)黎曼解結(jié)構(gòu)中由密度、速度、壓力組成的原始變量從左至右依次用 U_L， U_L^* 、 U_R^* 和 U_R 表示，那么如圖1所示，黎曼解滿足[33]：

黎曼解中3個波兩側(cè)的變量狀態(tài)均滿足Rankine-Hugoniot關(guān)系式：

圖1近似黎曼解結(jié)構(gòu)示意圖

Fig.1Schematic diagram of approximate Riemann solution structure

F_L^*=F_L+S_L（U_L^*-U_L）

中間的波為接觸間斷，還滿足以下關(guān)系式：

綜合式 （10）～（13） ，可知接觸間斷兩側(cè)變量為：

其中最慢波速和最快波速能夠通過以下關(guān)系式近似地求得：

S_L=min{u_L-c_L，u_R-c_R}

S_R=max{u_L+c_L，u_R+c_R}

在本文中，虛擬流體變量狀態(tài)的求解過程可以分為以下4個步驟。

（1）計算接觸間斷區(qū)域的變量：

（2）令 λ_K^*=λ_K （K=L，R） ，根據(jù)守恒變量 ρ_K^*?ρ_K^*u_K^* 和 ρ_K^*E_K^* 計算原始變量 p_K^* 。對于爆炸混合物，這個過程需要聯(lián)立以下6個方程，即炸藥與爆轟產(chǎn)物各自的狀態(tài)方程、溫壓平衡條件、混合內(nèi)能和混合比容關(guān)系式，得到未反應(yīng)炸藥比容與爆炸產(chǎn)物比容滿足的非線性方程，再通過一些求解非線性方程的數(shù)值方法如牛頓迭代法處理該方程，計算出混合壓力。

（3）重置黎曼問題初始條件：

ρ_K=ρ_K^*，u_K=u^*，p_K=p_K^*K=L，R

（4）如果 |p_L-p_R| 大于可允許誤差，重新進(jìn)入步驟（1）;否則跳出循環(huán)， ρ_K^*?p^*?u^* 、 λ_K^*（K=L，R） 即為預(yù)測的虛擬流體狀態(tài)。

3離散方法

方程組（1）的半離散格式為：

式中： F_i-1/2，j 和 F_i+1/2，j 為沿著 x 方向的半節(jié)點(diǎn)處的數(shù)值通量， G_i，j-1/2 和 G_i，j+1/2 為沿著 y 方向的半節(jié)點(diǎn)處的數(shù)值通量， S_i，j 為節(jié)點(diǎn)處源項(xiàng)的數(shù)值近似。

采用三階TVDRunge-Kutta方法進(jìn)行時間離散，應(yīng)用五階WENO-LF有限差分方法[34]求解數(shù)值通量F_i±1/2，j 和 G_i，j±1/2 。以 x 方向?yàn)槔瑢ENO-LF有限差分方法進(jìn)行簡要描述。將節(jié)點(diǎn)上的通量 F_i，j 進(jìn)行Lax-Friedrichs分裂：

式中： α=max（|λ（U_i，j）|），λ（U_i，j） 為雅可比矩陣 U的特征值。假設(shè)矩陣L和R分別為矩陣A的左、右特征矩陣，令 W_i，j^±=L_AF_i，j^± ，先將數(shù)值通量投影到局部特征空間，對 W_i+1/2，j^± 分別進(jìn)行重構(gòu)，然后再將其變換回物理空間F+1/2，j=RAW+1/2，。

將反應(yīng)歐拉方程組進(jìn)行特征分解，計算獲得雅可比矩陣為：

由此可知，求解雅可比矩陣的左右特征向量以及特征值時，不可避免地需要計算出壓力關(guān)于守恒量的偏導(dǎo)數(shù)。在化學(xué)反應(yīng)區(qū)域，由于同時存在炸藥和爆轟產(chǎn)物2種物質(zhì)，無法直接給出統(tǒng)一形式的狀態(tài)方程。因此，設(shè)置溫壓平衡條件：

式中：下標(biāo)r、p分別表示炸藥與爆轟產(chǎn)物的對應(yīng)變量。混合內(nèi)能和混合比容關(guān)系式滿足：

基于方程組 （23）～（24） 以及炸藥和爆轟產(chǎn)物各自的狀態(tài)方程，能夠通過求解線性方程組計算出2種比容 τ_r，τ_p 關(guān)于各守恒量的偏導(dǎo)數(shù)，再獲得壓力關(guān)于各守恒量的偏導(dǎo)數(shù)[35]。

4沖擊波加載下炸藥內(nèi)空穴塌陷過程的數(shù)值模擬

硝基甲烷及其爆轟產(chǎn)物的狀態(tài)方程和反應(yīng)率函數(shù)的相關(guān)參數(shù)參考文獻(xiàn)[7]。空穴直徑為 1mm ，比熱比為1.4，初始密度、 x 方向與 y 方向的運(yùn)動速度和壓力分別為 1.205kg/m³ 、 0m/s 、 0m/s 和 101.3kPa 。

在空穴左側(cè)設(shè)置入射沖擊波，考察了4、6和 8GPa 這3種壓力加載下的空穴塌陷過程。計算區(qū)域左邊界為入流邊界條件，其余邊界均為出流邊界條件。計算網(wǎng)格尺寸為 Δx=Δy=1mm/90 。

圖 2～10 展示了 4GPa 沖擊壓力下的計算結(jié)果，其中圖 2～3 為空穴壓縮以及塌陷過程中不同時刻的密度和壓力分布。由于這段時間內(nèi)爆轟產(chǎn)物的質(zhì)量分?jǐn)?shù)一直很小，因此不作展示。如圖2所示，沖擊波與空穴開始相互作用后，空穴的前表面受壓縮后向右移動發(fā)生形變，空穴對沖擊波產(chǎn)生稀疏作用，使得左側(cè)的壓力和密度下降。當(dāng)前表面與后表面接近時，空穴周圍的密度和壓力開始上升。隨著空穴被進(jìn)一步壓縮發(fā)生塌陷，會出現(xiàn)局部高溫、高壓、高密度區(qū)域，如圖3所示，熱點(diǎn)附近區(qū)域的壓力達(dá)到6GPa以上。

圖24GPa沖擊壓力下不同時刻的密度分布

Fig.2Density distributions under the impact pressure of4GPaat different times

圖34GPa沖擊壓力下不同時刻的壓力分布

Fig.3Pressure distributions under the impact pressure of 4GPaat different times

圖 4～6 展示的是空穴塌縮后3個時刻的密度、壓力和爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)分布。熱點(diǎn)附近區(qū)域的壓力更高，因此化學(xué)反應(yīng)速率更大，在熱點(diǎn)附近形成新的沖擊波后，向四周進(jìn)行擴(kuò)散。但是，新產(chǎn)生的沖擊波強(qiáng)度并不足以引發(fā)爆轟。 2.20μs 時，沖擊波波陣面處壓力峰值約為5GPa，遠(yuǎn)低于CJ爆轟壓力，最高的爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)也低于0.03，可以認(rèn)為幾乎沒有發(fā)生化學(xué)反應(yīng)。

圖44GPa沖擊壓力下不同時刻的密度分布

Fig.4Density distributions under the impact pressure of 4 GPa at different times

圖54GPa沖擊壓力下不同時刻的壓力分布

Fig.5Pressure distributionsunder the impact pressureof4GPaat different times

圖64GPa沖擊壓力下不同時刻的爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)分布

Fig.6Detonation product massfraction distributions under the impact pressure of 4 GPaat different times

圖 7～8 展示的是不同時刻中軸線處的密度、壓力和 x 方向運(yùn)動速度分布。從計算結(jié)果可以看出：起初受到稀疏波的影響， 0.17μs 時的密度和壓力從左至右下降。隨后空穴進(jìn)一步受到壓縮發(fā)生塌縮，0.37μs 時在 7.5～7.6mm 處形成能量匯聚的高密度、高壓力區(qū)域，壓力峰值約為 7.5GPa ，速度峰值約為3.6km/s. ，此時的壓力峰值明顯低于CJ爆轟壓力，因此化學(xué)反應(yīng)強(qiáng)度較低。空穴塌陷引起的沖擊波往四周擴(kuò)散，沖擊波波陣面處的壓力峰值逐漸下降， 2.20μs 時壓力峰值約為 4.8GPa 。

圖74GPa沖擊壓力下0、0.17和 0.37μs 時的中軸線密度、壓力和速度分布 Fig.7Density，pressure and velocitydistributions on the centralaxis under the impact pressure of4 GPa at 0，0.17，and 0.37μs

圖84GPa沖擊壓力下0.74、1.35和 2.20μs 時的中軸線密度、壓力和速度分布 Fig.8Density，pressure and velocity distributions on the centralaxis under the impact pressure of 4 GPa at 0.74，1.35，and 2.20μs

圖 9～10 展示了不同網(wǎng)格尺寸下 0.37μs 時的計算結(jié)果。顯然，網(wǎng)格數(shù)量越大，計算區(qū)域內(nèi)流場分布以及氣泡閉合與熱點(diǎn)產(chǎn)生過程刻畫得越清晰。如圖10所示，將網(wǎng)格尺寸為 Δx=Δy=1mm/180 的計算結(jié)果當(dāng)作參考解，當(dāng)網(wǎng)格尺寸為 Δx=Δy=1mm/22.5 時，捕捉不到準(zhǔn)確的空穴塌陷過程。隨著計算網(wǎng)格的加密，中軸線處的密度和壓力分布曲線與參考解越吻合。特別是流場分布較平緩的區(qū)域，網(wǎng)格尺寸為 Δx=Δy=1mm/90 與 Δx=Δy=1mm/180 計算出的密度和壓力分布曲線幾乎重合。數(shù)值結(jié)果顯示：Δx=Δy=1mm/22.5 ， Δx=Δy=1mm/45 Δx=Δy=1mm/90 與 Δx=Δy=1mm/180 計算得到的密度峰值分別相差0.377、0.092和 0.0285g/cm³ ，壓力峰值分別相差4.856、1.781和 0.85GPa， ，能夠說明本文算法的有效性。

圖9不同網(wǎng)格尺寸下 0.37μs 時的密度分布

Fig.9Densitydistributionswith different mesh sizesat 0.37μs

圖10 0.37μs 時的中軸線密度和壓力分布

Fig.10 Density and pressure distributions on the central axis at 0.37μs

圖 11～17 展示了6GPa沖擊壓力下7個時刻的計算結(jié)果，其中圖 11～12 為空穴壓縮以及塌陷過程中不同時刻的密度和壓力分布。空穴被壓縮發(fā)生塌陷，會出現(xiàn)高溫、高壓、高密度的熱點(diǎn)區(qū)域，局部區(qū)域壓力達(dá)到 9GPa 以上。

圖116GPa沖擊壓力下不同時刻的密度分布

Fig.11Density distributionsunder the impact pressure of6GPaat different times

圖126GPa沖擊壓力下不同時刻的壓力分布

Fig.12 Pressure distributionsunder the impact pressure of6 GPaat different times

圖136GPa沖擊壓力下不同時刻的密度分布

Fig.13Density distributions under the impact pressure of 6 GPa at different times

圖146GPa壓力下不同時刻的壓力分布

Fig.14Pressure distributions under the impact pressure of6GPa at different times

圖156GPa沖擊壓力下不同時刻的爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)分布

Fig.15Detonation product mass fraction distributions under the impact pressure of6GPaat different times

圖166GPa沖擊壓力下0、0.15和 0.30μs 時的中軸線密度、壓力和速度分布

圖176GPa沖擊壓力下0.70、1.11、1.94和 2.49μs 時的中軸線密度、壓力和速度分布Fig.17Density，pressureand velocity distributions onthe central axis under the impact pressre of 6 GPaat 0.70，1.11，1.94，and 2.49μs （20號

圖 13～15 展示的是空穴塌縮后4個時刻的密度、壓力和爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)分布。能夠看出：在熱點(diǎn)附近化學(xué)反應(yīng)速率顯著提高，形成新的沖擊波向四周傳播。在沖擊波作用下，高壓區(qū)域內(nèi)的炸藥發(fā)生化學(xué)反應(yīng)，快速生成新的爆轟產(chǎn)物。如圖 14～15 所示，隨著時間演化，沖擊波波陣面處壓力上升， 1.94μs 時，壓力峰值為 10GPa 2.49μs 時，壓力峰值增長到遠(yuǎn)高于CJ爆轟壓力的 22GPa 。高壓區(qū)域爆轟產(chǎn)物的質(zhì)量分?jǐn)?shù)持續(xù)上升， 2.49μs 時， 2.7～16.4mm 區(qū)域內(nèi)的炸藥幾乎全部轉(zhuǎn)化為爆轟產(chǎn)物。

圖 16～17 展示的是不同時刻中軸線密度、壓力和 x 方向運(yùn)動速度的分布。起初隨著空穴被壓縮，炸藥受到稀疏作用，密度以及壓力從左至右下降。空穴進(jìn)一步受到壓縮發(fā)生塌陷， 0.30μs 時，在 10.5～ 10.6mm 處形成能量匯聚的高密度、高壓力區(qū)域，壓力峰值約為 10GPa ，速度峰值約為 4.4km/s ，熱點(diǎn)附近區(qū)域內(nèi)的化學(xué)反應(yīng)速率大幅提升，產(chǎn)生新的沖擊波。沖擊波壓力峰值起先略微下降，然后隨著化學(xué)反應(yīng)的進(jìn)行不斷提高。 1.94μs 時，新產(chǎn)生的沖擊波傳播到 5mm 處，該處壓力值為 10.5GPa，2.49μs 時，傳播到 3mm 處，該處壓力值為 21.7GPa 。

圖 18～24 展示了8GPa沖擊壓力下7個時刻的計算結(jié)果，其中圖 18～19 為空穴壓縮以及塌陷過程中不同時刻的密度和壓力分布。沖擊波與空穴相互作用后，空穴的前表面受壓縮后發(fā)生形變。隨著空穴被進(jìn)一步壓縮，熱點(diǎn)區(qū)域壓力達(dá)到 11GPa 以上。

圖 20～22 展示的是空穴塌縮后4個時刻的密度、壓力和爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)分布。空穴塌縮后，熱點(diǎn)附近的壓力和溫度升高，化學(xué)反應(yīng)速率迅速升高。隨后，由熱點(diǎn)引發(fā)新的沖擊波向四周傳播，沖擊波波后高壓區(qū)的炸藥發(fā)生劇烈反應(yīng)，如圖 21～22 所示， 0.78μs 時壓力峰值高達(dá) 25GPa ，部分區(qū)域內(nèi)的炸藥已完全反應(yīng)，轉(zhuǎn)化為爆轟產(chǎn)物。 0.90μs 時壓力峰值升至 27GPa ，在加載沖擊波以及新產(chǎn)生沖擊波的作用下，高壓區(qū)的剩余炸藥繼續(xù)劇烈反應(yīng)，轉(zhuǎn)化為爆轟產(chǎn)物。 1.50μs 時 14mm 左側(cè)區(qū)域的炸藥幾乎全部發(fā)生反應(yīng)。

圖188GPa沖擊壓力下不同時刻的密度分布

Fig.18 Density distributions under the impact pressure of 8GPa at different times

圖198GPa沖擊壓力下不同時刻的壓力分布

圖208GPa沖擊壓力下不同時刻的密度分布

Fig.20 Density distributions under the impact pressure of8GPaat different times

圖228GPa沖擊壓力下不同時刻的爆轟產(chǎn)物質(zhì)量分?jǐn)?shù)分布

Fig.22Detonation product mass fraction distributions under the impact pressureof8GPaat diferent times

圖238GPa沖擊壓力下 0，0.14 和 0.27μs 時的中軸線密度、壓力和速度分布

圖248GPa沖擊壓力下0.58、0.78、0.90和 1.50μs 時的中軸線密度、壓力和速度分布 Fig.24Density，pressure and velocity distributions on the central axis under the impact pressure of 8GPa at 0.58，0.78，0.90，and 1.50μs

圖 23～24 展示的是不同時刻的中軸線密度、壓力和 x 方向運(yùn)動速度分布。 0.27μs 時 7.50～7.67mm 處為能量匯聚的高密度、高壓力區(qū)，速度峰值約為 4.9km/s 。沖擊波波陣面往四周傳播， 0.78μs 時在5.50mm 處沖擊壓力已經(jīng)達(dá)到 25GPa， ，高于約 14GPa 的CJ爆轟壓力，隨后壓力峰值一直高于CJ爆轟壓力。

通過數(shù)值模擬3種不同強(qiáng)度沖擊波加載下的液體炸藥硝基甲烷與氣泡相互作用過程，可以得到結(jié)論：在熱點(diǎn)即空穴塌縮處能夠?qū)崿F(xiàn)化學(xué)反應(yīng)的加速，并且化學(xué)反應(yīng)的劇烈程度隨著沖擊波壓力的增加而增大。 4GPa 沖擊壓力加載下，熱點(diǎn)產(chǎn)生后爆轟產(chǎn)物的質(zhì)量分?jǐn)?shù)增長非常緩慢。而 6GPa 沖擊壓力加載下，炸藥內(nèi)空穴塌縮形成熱點(diǎn)，能夠?qū)崿F(xiàn)爆轟，引起化學(xué)反應(yīng)區(qū)域的擴(kuò)散。這與Turley等的實(shí)驗(yàn)結(jié)果[36]一致，即沖擊波強(qiáng)度不夠高時，化學(xué)反應(yīng)區(qū)域會停止增長，而當(dāng)沖擊壓力超過5GPa時，在空穴塌縮處發(fā)生爆轟，爆轟產(chǎn)物區(qū)域會逐漸向四周擴(kuò)大。

5結(jié)論

提出了一種處理炸藥與爆轟產(chǎn)物的混合物以及空穴間界面問題的虛擬流體方法，并且采用基于水平集技術(shù)的多介質(zhì)界面追蹤方法對炸藥內(nèi)空穴壓縮、塌陷等介觀尺度下的動力學(xué)過程開展了數(shù)值模擬。模擬結(jié)果表明：在沖擊波載荷下，含空穴液體炸藥硝基甲烷比純凈炸藥的化學(xué)反應(yīng)速率顯著提高。在4GPa的沖擊壓力加載下，由于空穴塌陷形成新沖擊波的強(qiáng)度不足以引發(fā)爆轟，生成的爆轟產(chǎn)物非常少。當(dāng)沖擊壓力升至6GPa時，空穴塌縮、閉合處炸藥的化學(xué)反應(yīng)速率大幅提高，成功實(shí)現(xiàn)爆轟，反應(yīng)區(qū)內(nèi)不斷產(chǎn)生爆轟產(chǎn)物。含空穴硝基甲烷在8GPa的沖擊壓力加載下能更快實(shí)現(xiàn)爆轟。

參考文獻(xiàn)：

[1] 彭亞晶，葉玉清.含能材料起爆過程“熱點(diǎn)”理論研究進(jìn)展[J].化學(xué)通報，2015，78（8）：693-701.DOI：10.14159/j.cnki. 0441-3776.2015.08.004. PENG Y J， YEYQ.Research progress of ‘hot-spot’ theory in energetic materials initiation[J].Chemistry，2015，78（8）： 693-701. DO1： 10.14159/j.cnki.0441-3776.2015.08.004.

[2] MADER CL.Numerical modeling of detonations [M].Berkeley，USA： University ofCalifornia Press，1979.

[3] LAUERE，HUXY，HCKELS，etal.Numericalinvestigationofcollapsingcavityarrays[J].hsicsofFids，1224（5）： 052104. DOI： 10.1063/1.4719142.

[4] OZLEM M， SCHWENDEMAN D W， KAPILA A K， et al. A numerical study of shock-induced cavity collpse [J]. Shock Waves，2012，22（2）： 89-117.DOI：10.1007/s00193-011-0352-9.

[5] KAPILA A K，SCHWENDEMANDW，GAMBINOJR，et al.A numerical studyof thedynamics of detonation initiated by cavity collapse [J]. Shock Waves，2015，25（6）： 545-572.DOI： 10.1007/s00193-015-0597-9.

[6] KAPAHI A，UDAYKUMAR H S.Dynamicsof void collapse inshockedenergetic materials：physics ofvoid-void interactions[J]. Shock Waves，2013，23（6）： 537-558. DOI： 10.1007/s00193-013-0439-6.

[7] MICHAEL L，NIKIFORAKIS N.A hybrid formulation for the numerical simulation of condensed phase explosives [J]. Journal of Computational Physics，2016，316： 193-217. DO1： 10.1016/j.jcp.2016.04.017.

[8] XIANG GM，WANGB.Numerical studyofa planar shock interacting withacylindrical watercolumn embedded withan air cavity[J]. Journal ofFluidMechanics，2017，825：825-852.DOI： 10.1017/jfm.2017.403.

[9] BETNEY MR，ANDERSON PA，DOYLE H， etal. Numerical and experimental studyof shock-driven cavitycollapse [C]/ 9th Intermational Symposium on Cavitation （CAV2015）.Lausanne： IOP Publishing， 2015： 01201. DOI： 10.1088/1742-6596/ 656/1/012011.

[10]SUNJ，YANGPF，MENGB Q，etal.Effectsofmicrometer-scalecavitiesontheshock-to-detonationtransitionia heterogeneous LX-17 energetic material[J]. Physics ofFluids，2023，35（12）： 126117.DOI： 10.1063/5.0174851.

[11]于明.固體炸藥爆轟與惰性介質(zhì)相互作用的一種擴(kuò)散界面模型[J].爆炸與沖擊，2020，40（10）：104202.DOI：10.11883/ bzycj-2019-0435. YUM.Animproved difuse interfacemodel for the numericalsimulation of interaction between solid explosive detonation and inert media[J]. Explosion and Shock Waves，2020， 40（10）： 104202. DOI： 10.11883/bzycj-2019-0435.

[12]SCHMIDMAYER K， BRYNGELSON S H， COLONIUS T. An assessment of multcomponent flow models and interface capturingschemes forsphericalbubbledynamics[J].JournalofComputational Physics，2020，402：109080.DOI：10.016/J. JCP.2019.109080.

[13]XIAO M，NI G X， WANG C，etal. Front capturing by level set method for thereactive Eulerequations[J]. Inteational Journal for Numerical Methods in Fluids， 2021， 93（8）： 2723-2743. DO1： 10.1002/fld.4995.

[14]KLOPSCHR，GARANN，BACHE，etal.Parametric influenceonrotatingdetonationcombustion：insights fromfstactive Euler simulations [J].AIAA Journal， 2024， 62（1）： 127-139. DOI： 10.2514/1.J063193.

[15]SHYUE K M.An efficient shock-capturing algorithm forcompressble multicomponent problems[J]. Journalof Computational Physics， 1998， 142（1）： 208-242. DOI： 10.1006/jcph.1998.5930.

[16]SHYUE KM.Afluid-mixture typealgorithmforcompresiblemulticomponentfow withvan der Waals equationof state[J]. Journal of Computational Physics，1999，156（1）： 43-88. DOI： 10.1006/jcph.1999.6349.

[1/」SHYUE K M. A Iua-mixture type aIgoritnm Ior compressipie muucomponent IIow win Me-Gruneisen equauon oI state[J]. Journal of Computational Physics， 2001，171（2）： 678-707. DOI： 10.1006/jcph.2001.6801.

[18]ZHANGF，CHENGJ.Abound-preservingand positivity-preserving high-order arbitraryLagrangian-Eulerian discontiuous Galerkin method for compressble multi-medium flows[J].SIAM Journal on Scientific Computing， 2024，46（3）： B254B279. DOI： 10.1137/23M1588810.

[19]HIRTCW，NICHOLS BD.Volumeofuid（VOF）method forthe dynamicsoffreeboundaries[J]. JournalofComputational Physics，1981，39（1）： 201-225. DO1： 10.1016/0021-991（81）90145-5.

[20]MULBAH C，KANG C，MAO N， et al. A reviewof VOF methods for simulating bubble dynamics[J]. Progressin Nuclear Energy，2022，154： 104478. DO1： 10.1016/j.pnucene.2022.104478.

[21]TRYGGVASON G，BUNNERB，ESMAEELI A，etal.A front-tracking method forthecomputations of multiphase flow[J]. Journal of Computational Physics， 2001，169（2）： 708-759. DOI： 10.1006/jcph.2001.6726.

[22]GIBOUF，F(xiàn)EDKIWR，OSHER S.Areviewof level-set methods andsomerecent applications [J]. Jourmal ofComputational Physics，2018，353： 82-109. DO1： 10.1016/j jcp.2017.10.006.

[23]姚成寶，王宏亮，浦錫鋒，等.空中強(qiáng)爆炸沖擊波地面反射規(guī)律數(shù)值模擬研究[J].爆炸與沖擊，2019，39（11）：112201.DOI： 10.11883/bzycj-2018-0287. YAO CB， WANG HL，PUXF，et al. Numerical simulation ofintense blast wave reflected onrigid ground[J].Explosionand Shock Waves， 2019，39（11）： 112201. DOI： 10.11883/bzycj-2018-0287.

[24]劉鐵鋼，許亮.模擬多介質(zhì)界面問題的虛擬流體方法綜述[J].氣體物理，2019，4（2）：1-16.DOI：10.19527/j.cnki.2096- 1642.0746. LIUT G， XUL.Areview of ghost fluid methods for multi-medium interface simulation [J]. Physics of Gases，2019，4（2）： 1-16. DOI： 10.19527/j.cnki.2096-1642.0746.

[25]LIUTG， KHOO B C，YEOK S. Ghost fluid method for strong shock impacting on material interface [J]. Journal of Computati0nal Physics，2003，190（2）： 651-681. DOI： 10.1016/S0021-9991（03）00301-2.

[26] WANG C W，LIUT G， KHOO B C.A real ghost fluid method for the simulation of multimedium compressble flow [J]. SIAM Journal on Scientific Computing， 2006， 28（1）： 278-302. DOI： 10.1137/030601363.

[27]XU L，F(xiàn)ENG CL，LIUTG.Practical techniques in ghost fluid method forcompressble multi-medium flows[J]. Communications in Computational Physics， 2016，20（3）： 619-659. DOI： 10.4208/cicp.190315.290316a.

[28]HUO ZX，JIAZP.AGRP-based tangential effects preserving，highresolution and eficient ghostfluid method forthe simulationof two-dimensional multi-medium compressible flows [J]. Computers and Fluids，2024，276：106261.DOI： 10.1016/j.compfluid.2024.106261.

[29]XUL，YANG WB，LIUTG.Aninterface treatment for two-material multi-species flows ivolving thermall perfect gases with chemical reactions [J]. Jourmal of Computational Physics，2022，448： 110707.DOI： 10.1016/jjcp.2021.110707.

[30]ZHAO Z T，RONGJL， ZHANG S X.Anumerical studyofunderwater explosions basedon the ghostfluid method[J].Ocean Engineering，2022， 247： 109796. DO1： 10.1016/j.0ceaneng.2021.109796.

[31]OSHER S，SETHIANJA.Frontspropagating withcurvature-dependent speed：algorithms basedon Hamilton-Jacobi formulations [J]. Jourmal ofComputational Physics，1988，79（1）： 12-49.DOI： 10.1016/0021-9991（88）90002-2.

[32]FEDKIWR P，ASLAMT，MERRIMANB，etal.Anon-oscilatory Eulerianaproach to interfaces in multimaterialflows （the ghost fluid method）[J]. Joumal ofComputational Physics，1999，152（2）： 457-492. DO1： 10.1006/jcph.1999.6236.

[33]TOROEF.Riemann solvers and numerical methods forfluid dynamics：a practical introduction[M].3rded.Heidelberg： Springer Scienceamp; Business Media，2009.

[34]JIANGGS，SHU CW.Eficient implementation ofweighted ENO schemes[J]. Joural of Computational Physics，996， 126（1）： 202-228. DO1： 10.1006/jcph.1996.0130.

[35]WANGC，LIUXQ.Highresolutionnumericalsimulationofdetonation difractionofcondensed explosives[J].Inteational Journal of Computational Methods，2015，12（2）： 1550005. DOI： 10.1142/S021987621550005X.

[36]TURLEYWD，LALONEBM，MANCEJG，etal.Experimentalobservationsofshock-wave-induced bubblecollapsand hot-spot formation innitromethane liqud explosive[J].Journal of Appied Physics，2021，129（14）：145102.DOI：10.1063/ 5 0039414

（責(zé)任編輯 張凌云）