沖擊作用下炸藥熱點形成的3維離散元模擬*

尚海林,趙 鋒,王文強,傅 華

（中國工程物理研究院流體物理研究所沖擊波物理與爆轟物理國防科技重點實驗室,四川 綿陽 621900）

沖擊波作用下熱點的形成和成長是非均質炸藥沖擊起爆的關鍵。研究炸藥在不同細觀結構和加載條件下熱點形成的機制,對于炸藥的安全評估和應用具有重要意義。

熱點是細觀尺度的問題,在炸藥的多尺度響應中處于承上（宏觀）啟下（微觀）的地位。在各種情況下,熱點形成的機制主要包括:氣泡的絕熱壓縮、粘塑性孔洞塌縮,剪切帶、晶體中位錯的運動,炸藥顆粒間摩擦以及流體動力學熱點等[1-3]。由于實驗上難以對熱點進行直接測量,目前對這一問題的研究還只能依靠數值模擬和理論分析。

近年來,研究人員對炸藥熱點形成機制進行了一系列的數值模擬研究。R.Menikoff等[4]用2維歐拉有限元程序模擬了顆粒HMX炸藥中的壓實波,結果表明在500 m/s活塞推動下熱點溫度可達到600 K,R.Menikoff[5]還研究了炸藥中氣泡塌縮與熱點的關系;J.E.Reaugh[6]用ALE-3D（3維彈塑性流體動力學任意拉格朗日-歐拉有限差分程序）在不同特征尺度下模擬了炸藥在沖擊加載和其他壓力加載情況下的力學、熱學和化學響應;M.R.Baer等[7-8]用3維歐拉有限體積程序CTH,以離散的HMX晶粒受碰撞為例研究了沖擊條件下多孔含能材料的細觀壓實、變形和反應過程,數值模擬結果表明熱點強烈依賴于多晶體之間的相互作用。許愛國等[9]用物質點法研究了沖擊作用下固體中孔洞的塌縮,分析了孔洞塌縮的對稱性和沖擊強度的關系,結果表明材料中熱點的分布很大程度上受塌縮歷史的影響。

本文中用3維離散元程序對HMX基PBX塑料粘接炸藥、含不同形狀不同尺寸孔洞的HMX炸藥在沖擊作用下熱點的形成過程進行數值模擬,探討非均質炸藥在沖擊作用下的細觀響應特性。考慮到顆粒尺度單質炸藥的化學反應復雜,僅計算到熱點生成,不包括化學反應。

1 離散元法簡介

1.1 基本原理

離散元法[10]把不連續體分離為剛性元素的集合,使各個剛性元素滿足運動方程,用時步迭代法求解各剛性元素的運動方程,繼而求得不連續體的整體運動形態。離散元法允許單元間的相對運動,不一定要滿足位移連續和變形協調條件,計算速度快,所需存儲空間小,適合求解大位移和非線性問題。研究連續介質時,離散元法是將介質直接離散為元,根據材料性質確定元間作用力模型,通過求解元的控制方程（牛頓第二定律）獲知整個系統的演化規律[11]。本文離散元模型以3維離散元方法[12-13]為基礎。

1.2 相鄰元間的作用力

相鄰元間的作用力包括中心勢作用力和法向粘性作用力。中心勢作用力表達式如下

法向粘性作用力為

式中:Cvn為法向粘性系數,對應連續介質力學上的粘性;vn,ij是在接觸點元i相對于j的法向速度。

1.3 能耗與溫升

粘性力引起能量耗散,轉化為熱量,導致離散元溫度升高。由于與相鄰元間的相互作用使元i在1個時間步長內的耗散能增量為

式中:θvn是元i得到的能耗的百分比。假設粘性耗散能90%轉化為熱,元i和其相鄰元j各得到一半。則有θvn=0.45。因此,在絕熱條件下,元i在1個時步內的溫升為

式中:cV,i為元i的定容熱容,mi為元i的質量。

2 熱點的數值模擬結果

2.1 材料參數

