冲击作用下炸药热点形成的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.