沖擊作用下非均質(zhì)炸藥熱點形成的離散元方法

劉 超， 石藝娜， 梁仙紅

（北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100088）

沖擊作用下非均質(zhì)炸藥熱點形成的離散元方法

劉 超， 石藝娜， 梁仙紅

（北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100088）

用離散元方法模擬以HMX為基的塑料粘結(jié)炸藥與HMX顆粒炸藥在沖擊載荷作用下的細觀響應(yīng)過程.結(jié)果表明對于HMX顆粒炸藥，炸藥顆粒界面處的溫度及壓力遠高于顆粒內(nèi)部；對于HMX為基的塑料粘結(jié)炸藥，粘結(jié)劑處的溫度及壓力亦高于顆粒內(nèi)部；兩種沖擊響應(yīng)過程的對比分析表明，粘結(jié)劑有效地降低了炸藥顆粒邊界處的溫度及壓力.

離散元方法；熱點；沖擊加載；炸藥

0 引言

非均質(zhì)炸藥的熱點形成機制研究是當前爆轟物理研究的熱點和難點問題之一.熱點形成機制研究涉及材料物理、化學(xué)以及力學(xué)等多個學(xué)科，同時涵蓋宏、細、微觀多個尺度.已有研究表明非均質(zhì)炸藥中的雜質(zhì)、氣泡、孔隙、粘結(jié)劑和炸藥的顆粒尺度、粘性、熱容、導(dǎo)熱系數(shù)、活化能、反應(yīng)熱、密度以及加載速率、應(yīng)力狀態(tài)、環(huán)境溫度等因素都會對熱點的形成與演化過程產(chǎn)生重要影響.因此，盡管提出了多種熱點形成機制［1］，但由于炸藥細觀結(jié)構(gòu)的多樣性以及熱點形成與演化的復(fù)雜性，目前尚無法通過實驗或唯象理論研究對炸藥內(nèi)部熱點的形成與演化過程進行細致描述，因而細觀尺度數(shù)值模擬就成為研究此類問題的重要手段.

對非均質(zhì)炸藥熱點機制的細觀數(shù)值模擬已成為當前爆轟物理研究的前沿課題.Los Alamos國家實驗室的Menikoff利用二維歐拉程序先后對顆粒炸藥及含缺陷炸藥，在沖擊載荷作用下的熱點形成過程進行了模擬［2－5］；Livermore國家實驗室的Reaugh采用三維ALE（任意拉氏－歐拉有限差分）程序，模擬了不同尺度下受沖擊炸藥的力熱響應(yīng)及化學(xué)反應(yīng)［6－8］；Sandia國家實驗室的Baer用三維歐拉程序CTH模擬了多孔炸藥的沖擊響應(yīng)過程［9－10］；2007年Baer在沖擊波科學(xué)與技術(shù)參考叢書中專門用一章，對于沖擊載荷作用下非均質(zhì)炸藥的細觀模擬研究進行了總結(jié)［11］.

近年來，多尺度算法研究發(fā)展迅猛，但多集中于宏觀尺度與微觀尺度計算方法耦合的搭橋類算法研究，細觀尺度的計算方法相對較少.目前，研究多晶材料沖擊響應(yīng)的細觀尺度數(shù)值模擬方法，原則上可分為兩類，一類是以離散元為代表的粒子類方法，另一類是以連續(xù)介質(zhì)力學(xué)為基礎(chǔ)的傳統(tǒng)數(shù)值模擬方法.與傳統(tǒng)數(shù)值模擬方法相比離散元具有構(gòu)建非均質(zhì)模型方法易行，顆粒間取向分布特征特性表征便捷，處理沖擊載荷作用下炸藥顆粒間常見的大變形、摩擦、斷裂等問題方便，算法實現(xiàn)簡單等優(yōu)點.