計算過程中用到的物理參數見表1,其中HMX炸藥和粘結劑間的相互作用參數的確定方法如下:兩者相互作用某一參數的倒數等于2種材料對應倒數和的平均值,如對于參數λ,設HMX的為λ1,粘結劑的為λ2,則HMX與粘結劑間的對應參數λ滿足關系1/λ=（1/λ1+1/λ2）/2;粘性系數Cvn與壓力、溫度、應變率等因素有關,并不是一個常數,文獻[15]中給出的粘性系數Cvn=10～100 Pa·s,本文計算中Cvn=11 Pa·s。

表1 計算所用物理參數Table 1 Physical parameters for the calculation

2.2 塑料粘接炸藥中的熱點

本文中基于3維Voronoi拼圖方法[16]建立能夠基本反映炸藥細觀構型的幾何模型,如圖1所示,其中淺色離散元代表HMX炸藥,深色離散元代表粘結劑。模型為邊長0.4 mm的立方體炸藥樣品,含有64個HMX炸藥顆粒、77 904個HMX離散元（離散元半徑為5 μm）、123 072個粘結劑離散元（離散元半徑為2 μm）,離散元按照類似于晶格點陣中立方密堆積的方式進行排列,HMX炸藥的質量分數為94.5%,粘結劑的質量分數為5.5%,接近于PBX9501炸藥。模型在x和y方向上的邊界皆為滑移固壁,即離散元到達邊界的時候只允許上下運動,不能水平運動;上邊界為自由面;下邊界為以500 m/s的速度向上運動的活塞。樣品初始溫度為293 K,計算時間步長為0.2 ns,其余參數見表1。

為了更清楚地看見炸藥內部溫度的變化和熱點的形成過程,在垂直于z方向1/3樣品厚度處選取1個截面,圖2為該截面在沖擊波到達后2個不同時刻PBX炸藥的溫度分布圖,圖3為該截面HMX晶粒的溫度分布圖。模擬結果清楚地顯示了隨著沖擊波自下而上傳播,晶體與粘結劑間的相互作用使炸藥局部溫度不斷升高并且形成熱點的過程。計算得到的熱點溫度為幾百K,與 R.Menikoff等[4]給出的結果基本符合。從計算結果中可以看出熱點集中在炸藥晶體與粘結劑的結合部位,且HMX晶體溫升明顯低于粘結劑,晶體邊界溫升高于內部,這與P.A.Conley等[14]的模擬結果一致。

2.3 含孔洞HMX炸藥中的熱點

炸藥樣品為邊長0.4 mm的立方體,炸藥為HMX單晶炸藥,在樣品中心處挖掉1個孔洞,孔洞內部為真空,邊界條件、初始條件以及計算所用物理參數與上節一樣。為了確定孔洞形狀和尺寸對塌縮過程以及塌縮后形成熱點溫度的影響,分別選擇了球形和立方體孔洞進行計算,球形孔洞直徑取0.08和0.12 mm,立方體孔洞邊長取0.08和0.12 mm。

圖1 PBX炸藥計算模型Fig.1 A calculation model for the PBX explosive

圖2 z平面PBX炸藥溫度分布Fig.2 Temperature distribution of PBX explosive in the z-plane

圖3 z平面HMX晶體溫度分布Fig.3 Temperature distribution of HMX crystal in the z-plane

2.3.1 球形孔洞的塌縮

對于含球形孔洞的炸藥,模擬結果記錄了孔洞塌縮過程中炸藥的速度場以及孔洞周圍炸藥的溫度變化,如圖4～5所示,3維模型中看不到內部孔洞,圖中顯示的是經過孔洞中心垂直于x方向具有一定厚度的薄片。

沖擊波到達孔洞下表面后,該表面的速度迅速上升到沖擊波后粒子速度的2倍,同時沖擊波在該自由面上反射1個稀疏波進入波后炸藥中。隨著球形孔洞逐漸閉合,初始時刻位于孔洞周圍的炸藥在孔洞內產生一定的會聚作用。此后經過加速的高速自由面與孔洞對面的靜止界面發生碰撞,碰撞后自由面的速度逐漸下降,其動能轉化為內能加熱周圍的炸藥,使孔洞周圍炸藥的溫度進一步上升,形成熱點。

