炸药颗粒床冲击加热多尺度数值计算

张新明，吴艳青，黄风雷

(北京理工大学 爆炸科学与技术国家重点实验室，北京100081)

0 引言

塑料黏结炸药(PBXs)可看作是一种特殊的颗粒复合材料，其主要组成部分包括HMX 含能颗粒(质量分数约占95%)、黏结剂、增塑剂及钝感剂等(质量分数约占5%).PBXs 在受到冲击加载时，HMX 炸药颗粒为主要的承载对象，因此，在进行PBXs 冲击起爆研究时，为了简化模型，有时忽略黏结剂、增塑剂等的作用，将PBXs 看成由HMX 炸药颗粒构成的颗粒床进行研究。另外，对于无空隙的均质含能材料，由于化学构成的不稳定性、老化、损伤等因素会出现孔隙，研究这类材料的性能，也经常将它们看成颗粒炸药床进行研究［1］。

燃烧转爆轰(DDT)的研究已有几十年的历史［2-10］，也是颗粒含能材料受冲击加载时常常产生的一个现象。文献［3 -8］均表明颗粒含能材料的燃烧转爆轰是颗粒材料的压实(Material Compaction)和热量释放(Heat Release)的相互作用，而多数情况下，燃烧转爆轰都与致密波(Compaction Wave)作用密切相关［11］。所谓致密波，即为使颗粒床孔隙率发生变化的行进波［11］，Sandusky 等［12］最早在试验中观察到稳态致密波。在致密区内，由于颗粒与颗粒之间的摩擦、应力集中、塑性变形等原因会导致材料的局部温度大幅升高，形成所谓的热点(Hot-spot)，在整体平均温度很低，甚至还没达到熔点的情况下，就有可能使材料局部熔化，发生相变，激起剧烈的化学反应，并最终导致爆燃甚至是爆轰。自上世纪80年代始，国内外很多学者开展了这方面的实验、理论和数值研究［13-19］。

但宏观尺度的整体平均模型［13-17］不能有效地反映出颗粒尺度的响应，而基于单个颗粒的统计模型［18-19］由于计算量太大，不适合工程计算。近年来，Gonthier 等［20］根据Hertz 接触力学建立了颗粒材料冲击加热的多尺度模型，该模型能够有效地反映出颗粒尺度的响应，且适合大量工程计算。但是，Gonthier［20］没有考虑颗粒可压缩性和固液相变对压缩过程的影响，并且颗粒构形应力函数由准静态压缩试验测得，工作量大。本文选取Tait 状态方程，由颗粒构形应力的定义式来确定该函数，避免进行大量试验，从而减小工作量;同时引入化学反应进程变量，并考虑炸药颗粒的可压缩性以及固液相变的影响，对活塞冲击颗粒床进行多尺度数值计算。

1 计算模型

Baer 等［7］的两相流模型，通常称为BN 模型，成为国内外学者研究致密波及含能材料燃烧转爆轰基础，本文在BN 模型的基础上进行简化，并结合颗粒尺度热点模型和阿伦尼乌斯反应速率方程，对活塞冲击颗粒床进行数值计算。

1.1 一维单相流模型

根据Baer 等［13］和Lovey 等［21］的工作，气相对致密波的影响很小，本文忽略气相作用，一维单相流模型如下:

式中:t 为时间;x 为位置;ρs为固相密度;φs为固相体积分数;vs为颗粒速度;ps为压强;E 为单位体积炸药的总能量，其中E =es+B，es为单位体积炸药的内能，B 为单位体积炸药的压缩潜能;β 为颗粒构形应力函数，表示颗粒床受到挤压后，颗粒之间的相互作用力;常数μc为压缩粘性。(1)～(3)式分别为质量、动量和能量守恒方程;(4)式表示由于固相压力与颗粒构形应力存在压差而引起的固相体积分数的变化。

1.2 颗粒尺度的热点模型

图1为Gonthier［20］根据Hertz 接触力学建立的热点模型。图中，rc为接触球半径，r0为接触球半径的最大值。当颗粒床受到外载荷时，颗粒与颗粒之间由于摩擦、应力集中、塑性变形等因素使得颗粒接触面附近温度升高，形成所谓的“热点”，在整体平均温度很低的情况下有可能发生燃烧转爆轰。

图1 颗粒与颗粒间接触的几何模型Fig.1 Illustration of the intergranular contact geometry used for the thermal energy localization strategy

为简化模型，假设HMX 颗粒不在发生破碎和聚并，颗粒半径为R，单位体积内的颗粒数n，颗粒数守恒方程为

本文忽略颗粒与颗粒之间的摩擦作用，单位体积内塑性功功率等于内能的变化率，