離散元法是一種無網(wǎng)格方法，于20世紀70年代初由美國學(xué)者Cundall首先提出［12］，最初主要應(yīng)用于巖石力學(xué)、顆粒態(tài)群體及土壤力學(xué)問題分析中.20世紀90年代初，日本學(xué)者Sawamoto［13］首先將離散元方法成功地用于混凝土動態(tài)沖擊破壞等非線性大變形問題的數(shù)值模擬研究.北京大學(xué)的劉凱欣教授，在這一領(lǐng)域做了大量的研究工作［14－15］.1999年以來，Horie等利用離散元方法研究了銅、鐵等多晶金屬的沖擊響應(yīng)［16－17］. 2000年以后，中國科學(xué)技術(shù)大學(xué)的唐志平、王文強將離散元方法應(yīng)用到非均質(zhì)材料炸藥在沖擊作用下細觀變形及損傷研究中［18］.中國工程物理研究院流體物理研究所的傅華、趙峰等利用該方法研究了HMX晶體在沖擊載荷作用下，由于空洞塌縮形成的熱點［19］.這些工作展示了離散元方法模擬細觀非均質(zhì)材料沖擊動力學(xué)問題的能力.由于可以方便地表征顆粒間的取向分布特征，描述顆粒間的相互作用，離散元方法在模擬細觀非均質(zhì)材料的沖擊響應(yīng)方面具有獨特的優(yōu)勢.

本文利用離散元方法對以HMX為基的塑料粘結(jié)炸藥與HMX顆粒炸藥沖擊作用下的熱點形成與演化過程進行模擬，探討非均質(zhì)炸藥在沖擊載荷作用下的細觀響應(yīng)特征，分析沖擊過程中塑料粘結(jié)劑所起的作用.計算中僅考慮未反應(yīng)炸藥在沖擊作用下的力、熱響應(yīng)，未考慮固固、固液相變及化學(xué)反應(yīng)等復(fù)雜因素.

1 計算模型與方法

離散元方法通過求解多體運動的牛頓力學(xué)方程組，跟蹤全部個體（即單元）的運動軌跡，來揭示系統(tǒng)與外界的相互作用和自身響應(yīng)、演化規(guī)律.與傳統(tǒng)數(shù)值模擬方法相比離散元具有非均質(zhì)材料模型構(gòu)建方法易行，顆粒間取向分布特征表征便捷，處理沖擊載荷作用下炸藥顆粒間常見的大變形、摩擦、斷裂等問題方便，算法實現(xiàn)簡單等優(yōu)點.

1.1 單元間相互作用力模型

離散元方法的核心是根據(jù)所研究的問題，選取合適的單元間相互作用力模型.一般單元間可能包括的相互作用力，見圖1［20］.

圖1 DM2中單元間的相互作用力的示意圖Fig.1 Schematic of interaction forces between elemonts in DM2

圖1（a）中單元間的相互作用力包括：①中心勢力，②中心阻尼，③彈塑性剪切力，④切向阻尼，⑤干摩擦力.圖中V為線速度，ω為角速度.圖1（b）為由中心勢力構(gòu)成的應(yīng)力（σcp）與單元間距（d）的關(guān)系；當0＜d ＜ds時，材料的響應(yīng)由Lennard-Johns勢函數(shù)描述；當d＝ds時，材料發(fā)生屈服，屈服強度與溫度有關(guān)；當d≥dmax時單元間將失去連接.如圖1所示，單元i與單元j間的相互作用力合力可表示為

表1中，αij、m、n為中心勢力參數(shù)，Cn為阻力系數(shù)，Cd為摩擦系數(shù).

表1 單元間作用力模型參數(shù)Table1 Interaction force parameters

1.2 單元溫度的計算

將影響單元溫度的力學(xué)過程分為可逆過程與不可逆過程，可逆的力學(xué)過程如中心勢力的作用過程，不可逆的力學(xué)過程為耗散力的作用過程

溫度的可逆部分可表示為

其中，T0＝300 K為初始溫度，V0為初始比體積，V為比體積，Γ為Gruneisen系數(shù).

不可逆部分考慮了由熱傳導(dǎo)、粘性力和炸藥顆粒間的摩擦力帶來的能量耗散過程，不可逆過程帶來的溫升是上述過程的累計效應(yīng)：

接觸單元間考慮了基于傅立葉定律的熱傳導(dǎo)過程

ΔQ為在Δt時間內(nèi)由單元j向單元i傳導(dǎo)的熱，κ為熱傳導(dǎo)系數(shù)，Ti和Tj分別為單元i與單元j的溫度，dij為單元i與單元j間的距離，Aij為單元i與單元j間的接觸面積.

由熱傳導(dǎo)導(dǎo)致的單元溫度變化，表示為

由粘性力帶來的溫升

其中，M為單元質(zhì)量，cV為等容比熱，Cn為粘性系數(shù)，vijn為單元間法向速度.