比較2個不同半徑孔洞塌縮過程可以發現,在同樣的樣品以及加載條件下,大尺寸孔洞塌縮后得到的熱點溫度高于小孔洞塌縮得到的熱點溫度。這是因為大孔洞在塌縮過程中形成的會聚流動比小孔洞更劇烈,在高速自由面與孔洞對面靜止界面碰撞后有更多的動能轉化為內能來加熱周圍炸藥。孫錦山等[2]用ZDE程序研究熱點形成機理的過程中也得到過同樣的結果。

圖4 沖擊波剛到達炸藥上自由面時炸藥的速度矢量場Fig.4 Velocity field of the explosive upon the shock wave reaching the free surface

圖5 球形孔洞完全塌縮后炸藥整體的溫度分布Fig.5 Temperature distribution of explosive after the collapse of a spherical void

2.3.2 立方體孔洞的塌縮

對于含立方體孔洞的炸藥,模擬結果同樣記錄了孔洞塌縮過程中炸藥的速度場以及孔洞周圍炸藥的溫度變化,如圖6～7所示。3維模型中看不到內部孔洞,圖中顯示的是經過孔洞中心垂直于x方向具有一定厚度的薄片。

對于立方體形狀的孔洞,除了孔洞下表面受沖擊后高速向上運動并與上表面發生碰撞將動能轉化為內能外,沖擊壓力還會將孔洞的側面擠壓變形,這些被擠壓變形的面達到斷裂極限后發生斷裂,動量和能量都聚到孔洞中心加速粒子,被加速的粒子與孔洞對面靜止界面產生高速碰撞,形成熱點。同樣,圖7顯示大孔洞塌縮后的熱點溫度更高。H.Takahiro[17]曾經用3維分子動力學方法模擬了炸藥晶體中孔洞的塌縮過程,得到同樣的結果,他的研究結果也顯示熱點溫度隨著縱向（沿沖擊波傳播的方向）尺寸的增大而增加,若固定縱向尺寸而變化橫向尺寸,熱點溫度會在橫向尺寸與縱向尺寸相等即孔洞為立方體形狀的情況下達到最高。

圖6 沖擊波剛到達炸藥上自由面時炸藥的速度矢量場Fig.6 Velocity field of the explosive upon the shock wave reaching the free surface

圖7 立方體孔洞完全塌縮后炸藥整體的溫度分布Fig.7 Temperature distribution of explosive after the collapse of a cubic void

2.3.3 2種孔洞的比較

本文模擬結果也記錄了最初位于孔洞下表面附近粒子z方向平均速度隨時間的變化過程,如圖8所示。可以看到,不管對于何種形狀何種尺寸的孔洞,當沖擊波到達該自由面后,沖擊波速度迅速上升到波后粒子速度的2倍,然后受到旁側稀疏波等的影響,沖擊波速度逐漸降低。

圖9顯示孔洞周圍一薄層炸藥在整個過程中平均溫度隨時間變化的曲線。從圖中可以看到,對于同種形狀的孔洞,其溫度曲線相互交叉,完全塌縮之前小孔洞周圍炸藥平均溫度較高,這是因為小孔洞塌縮所需要的時間較短,溫度上升快,完全塌縮之后大孔洞周圍平均溫度高;另外,比較同樣尺寸的球形空洞和立方體孔洞塌縮的模擬結果可以看到,球形孔洞塌縮最終形成的熱點溫度更高,這是由于球形孔洞有一定的聚心作用而引起的[2]。

圖8 孔洞下表面z向平均速度隨時間的變化Fig.8 Average z-velocity evolution on the lower wall of the void

圖9 孔洞周圍炸藥平均溫度隨時間的變化Fig.9 Average temperature evolution of the explosive around the void

3 結束語

本文中用3維離散元方法模擬了沖擊波作用下PBX炸藥和含孔洞的HMX炸藥中熱點的形成過程。模擬結果表明,PBX炸藥受沖擊作用后熱點集中在炸藥晶體與粘結劑的結合部位,且HMX晶體溫升低于粘結劑,晶體邊界溫升高于內部。對于含孔洞的HMX炸藥,大尺寸孔洞塌縮形成的熱點溫度高于小尺寸孔洞塌縮形成的熱點溫度,這是因為大孔洞在塌縮過程中形成的會聚流動比小孔洞更劇烈,在高速自由面與孔洞對面靜止界面碰撞后有更多的動能轉化為內能來加熱周圍炸藥;球形孔洞塌縮形成的熱點溫度比立方體孔洞塌縮形成的熱點溫度高,這是由于球形孔洞有一定的聚心作用。

