含空穴炸藥硝基甲烷沖擊轉爆轟過程的數值模擬

中圖分類號：0381 國標學科代碼：13035 文獻標志碼：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

炸藥在澆鑄、壓裝過程中很難避免在內部形成空穴。根據熱點理論，非均質炸藥在沖擊波作用下，空穴受壓塌陷產生局部高溫高壓區域，當熱點區域足夠大或者高溫持續時間足夠長時，會引起炸藥局部發生點火，實現爆轟，點火時間會明顯短于點燃純凈炸藥所需的時間。因此，考察炸藥對于空穴的敏感性具有重要的研究意義和應用價值。

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

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

1控制方程組

1.1 反應歐拉方程組

考察如下形式的可壓縮反應歐拉方程組：

其中：

式中： ρ，u，ν，p，E 和 λ 分別為密度、 x 方向的運動速度、 y 方向的運動速度、壓力、比總能以及爆轟產物的質量分數， R 為化學反應率函數。比總能 E 滿足：

式中： e 為比內能， 為熱量參數。

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

對于水平集函數大于零的介質，即炸藥與爆轟產物的混合物，化學反應率函數取為簡化的點火增長模型，它是依賴壓力項的單步反應模型[：

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

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

1.2 狀態方程

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

式中： ω，ρ₀，A，B，R₁ 和 R₂ 為常數。用式（4）分別描述未反應炸藥和爆轟產物時需要不同的狀態方程參數。

空氣滿足理想氣體狀態方程：

p=（γ-1）ρe

式中： γ 為比熱比。

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

2.1 水平集方法

通過水平集方法確定化學反應混合物與空氣的界面，水平集函數的控制方程為對流方程[3]：

式中： ? 為水平集函數。零水平集標志著界面所在處。

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

對水平集函數進行重新初始化，維持距離函數的性質，式中 S 為符號函數。

2.2 虛擬流體的構造方法

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

首先通過延拓方程，將當地流體的變量外推到虛擬流體所在計算單元：

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

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

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

黎曼解中3個波兩側的變量狀態均滿足Rankine-Hugoniot關系式：

圖1近似黎曼解結構示意圖

Fig.1Schematic diagram of approximate Riemann solution structure

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

中間的波為接觸間斷，還滿足以下關系式：

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

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

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

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

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

（1）計算接觸間斷區域的變量：

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

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

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

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

3離散方法

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

式中： F_i-1/2，j 和 F_i+1/2，j 為沿著 x 方向的半節點處的數值通量， G_i，j-1/2 和 G_i，j+1/2 為沿著 y 方向的半節點處的數值通量， S_i，j 為節點處源項的數值近似。

采用三階TVDRunge-Kutta方法進行時間離散，應用五階WENO-LF有限差分方法[34]求解數值通量F_i±1/2，j 和 G_i，j±1/2 。以 x 方向為例，對WENO-LF有限差分方法進行簡要描述。將節點上的通量 F_i，j 進行Lax-Friedrichs分裂：

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

將反應歐拉方程組進行特征分解，計算獲得雅可比矩陣為：

由此可知，求解雅可比矩陣的左右特征向量以及特征值時，不可避免地需要計算出壓力關于守恒量的偏導數。在化學反應區域，由于同時存在炸藥和爆轟產物2種物質，無法直接給出統一形式的狀態方程。因此，設置溫壓平衡條件：

式中：下標r、p分別表示炸藥與爆轟產物的對應變量。混合內能和混合比容關系式滿足：

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

4沖擊波加載下炸藥內空穴塌陷過程的數值模擬

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

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

圖 2～10 展示了 4GPa 沖擊壓力下的計算結果，其中圖 2～3 為空穴壓縮以及塌陷過程中不同時刻的密度和壓力分布。由于這段時間內爆轟產物的質量分數一直很小，因此不作展示。如圖2所示，沖擊波與空穴開始相互作用后，空穴的前表面受壓縮后向右移動發生形變，空穴對沖擊波產生稀疏作用，使得左側的壓力和密度下降。當前表面與后表面接近時，空穴周圍的密度和壓力開始上升。隨著空穴被進一步壓縮發生塌陷，會出現局部高溫、高壓、高密度區域，如圖3所示，熱點附近區域的壓力達到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個時刻的密度、壓力和爆轟產物質量分數分布。熱點附近區域的壓力更高，因此化學反應速率更大，在熱點附近形成新的沖擊波后，向四周進行擴散。但是，新產生的沖擊波強度并不足以引發爆轟。 2.20μs 時，沖擊波波陣面處壓力峰值約為5GPa，遠低于CJ爆轟壓力，最高的爆轟產物質量分數也低于0.03，可以認為幾乎沒有發生化學反應。

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

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

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

Fig.5Pressure distributionsunder the impact pressureof4GPaat different times

圖64GPa沖擊壓力下不同時刻的爆轟產物質量分數分布

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

圖 7～8 展示的是不同時刻中軸線處的密度、壓力和 x 方向運動速度分布。從計算結果可以看出：起初受到稀疏波的影響， 0.17μs 時的密度和壓力從左至右下降。隨后空穴進一步受到壓縮發生塌縮，0.37μs 時在 7.5～7.6mm 處形成能量匯聚的高密度、高壓力區域，壓力峰值約為 7.5GPa ，速度峰值約為3.6km/s. ，此時的壓力峰值明顯低于CJ爆轟壓力，因此化學反應強度較低。空穴塌陷引起的沖擊波往四周擴散，沖擊波波陣面處的壓力峰值逐漸下降， 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 展示了不同網格尺寸下 0.37μs 時的計算結果。顯然，網格數量越大，計算區域內流場分布以及氣泡閉合與熱點產生過程刻畫得越清晰。如圖10所示，將網格尺寸為 Δx=Δy=1mm/180 的計算結果當作參考解，當網格尺寸為 Δx=Δy=1mm/22.5 時，捕捉不到準確的空穴塌陷過程。隨著計算網格的加密，中軸線處的密度和壓力分布曲線與參考解越吻合。特別是流場分布較平緩的區域，網格尺寸為 Δx=Δy=1mm/90 與 Δx=Δy=1mm/180 計算出的密度和壓力分布曲線幾乎重合。數值結果顯示：Δ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不同網格尺寸下 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個時刻的計算結果，其中圖 11～12 為空穴壓縮以及塌陷過程中不同時刻的密度和壓力分布。空穴被壓縮發生塌陷，會出現高溫、高壓、高密度的熱點區域，局部區域壓力達到 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沖擊壓力下不同時刻的爆轟產物質量分數分布

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個時刻的密度、壓力和爆轟產物質量分數分布。能夠看出：在熱點附近化學反應速率顯著提高，形成新的沖擊波向四周傳播。在沖擊波作用下，高壓區域內的炸藥發生化學反應，快速生成新的爆轟產物。如圖 14～15 所示，隨著時間演化，沖擊波波陣面處壓力上升， 1.94μs 時，壓力峰值為 10GPa 2.49μs 時，壓力峰值增長到遠高于CJ爆轟壓力的 22GPa 。高壓區域爆轟產物的質量分數持續上升， 2.49μs 時， 2.7～16.4mm 區域內的炸藥幾乎全部轉化為爆轟產物。

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

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

圖 20～22 展示的是空穴塌縮后4個時刻的密度、壓力和爆轟產物質量分數分布。空穴塌縮后，熱點附近的壓力和溫度升高，化學反應速率迅速升高。隨后，由熱點引發新的沖擊波向四周傳播，沖擊波波后高壓區的炸藥發生劇烈反應，如圖 21～22 所示， 0.78μs 時壓力峰值高達 25GPa ，部分區域內的炸藥已完全反應，轉化為爆轟產物。 0.90μs 時壓力峰值升至 27GPa ，在加載沖擊波以及新產生沖擊波的作用下，高壓區的剩余炸藥繼續劇烈反應，轉化為爆轟產物。 1.50μs 時 14mm 左側區域的炸藥幾乎全部發生反應。

圖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沖擊壓力下不同時刻的爆轟產物質量分數分布

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 方向運動速度分布。 0.27μs 時 7.50～7.67mm 處為能量匯聚的高密度、高壓力區，速度峰值約為 4.9km/s 。沖擊波波陣面往四周傳播， 0.78μs 時在5.50mm 處沖擊壓力已經達到 25GPa， ，高于約 14GPa 的CJ爆轟壓力，隨后壓力峰值一直高于CJ爆轟壓力。

通過數值模擬3種不同強度沖擊波加載下的液體炸藥硝基甲烷與氣泡相互作用過程，可以得到結論：在熱點即空穴塌縮處能夠實現化學反應的加速，并且化學反應的劇烈程度隨著沖擊波壓力的增加而增大。 4GPa 沖擊壓力加載下，熱點產生后爆轟產物的質量分數增長非常緩慢。而 6GPa 沖擊壓力加載下，炸藥內空穴塌縮形成熱點，能夠實現爆轟，引起化學反應區域的擴散。這與Turley等的實驗結果[36]一致，即沖擊波強度不夠高時，化學反應區域會停止增長，而當沖擊壓力超過5GPa時，在空穴塌縮處發生爆轟，爆轟產物區域會逐漸向四周擴大。

5結論

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

參考文獻：

[1] 彭亞晶，葉玉清.含能材料起爆過程“熱點”理論研究進展[J].化學通報，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]于明.固體炸藥爆轟與惰性介質相互作用的一種擴散界面模型[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，FEDKIWR，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]姚成寶，王宏亮，浦錫鋒，等.空中強爆炸沖擊波地面反射規律數值模擬研究[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]劉鐵鋼，許亮.模擬多介質界面問題的虛擬流體方法綜述[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，FENG 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

（責任編輯 張凌云）