炸藥顆粒間的摩擦力僅存在于炸藥顆粒邊界處的單元內(nèi)，炸藥顆粒內(nèi)部的單元間不考慮摩擦力.由摩擦力帶來的溫升表示為

表2 熱力學(xué)參數(shù)Table 2 Thermodynamic parameters

其中，Cd為摩擦系數(shù)，Nij為單元間的正壓力，vijt為單元間切向速度.單元溫度計算中熱力學(xué)參數(shù)的取值參見表2.

2 均質(zhì)HMX炸藥的數(shù)值模擬及結(jié)果分析

計算模型中飛片與靶板分別取300 μm×1.2 μm與900 μm×1.2 μm的均質(zhì)HMX炸藥，單元直徑為0.4 μm，計算單元總數(shù)10 392個.如圖2所示，計算模型的上、下邊界采用周期性邊界條件，左、右邊界采用自由邊界條件.模型初始溫度300 K，初始時刻飛片不同的速度從左端撞擊靶板.

圖3給出了HMX的P-V／V0Hugonoit曲線，圖中實心方塊為Yoo等人的實驗數(shù)據(jù)［21］，實心三角為LASL實驗室的實驗數(shù)據(jù)［22］，實線為本文計算結(jié)果.圖4中點劃線為Menikoff［3］等利用Hayes狀態(tài)方程得到的結(jié)果，實線為本文采用上述溫度計算模型得到的計算結(jié)果.計算結(jié)果與文獻中的實驗結(jié)果或數(shù)值模擬結(jié)果符合較好，證明計算所用單元間作用力模型及溫度計算模型能夠正確反應(yīng)HMX的沖擊壓縮特性.

圖2 計算模型示意圖Fig.2 Initial model sketch

3 非均質(zhì)HMX炸藥的數(shù)值模擬及結(jié)果分析

從細觀尺度看，固體炸藥是典型的非均質(zhì)材料，其細觀結(jié)構(gòu)特征將顯著影響其動態(tài)響應(yīng)過程：王金英等人在落錘實驗中發(fā)現(xiàn)，粘結(jié)劑可以有效降低炸藥感度［23］.傳統(tǒng)的數(shù)值模擬方法難以表征炸藥顆粒間的取向分布特征，及其在沖擊載荷作用下可能出現(xiàn)的大變形、摩擦、斷裂等問題，因而在模擬細觀非均質(zhì)材料響應(yīng)方面遇到較大困難，而離散元方法在模擬此類問題時具有獨特的優(yōu)勢.本文采用離散元方法對于以HMX為基的塑料粘結(jié)炸藥與HMX顆粒炸藥，在沖擊作用下的熱點形成與演化過程進行了數(shù)值模擬研究.

圖3 HMX的P-V／V0Hugoniot關(guān)系Fig.3 P-V／V0Hugonoit of HMX

圖4 HMX的P-T關(guān)系Fig.4 P-T relation of HMX

為避免邊界效應(yīng)的影響，計算模型的左、右邊界為周期型邊界條件，下邊界為固壁邊界，上邊界為自由邊界，初始時刻計算模型（見圖5）以1 km·s－1的速度撞擊固壁.

圖5 HMX顆粒炸藥與以HMX為基的PBX炸藥初始計算模型Fig.5 Initial model of granule HMX and PBX

圖5中不同顏色的顆粒代表顆粒間不同的晶粒取向，計算模型中的每個炸藥顆粒的晶粒取向呈隨機分布，以此表征炸藥顆粒間的取向分布特性.圖6為計算模型的放大圖，可見顆粒內(nèi)部單元為規(guī)則的密排六邊形，顆粒處的間隙由小尺度單元填滿，各顆粒的晶粒取向為（0°～60°）之間的隨機數(shù).圖5與圖6中PBX炸藥計算模型的顆粒邊界處可見黑色的單元，代表塑料粘接劑Estane.

圖6 HMX顆粒炸藥與以HMX為基的PBX炸藥模型局部放大圖Fig.6 An enlarged view of granule HMX and PBX