式中:nc为单位体积内的接触球个数;Vc=4/3r3c为单个接触球的体积;σY=3.0Y 为塑形应力，Y 为HMX的屈服强度;des/dt 为内能变化率。

阿伦尼乌斯反应速率方程［20］

式中:λ(x，r，t)(0≤λ≤1)为沿接触球径向(r 方向，见图1)变化的化学反应过程变量，λ =0 表示无化学反应发生，λ=1 表示完全反应;z 为阿伦尼乌斯因子;Tα为活化温度。

内能变化率

式中:q 为传导热通量，q = -k ∂T(x，r，t)/∂r，k 为热传导率;ρs0为固相初始密度。(8)式右端第1 项为由热传导引起的内能变化率，第2 项为由化学反应引起的内能变化率，第3 项为由颗粒材料孔隙率变化引起的内能变化率Sφ以及由颗粒密度变化引起的内能变化率Sρ.

HMX 在常温常压下的熔点T0m为552 K，根据Kraut-Kennedy［22］关系，熔点Tm= T0m(1 +αΔV/V0)，其中，常数α =1.53，V 为比容，V0为初始比容;在ps≈5 GPa 的条件下，ΔV≈0.1，熔点升高到约636 K，而本文条件下ps≪5 GPa，因此假设Tm=T0m，颗粒床发生等温相变:

式中:cvs为定容比热;q0m为熔化潜热;0≤χ≤1 为液相质量分数。

1.3 本构关系

为使方程组封闭，还需要给出固相状态方程和颗粒构形应力函数β.本文对Powers 等［23］给出的状态方程进行修改，引入化学反应进程变量λs，Helmholtz 自由能的表达式为:

式中:Ts0为固相初始温度;ps0为固相初始压强;pg0为气相初始压强;φs0为固相初始体积分数。

状态方程

颗粒构形应力函数β 的表达式如下:

Tait 状态方程中的参数如表1所示。

表1 Tait 状态方程中的参数Tab.1 Parameters in Tait EOS for solid HMX

2 计算模型的无量纲化

本文以HMX 颗粒为例，对其致密波传播特性进行分析，表2为HMX 的物性参数及初始条件。

表2 HMX 物理性质参数及初始条件Tab.2 HMX material properties and initial conditions

假设致密波以波速D 向右传播，利用伽利略坐标变换公式:相对位置ξ=x-D·t 和相对速度v=vs-D 以及定常假设，∂/∂t ≡0.定义如下无量纲量:

将计算模型写成无量纲形式

3 计算结果

为确定压缩过程的终态，令β =ps，即颗粒床达到压力平衡。图2为颗粒构形应力函数β 与固相体积分数φs之间的关系。由图可知，β 为体积分数φs的单调递增函数。图3为不同致密波速D，固相压力Ps与体积分数φs之间的关系。由图可知，随着致密波速度的增加，终态体积分数φs逐渐增加，当D ＞500 m/s 时，终态体积分数φs≈1，即颗粒床处于完全压实状态。

图2 颗粒构形应力β 与固相体积分数φs之间的关系Fig.2 Intergranular stress β vs.solid volume fraction φs

图3 固相压力ps与固相体积分数φs之间的关系Fig.3 Solid pressure ps vs.solid volume fraction φs

3.1 惰性颗粒药床

图4 初始体积分数φs0 =0.76，活塞速度up =100 m/s 时模型物理量在ξ 方向上的分布Fig.4 Variation of model parameters in ξ direction for φs0 =0.76 and up =100 m/s

图4为活塞速度up=100 m/s 撞击颗粒药床，产生亚音速致密波(D=418.7 m/s ＜C0)，固相体积分数φs、压力ps、颗粒间的应力β、整体平均温度Ts和颗粒密度ρs的分布。由图4可知，各物理量由初始状态光滑地过渡到平衡状态，终态值φs=0.995，ps=67.73 MPa，Ts=304.08 K，ρs=1 906.62 kg/m3.由图4(b)可知，在初始阶段颗粒床处于平衡状态，即Ps0=β(φs0)，随着活塞向右运动，压力增加，颗粒间构形应力增加，当ps=β(φs)时，颗粒床再次达到压力平衡状态。

图5为up=100 m/s 撞击颗粒床时，颗粒半径R、接触球半径rc和颗粒温度T 的分布。由图5(a)可知，压缩过程中颗粒半径R 几乎保持不变，接触球半径rc由初始值rc0≈1.859 8 μm 光滑地过渡到稳态值rc≈4.849 4 μm;由图5(b)可知，由于颗粒之间的挤压，发生塑性变形，在颗粒接触面附近，颗粒温度明显升高，最大值T≈973 K，超过了HMX 的熔点T0m=552 K，并可以观察到相变发生。而由图4(c)可知，相同的加载条件下，整体平均温度Ts≈304 K，远低于颗粒尺度温度的最大值，即整体平均模型不能有效地反应颗粒尺度的响应。