利用離散元方法可以方便地建立炸藥復雜的細觀模型,處理大變形和破碎等問題。但是由于離散元方法的離散特性,該方法在定量描述剪切變形及塑性等現象時不夠理想,因而本文中只考慮了離散元間的粘彈性法向作用力,這適用于流體力學行為起主要作用的較高沖擊壓力加載情況,對于低壓情況則有較大偏差。有必要進一步研究離散元和有限元相結合的方法,利用2種計算方法的優勢,更好地模擬在低沖擊壓力作用下粘塑性對炸藥熱點形成的影響。

[1] 章冠人,陳大年.凝聚炸藥起爆動力學[M].北京:國防工業出版社,1991:89-128.

[2] 孫錦山,朱建士.理論爆轟物理[M].北京:國防工業出版社,1995:327-341.

[3] 孫承緯,衛玉章,周之奎.應用爆轟物理[M].北京:國防工業出版社,2000:430-433.

[4] Menikoff R,Kober E.Compaction waves in granular HM X[C]//Furnish M D,Chhabildas L C,Hixson R S.Proceedings of the APS Topical Conference on Shock Compression of Condensed Matter.Snowbird,CO,1999.

[5] Menikoff R.Pore collapse and hot spots in HMX[C]//Proceedings of the Conference of the American Physical Society Topical Group on Shock Compression of Condensed M atter.Portland,Oregon,United States:AIP Press,2004,706:393-396.

[6] Reaugh J E.M ulti-scale computer simulations to study the reaction zone of solid explosives[C]//Proceedings of the 13th International Detonation Symposium.Norfolk,VA,United States,2006.

[7] Baer M R.Modeling heterogeneous energetic materials at the mesoscale[J].Thermochimica Acta,2002,384:351-367.

[8] Baer M R,Kipp M E,Van Swol F.Micromechanical modeling of heterogeneous energetic materials[C]//Proceedings of the 11th International Detonation Symposium.Snowmass,CO,United States,1998:788-797.

[9] Xu A G,Pan X F,Zhang G C.M aterial-point simulation of cavity collapse under shock[J].Journal of Physics:Condensed Matter,2007,19:1-11.

[10] Cundall P A.A computer model for simulating progressive,large-scale movements in blocky rock systems[C]//Proceedings of ISRM Symposium on Rock Fracture.Nancy,1971:11-18.

[11] 劉凱欣,高凌天.離散元法研究的評述[J].力學進展,2003,33（4）:483-490.

LIU Kai-xin,GAO Ling-tian.A review on the discrete element method[J].Advances in Mechanics,2003,33（4）:483-490.

[12] Tang Z P,Horie Y,Psakhic S G.Discrete meso-element modeling of shock processes in porous materials:Theory and model calculations[R].Technical Report,North Carolina State University,1995.

[13] 唐志平.三維離散元方法及其在沖擊動力學中的應用[J].中國科學:E輯,2003,33（11）:989-998.

TANG Zhi-ping.3D discrete element method and its application in impact dynamics[J].Science in China:E,2003,33（11）:989-998.

[14] Conley P A,Benson D J.Microstructural effects in shock initiation[C]//Proceedings of the 11th International Detonation Symposium.Snowmass,CO,United States,1998:768-780.

[15] M enikoff R,Sewell T D.Constituent properties of HMX needed for mesoscale simulations[J].Combustion Theory and M odeling,2002,6:103-125.

[16] 劉雪娜.三維點集Voronoi圖的算法實現[J].計算機輔助工程,2006,15（1）:1-4.

LIU Xue-na.Implementation of Voronoi diagram algorithms for 3D point set[J].Computer Aided Engineering,2006,15（1）:1-4.

[17] Takahiro H.Spatiotemporal behavior of void collapse in shocked solids[J].Physical Review Letters,2004,92（1）:015503-1-4.