圖7為HMX顆粒炸藥與以HMX為基的PBX炸藥模型分別以1 km·s－1速度撞擊固壁后，某時刻流場內(nèi)的壓力分布圖.圖8著重顯示了同一時刻，兩種炸藥內(nèi)部壓力高于10 GPa的區(qū)域.可以看到由于炸藥顆粒間取向差異及顆粒邊界效應(yīng)的影響，壓力場的分布很不均勻，應(yīng)力遠高于或低于波后平均應(yīng)力的區(qū)域集中于顆粒邊界附近.并且雖然構(gòu)建的PBX炸藥模型中粘結(jié)劑的組分含量很低，其質(zhì)量百分比約4%，但其緩沖作用顯著，從圖8可見以HMX為基的PBX炸藥受沖擊后的較高壓力區(qū)明顯少于HMX顆粒炸藥.

圖7 兩種炸藥受沖擊后的壓力分布Fig.7 Pressure fields in two models

圖8 兩種炸藥受沖擊后的高壓位置分布Fig.8 High pressure positions in two models

圖9為兩種炸藥同一位置的壓力放大圖，圖9（a）為HMX顆粒炸藥的顆粒取向與壓力放大圖，圖9（b）為PBX炸藥的物質(zhì)與壓力放大圖.可見應(yīng)力遠高于或低于波后平均應(yīng)力的區(qū)域集中于粘結(jié)劑或炸藥顆粒邊界附近，并且PBX炸藥粘結(jié)劑處的壓力明顯低于HMX炸藥顆粒邊界處的壓力.這是由于粘結(jié)劑較軟，起到一定的緩沖作用，有效地降低了PBX炸藥顆粒邊界處的壓力.

圖10為HMX顆粒炸藥與以HMX為基的PBX炸藥模型分別以1 km·s－1速度撞擊固壁后，某時刻波后流場內(nèi)的溫度分布圖.圖11著重顯示了同一時刻，兩種炸藥內(nèi)部溫度高于700 K的區(qū)域.可以看到波后的溫度場的分布很不均勻，由于炸藥顆粒間的取向差異及顆粒邊界的影響，溫度遠高于或低于波后平均溫度的區(qū)域主要集中于炸藥顆粒邊界或粘結(jié)劑區(qū)域附近.盡管PBX模型中粘結(jié)劑所占百分比很少，但由于粘結(jié)劑較軟易于變形，并在變形過程中吸收沖擊能量，使得PBX炸藥的熱點溫度明顯低于顆粒炸藥，如圖11.

圖9 兩種炸藥受沖擊后的壓力放大圖Fig.9 Amplified pressure fields in two models

圖10 兩種炸藥受沖擊后的溫度分布Fig.10 Temperature fields in two models

圖11 兩種炸藥受沖擊后高溫位置分布Fig.11 High temperature positions in two models

4 結(jié)論

用離散元方法模擬了HMX炸藥的沖擊響應(yīng)過程，得到了HMX的沖擊Hugoniot關(guān)系與P-T關(guān)系，計算結(jié)果與文獻中的實驗結(jié)果及數(shù)值模擬結(jié)果符合較好，證明計算所用單元間作用力模型及溫度計算模型能夠正確反應(yīng)HMX的沖擊壓縮特性.對于以HMX為基的塑料粘結(jié)炸藥與HMX顆粒炸藥在沖擊載荷作用下的細觀響應(yīng)過程進行了數(shù)值模擬，結(jié)果表明對于HMX顆粒炸藥，炸藥顆粒界面處的溫度及壓力明顯高于顆粒內(nèi)部；對于HMX為基的PBX炸藥，粘結(jié)劑處的溫度及壓力亦高于顆粒內(nèi)部；對于上述兩種沖擊響應(yīng)過程的比較表明粘結(jié)劑較軟，能起到一定的緩沖作用，并在變形過程中吸收沖擊能量，從而有效的降低了炸藥顆粒邊界處的溫度與壓力.

［1］Field J E.Hot spot ignition mechanisms for explosives［J］.Acc Chenm Res，1992，25：489－496.

［2］Menikoff R.Pore collapse and hot spots in HMX［C］∥Conferences on Shock Compression of Condensed Matter，2004，706：393－396.

［3］Menikoff R，Sewell T D.Constituent properties of HMX needed for mesoscale simulations［J］.Combust Theory Model，2002，6：103－125.

［4］Menikoff R.Granular explosives and initiation sensitivity［C］∥Conferences on Shock Compression of Condensed Matter，2002，979－983.

［5］Menikoff R.Hot spot formation from shock reflections［J］.Shock Wave，2011，21：141－148.