图5 颗粒半径和颗粒尺度温度的分布(up =100 m/s)Fig.5 Variation of grain radius in ξ direction and grainscale temperature in ξ and radial directions at up =100 m/s

3.2 考虑化学反应

图6为考虑固相可压缩性和相变的情况下，化学反应进程变量λ、颗粒温度T、整体平均反应进程变量λs和整体平均温度Ts的分布。由(19)式可知，由于考虑了固相可压缩性，接触球半径rc会增大，这就意味着内能将会分布在颗粒接触面附近更大的区域上，从而使颗粒温度降低;另外由于相变会消耗掉一部分能量，也将使颗粒温度降低，因此在相同的初始条件下要激发化学反应必须增大活塞速度。当活塞速度up=79.0 m/s 时，致密波速度D=348.75 m/s，由图6(a)和图6(b)可知，在颗粒接触面附近化学反应进程变量λ =1，即完全反应，相应地温度出现了突跳，颗粒温度的最大值为1 638 K，并且可以观察到相变的发生。由图6(c)和图6(d)可知，整体平均化学反应进程变量λ≈7.5 ×10-7，整体平均温度Ts≈302.5 K.结果再次表明整体平均模型不能有效地反应颗粒尺度的响应。

3.3 参数研究

图7为活塞速度up、初始体积分数φs0、颗粒半径R 等模型参数对颗粒床压缩过程的影响。由图7(a)可知，随着活塞速度up的增加，固相体积分数φs逐渐增加，up～φs曲线的斜率逐渐减小，表明颗粒床更难被压实;当活塞速度up＞100 m/s 时，固相体积分数φs≈1，即颗粒药床被完全压实。图7(b)为活塞速度up与致密波速度D 之间的关系，由图可知，随着活塞速度up的增加，致密波速度D 逐渐增加，且两者之间近似成线性关系。由图7(c)可知，随着初始体积分数φs0的增加，要想激起化学反应，所需压力逐渐增加;当φs0＞0.90 时，所需压力急剧增加，也即致密的颗粒床更难发生化学反应。由图7(d)可知，随着颗粒半径R 的增加，激起化学反应所需的活塞速度逐渐减小，即粗颗粒床更容易发生化学反应，这一结果与已有的试验和数值模拟结果［25-28］相一致。

4 结论

本文对活塞冲击HMX 颗粒床进行多尺度数值模拟。结果表明，对于亚音速致密波，固相体积分数、压力、密度和整体平均温度光滑地从初始值过渡到平衡状态，致密波速度随着活塞速度的增加而逐渐增加，颗粒床的致密程度随着致密波速度的增加而增加，直到药床被完全压实(φs≈1).颗粒尺度的温度分布表明，在相同的加载条件下(up=100 m/s)，颗粒尺度的温度明显高于整体平均温度，并超过了熔点，同时可以观察到固液相变的发生，结果表明整体平均模型不能有效地反映出材料颗粒尺度的响应。在考虑化学反应的情况下，由于颗粒材料可压缩性和固液相变的影响，在相同的初始条件下，必须增大活塞速度才能激起化学反应。计算结果还表明，随着固相初始体积分数的增大，发生化学反应所需的外界载荷越大，即致密床更难发生化学反应;随着颗粒半径的增大，发生化学反应所需的活塞速度越小，即粗颗粒床更容易发生化学反应。

图6 化学反应进程变量和温度的分布Fig.6 Predicted chemical reaction progress variable and temperature through the compaction zone

图7 模型参数研究Fig.7 Parametric sensitivity analysis results

References)

［1］ Greenaway M W.Measurement of intergranular stress and porosity during dynamic compaction of porous beds of HMX［J］.Journal of Applied Physics，2005，97(9):093521 -093521-7.

［2］ Gonthier K A，Powers J M.A high resolution numerical method for a two-phase model of deflagration-to-detonation transition［J］.Journal of Computational Physics，2000，163(2):376 -433.

［3］ Cochran M T，Powers J M.Computation of compaction in compressible granular material［J］.Mechanics Research Communications，2008，35(1 -2):96 -103.

［4］ Bernecker R R，Price D.Studies in the transition from deflagration to detonation in granular explosives-Ⅰ.Experimental arrangement and behavior of explosives which fail to exhibit detonation［J］.Combust and Flame，1974，22(1):111 -118.

