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

尚海林,趙 鋒,王文強(qiáng),傅 華

（中國(guó)工程物理研究院流體物理研究所沖擊波物理與爆轟物理國(guó)防科技重點(diǎn)實(shí)驗(yàn)室,四川 綿陽(yáng) 621900）

沖擊波作用下熱點(diǎn)的形成和成長(zhǎng)是非均質(zhì)炸藥沖擊起爆的關(guān)鍵。研究炸藥在不同細(xì)觀結(jié)構(gòu)和加載條件下熱點(diǎn)形成的機(jī)制,對(duì)于炸藥的安全評(píng)估和應(yīng)用具有重要意義。

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

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

本文中用3維離散元程序?qū)MX基PBX塑料粘接炸藥、含不同形狀不同尺寸孔洞的HMX炸藥在沖擊作用下熱點(diǎn)的形成過(guò)程進(jìn)行數(shù)值模擬,探討非均質(zhì)炸藥在沖擊作用下的細(xì)觀響應(yīng)特性?？紤]到顆粒尺度單質(zhì)炸藥的化學(xué)反應(yīng)復(fù)雜,僅計(jì)算到熱點(diǎn)生成,不包括化學(xué)反應(yīng)。

1 離散元法簡(jiǎn)介

1.1 基本原理

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

1.2 相鄰元間的作用力

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

法向粘性作用力為

式中:Cvn為法向粘性系數(shù),對(duì)應(yīng)連續(xù)介質(zhì)力學(xué)上的粘性;vn,ij是在接觸點(diǎn)元i相對(duì)于j的法向速度。

1.3 能耗與溫升

粘性力引起能量耗散,轉(zhuǎn)化為熱量,導(dǎo)致離散元溫度升高。由于與相鄰元間的相互作用使元i在1個(gè)時(shí)間步長(zhǎng)內(nèi)的耗散能增量為

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

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

2 熱點(diǎn)的數(shù)值模擬結(jié)果

2.1 材料參數(shù)

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

表1 計(jì)算所用物理參數(shù)Table 1 Physical parameters for the calculation

2.2 塑料粘接炸藥中的熱點(diǎn)

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

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

2.3 含孔洞HMX炸藥中的熱點(diǎn)

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

圖1 PBX炸藥計(jì)算模型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 球形孔洞的塌縮

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

沖擊波到達(dá)孔洞下表面后,該表面的速度迅速上升到?jīng)_擊波后粒子速度的2倍,同時(shí)沖擊波在該自由面上反射1個(gè)稀疏波進(jìn)入波后炸藥中。隨著球形孔洞逐漸閉合,初始時(shí)刻位于孔洞周圍的炸藥在孔洞內(nèi)產(chǎn)生一定的會(huì)聚作用。此后經(jīng)過(guò)加速的高速自由面與孔洞對(duì)面的靜止界面發(fā)生碰撞,碰撞后自由面的速度逐漸下降,其動(dòng)能轉(zhuǎn)化為內(nèi)能加熱周圍的炸藥,使孔洞周圍炸藥的溫度進(jìn)一步上升,形成熱點(diǎn)。

比較2個(gè)不同半徑孔洞塌縮過(guò)程可以發(fā)現(xiàn),在同樣的樣品以及加載條件下,大尺寸孔洞塌縮后得到的熱點(diǎn)溫度高于小孔洞塌縮得到的熱點(diǎn)溫度。這是因?yàn)榇罂锥丛谒s過(guò)程中形成的會(huì)聚流動(dòng)比小孔洞更劇烈,在高速自由面與孔洞對(duì)面靜止界面碰撞后有更多的動(dòng)能轉(zhuǎn)化為內(nèi)能來(lái)加熱周圍炸藥。孫錦山等[2]用ZDE程序研究熱點(diǎn)形成機(jī)理的過(guò)程中也得到過(guò)同樣的結(jié)果。

圖4 沖擊波剛到達(dá)炸藥上自由面時(shí)炸藥的速度矢量場(chǎng)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 立方體孔洞的塌縮

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

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

圖6 沖擊波剛到達(dá)炸藥上自由面時(shí)炸藥的速度矢量場(chǎng)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種孔洞的比較

本文模擬結(jié)果也記錄了最初位于孔洞下表面附近粒子z方向平均速度隨時(shí)間的變化過(guò)程,如圖8所示?？梢钥吹?不管對(duì)于何種形狀何種尺寸的孔洞,當(dāng)沖擊波到達(dá)該自由面后,沖擊波速度迅速上升到波后粒子速度的2倍,然后受到旁側(cè)稀疏波等的影響,沖擊波速度逐漸降低。

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

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

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

3 結(jié)束語(yǔ)

本文中用3維離散元方法模擬了沖擊波作用下PBX炸藥和含孔洞的HMX炸藥中熱點(diǎn)的形成過(guò)程。模擬結(jié)果表明,PBX炸藥受沖擊作用后熱點(diǎn)集中在炸藥晶體與粘結(jié)劑的結(jié)合部位,且HMX晶體溫升低于粘結(jié)劑,晶體邊界溫升高于內(nèi)部。對(duì)于含孔洞的HMX炸藥,大尺寸孔洞塌縮形成的熱點(diǎn)溫度高于小尺寸孔洞塌縮形成的熱點(diǎn)溫度,這是因?yàn)榇罂锥丛谒s過(guò)程中形成的會(huì)聚流動(dòng)比小孔洞更劇烈,在高速自由面與孔洞對(duì)面靜止界面碰撞后有更多的動(dòng)能轉(zhuǎn)化為內(nèi)能來(lái)加熱周圍炸藥;球形孔洞塌縮形成的熱點(diǎn)溫度比立方體孔洞塌縮形成的熱點(diǎn)溫度高,這是由于球形孔洞有一定的聚心作用。

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

[1] 章冠人,陳大年.凝聚炸藥起爆動(dòng)力學(xué)[M].北京:國(guó)防工業(yè)出版社,1991:89-128.

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

[3] 孫承緯,衛(wèi)玉章,周之奎.應(yīng)用爆轟物理[M].北京:國(guó)防工業(yè)出版社,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] 劉凱欣,高凌天.離散元法研究的評(píng)述[J].力學(xué)進(jìn)展,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] 唐志平.三維離散元方法及其在沖擊動(dòng)力學(xué)中的應(yīng)用[J].中國(guó)科學(xué):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] 劉雪娜.三維點(diǎn)集Voronoi圖的算法實(shí)現(xiàn)[J].計(jì)算機(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.