［6］Reaugh J E.Grain-scale dynamics in explosive［R］.Lawrence Livermore National Laboratory，UCRL-ID－150388.2002.

［7］Reaugh J E.First results of reaction propagation rates in HMX at high pressure［C］∥Conferences on Shock Compression of Condensed Matter，2003：401－404.

［8］Reaugh J E.Multi-scale computer simulations to study the reaction zone of solid explosives［C］∥Proceedings of the 13th Int Symp on Detonation，2006：1276－1285.

［9］Baer M R，Trott W M.Mesoscale descriptions of shock-loaded heterogeneous porous materials［C］∥Conferences on Shock Compression of Condensed Matter，2001：713－716.

［10］Baer M R.Modeling heterogeneous energetie materials at the mesoseale［J］.Thermochimic Acta，2002，384：351－367.

［11］Baer M R，Trott W M.Mesoscale studies of shock loaded tin spheres［C］∥Conferences on Shock Compression of Condensed Matter，2003：517－520.

［12］Cundall P A.A computer model for simulating progressive large scale movement in blockrock system［J］.Symposium ISRM，1971，（2）：129－136.

［13］Sawamoto Y，Tsubota H，Kasai Y，Koshika N，Morikawa H.Analytical studies on local damage to reinforced concrete structures under impact loading by discrete element method［J］.Nucl Eng Des，1998，179：157－177.

［14］劉凱欣，高凌天，鄭文剛.混凝土動態(tài)破壞過程的數(shù)值模擬［J］.工程力學(xué)（增刊），2000：470－474.

［15］劉凱欣，高凌天.離散元法在求解三維沖擊動力學(xué)問題中的應(yīng)用［J］.固體力學(xué)學(xué)，2004，25（2）：181－185.

［16］Yano K，Horie Y.Discrete-element modeling of shock compression of polycrystalline copper［J］.Phys Rev B，1999，59：13672－13680.

［17］Case S，Horie Y.Mesomechanics of the α-d transition in Iron［J］.J Mech Phys Solids，2007，55：589－614.

［18］王文強.離散元方法及其在材料和結(jié)構(gòu)力學(xué)響應(yīng)分析中的應(yīng)用［D］.合肥：中國科學(xué)技術(shù)大學(xué)，2000.

［19］傅華，趙峰，譚多望，王文強，尚海林.沖擊作用下HMX晶體孔洞塌縮熱點生成機制的細觀數(shù)值模擬［J］.高壓物理學(xué)報，2011，25（1）：8－14.

［20］Tang Z P，Horie Y，Psakhie.Diserete meso-element modeling of shock processes in powders［M］∥Graham R A，et al.Shock Compression of Solids，IV.Springer-Verlag，1997：143－176.

［21］Yoo C S，Cynn H.Equation of state，phase transition，decomposition of β-HMX［J］.J Chem Phys，1999，111：10229－10235.

［22］Marsh S.LASL shock Hugoniot data［M／OL］.Univ Calif Press，1980，http：／／lib-www.lanl.gov／books／shd.pdf.

［23］王金英，柴濤，張景林.PBX傳爆藥撞擊感度影響因素的研究［J］.華北工學(xué)院學(xué)報，2004，25（4）：289－292.

DEM Study on Hot Spots Formation of Heterogeneous Explosives Under Shock Loading

LIU Chao，SHI Yina，LIANG Xianhong
（Institute of Applied Physics and Computational Mathematics，Beijing 100088，China）

Mesoscale responses of plastic boned explosives（PBX）and grannular HMX under shock loading were investigated with discrete element method（DEM）.It shows that for shocked grannular HMX temperature and pressure in grain boundary are much higher than those of inner HMX crystals.In the case of shocked PBX high temperature and pressure regions are mostly located in binder and near parts.Comparisons show that temperature and pressure in grain boundary are effectively reduced by binder.

discrete element method；hot spot；shock loading；explosive

date：2013－10－28；Revised date：2014－03－20

O347.1

A

2013－10－28；

2014－03－20

國家自然科學(xué)基金（11202034）；中物院科學(xué)技術(shù)發(fā)展基金（2011B0101028，2013B0101013）及中物院科學(xué)技術(shù)重點基金（2012A0101004）資助項目

劉超（1976－），女，碩士，副研究員，主要從事計算爆炸力學(xué)研究，E-mail：liu_chao＠iapcm.ac.cn

1001-246X（2014）05-0523-07