［5］ Bernecker R R，Price D.Studies in the transition from deflagration to detonation in granular explosives-Ⅱ.Transitional characteristics and mechanics observed in 91/9 RDX/wax［J］.Combust and Flame，1974，22(1):119 -129.

［6］ Butler P B，Krier H.Analysis of deflagration to detonation transition in high-energy solid propellant［J］.Combust.Flame，1986，63(1 -2):31 -48.

［7］ Bear M R，Nunziato J W.A two-phase mixture theory for the deflagration to detonation transition (DDT)in reactive granular materials［J］.International Journal of Multiphase Flow，1986，12(6):861 -889.

［8］ Xu S J，Stewart D S.Deflagration-to-detonation transition in porous energetic materials:a comparative model study［J］.Journal of Engineering Mechanics，1997，31(2 -3):143 -172.

［9］ Kapila A K，Menikoff R，Bdzil J B，et al.Two-phase modeling of deflagration-to-detonation transition in granular materials:reduced equations［J］.Physics of Fluids，2001，13(10):3002 -3024.

［10］ Butler P B，Krier H.Modeling deflagration to detonation transition in granular explosives pentaerythritol tetranitrate［J］.Journal of Applied Physics，2008，104(4):043519-(043519-14).

［11］ Powers J M，Stewart D S，Krier H.Analysis of steady compaction waves in porous materials［J］.Journal of Applied Mechanics，1989，56(1):15 -24.

［12］ Sandusky H W，Liddiard T P.Dynamic compaction of porous beds［R］.US:Naval Surface Warfare Center Silver Spring MD，1985.

［13］ Baer M R.Numerical studies of dynamic compaction of inert and energetic granular materials［J］.Journal of Applied Mechanics，1988，55:36 -43.

［14］ Gonthier K A，Menikoff R，Son S E，et al.Molding compactioninduced energy dissipation of granular HMX［C］∥11th Symposium (International)on Detonation.1998:153 -161.

［15］ 孙锦山，朱建士，贾祥瑞.颗粒材料中致密波结构研究［J］.力学学报，1999，31(4):423 -433.SUN Jin-shan，ZHU Jian-shi，JIA Xiang-rui.An analysis of compaction wave in granular material［J］.Journal of Theoretical and Applied Mechanics，1999，31(4):423 -433.(in Chinese)

［16］ Low C A，Greenaway M W.Compaction process in granular beds composed of different particle sizes［J］.Journal of Applied Physics，2005，98(12):123519-(123519 -12).

［17］ Low C A，Longbottom A W.Effect of particle distribution on the compaction behavior of granular beds［J］.Physics of Fluids，2006，18(6):066101-(066101 -18).

［18］ Menikoff R，Kober E.Compaction waves in granular HMX［J］.Shock Compression of Condensed Matter，1999，505:397 -400.

［19］ Bardenhagen S G，Brackbill J U.Dynamic stress bridging in granular material［J］.Journal of Applied Physics，1998，83(11):5732 -5740.

［20］ Gonthier K A.Modeling and analysis of reactive compaction for granular energetic solids［J］.Combustion Science and Technology，2003，175:1679 -1709.

［21］ Lovey A，Ben-Dor G，Sorek S，et al.Jump conditions across strong compaction waves in gas saturated rigid porous media［J］.Journal of Shock Waves，1993，(3):105 -111.

［22］ Menikoff R，Sewell T D.Constituent properties of HMX needed for mesoscale simulations［J］.Combustion Theory and Modeling，2002，(6):103 -125.

［23］ Powers J M.Two-phase viscous modeling of compaction of granular materials［J］.Physics of Fluids，2004，16(8):2975 -2990.

［24］ Elban W L，Chiarito M A.Quasi-static compaction study of coarse HMX explosive［J］.Power Technology，1986，46(2 -3):181 -193.

［25］ Moulard H，Kury J W，Declos A.The effect of RDX particle size on the shock sensitivity of cast PBX formulations［C］∥8th Symposium (International)on Detonation.1985:902 -913.

［26］ Moulard H.Particular Aspect of the explosive particle size effect on shock sensitivity of cast PBX formulations［C］∥9th Symposium (International)on Detonation.1989:18 -24.

［27］ Simpson R L，Helm F H，Crawford P C，et al.Particle size effects in the initiation of explosives containing reactive and nonreactive continuous phases［C］∥9th Symposium (International)on Detonation.1989:25 -38.

［28］ Gustavsen R L，Sheffield S A，Alcon R R.Low pressure shock initiation of porous HMX for two grain size distribution and two densities［J］.American Institute of Physics，1996，370:851 -854.