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

劉 超， 石藝娜， 梁仙紅

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

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

劉 超， 石藝娜， 梁仙紅

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

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

離散元方法；熱點(diǎn)；沖擊加載；炸藥

0 引言

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

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

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

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

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

1 計(jì)算模型與方法

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

1.1 單元間相互作用力模型

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

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

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

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

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

1.2 單元溫度的計(jì)算

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

溫度的可逆部分可表示為

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

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

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

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

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

由粘性力帶來(lái)的溫升

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

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

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

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

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

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

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

圖2 計(jì)算模型示意圖Fig.2 Initial model sketch

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

從細(xì)觀尺度看，固體炸藥是典型的非均質(zhì)材料，其細(xì)觀結(jié)構(gòu)特征將顯著影響其動(dòng)態(tài)響應(yīng)過(guò)程：王金英等人在落錘實(shí)驗(yàn)中發(fā)現(xiàn)，粘結(jié)劑可以有效降低炸藥感度［23］.傳統(tǒng)的數(shù)值模擬方法難以表征炸藥顆粒間的取向分布特征，及其在沖擊載荷作用下可能出現(xiàn)的大變形、摩擦、斷裂等問(wèn)題，因而在模擬細(xì)觀非均質(zhì)材料響應(yīng)方面遇到較大困難，而離散元方法在模擬此類問(wèn)題時(shí)具有獨(dú)特的優(yōu)勢(shì).本文采用離散元方法對(duì)于以HMX為基的塑料粘結(jié)炸藥與HMX顆粒炸藥，在沖擊作用下的熱點(diǎn)形成與演化過(guò)程進(jìn)行了數(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)的影響，計(jì)算模型的左、右邊界為周期型邊界條件，下邊界為固壁邊界，上邊界為自由邊界，初始時(shí)刻計(jì)算模型（見(jiàn)圖5）以1 km·s－1的速度撞擊固壁.

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

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

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

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

圖10為HMX顆粒炸藥與以HMX為基的PBX炸藥模型分別以1 km·s－1速度撞擊固壁后，某時(shí)刻波后流場(chǎng)內(nèi)的溫度分布圖.圖11著重顯示了同一時(shí)刻，兩種炸藥內(nèi)部溫度高于700 K的區(qū)域.可以看到波后的溫度場(chǎng)的分布很不均勻，由于炸藥顆粒間的取向差異及顆粒邊界的影響，溫度遠(yuǎn)高于或低于波后平均溫度的區(qū)域主要集中于炸藥顆粒邊界或粘結(jié)劑區(qū)域附近.盡管PBX模型中粘結(jié)劑所占百分比很少，但由于粘結(jié)劑較軟易于變形，并在變形過(guò)程中吸收沖擊能量，使得PBX炸藥的熱點(diǎn)溫度明顯低于顆粒炸藥，如圖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)過(guò)程，得到了HMX的沖擊Hugoniot關(guān)系與P-T關(guān)系，計(jì)算結(jié)果與文獻(xiàn)中的實(shí)驗(yàn)結(jié)果及數(shù)值模擬結(jié)果符合較好，證明計(jì)算所用單元間作用力模型及溫度計(jì)算模型能夠正確反應(yīng)HMX的沖擊壓縮特性.對(duì)于以HMX為基的塑料粘結(jié)炸藥與HMX顆粒炸藥在沖擊載荷作用下的細(xì)觀響應(yīng)過(guò)程進(jìn)行了數(shù)值模擬，結(jié)果表明對(duì)于HMX顆粒炸藥，炸藥顆粒界面處的溫度及壓力明顯高于顆粒內(nèi)部；對(duì)于HMX為基的PBX炸藥，粘結(jié)劑處的溫度及壓力亦高于顆粒內(nèi)部；對(duì)于上述兩種沖擊響應(yīng)過(guò)程的比較表明粘結(jié)劑較軟，能起到一定的緩沖作用，并在變形過(guò)程中吸收沖擊能量，從而有效的降低了炸藥顆粒邊界處的溫度與壓力.

［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］劉凱欣，高凌天，鄭文剛.混凝土動(dòng)態(tài)破壞過(guò)程的數(shù)值模擬［J］.工程力學(xué)（增刊），2000：470－474.

［15］劉凱欣，高凌天.離散元法在求解三維沖擊動(dòng)力學(xué)問(wèn)題中的應(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］王文強(qiáng).離散元方法及其在材料和結(jié)構(gòu)力學(xué)響應(yīng)分析中的應(yīng)用［D］.合肥：中國(guó)科學(xué)技術(shù)大學(xué)，2000.

［19］傅華，趙峰，譚多望，王文強(qiáng)，尚海林.沖擊作用下HMX晶體孔洞塌縮熱點(diǎn)生成機(jī)制的細(xì)觀數(shù)值模擬［J］.高壓物理學(xué)報(bào)，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é)報(bào)，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

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

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

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