白志玲，段卓平，黄风雷

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

0 引言

爆轰反应流模型是当前国内外爆轰物理学前沿研究的热点，构建具有广泛适应度的爆轰反应流模型，确定相互匹配的模型参数，实现炸药热-力学-化学响应的冲击起爆爆轰成长全过程的高保真计算，是炸药威力与冲击波感度匹配设计的基础。高聚物粘结炸药(PBX)作为一种典型的非均质固体炸药，其细观物理结构具有高度非均匀性，炸药组分颗粒的尺寸和形态各异，大小不同的孔洞缺陷随机分布，以致PBX冲击起爆爆轰成长过程非常复杂，受多种因素影响，除炸药的组成配方[1-3]、外界约束[4]、温度[5-7]和加载压力[8-10]外，与炸药自身的细观结构特征如炸药组分颗粒尺寸分布[11-13]、孔洞尺寸分布[14]等也密切相关，但现有的宏观和细观反应速率模型大多无法适应上述多种因素的变化，迫切需要发展高保真反映物理机制的冲击起爆模型。

目前普遍认为，非均质固体炸药的冲击起爆是由冲击波与炸药细观结构缺陷相互作用导致能量耗散引起局部高温造成的，即所谓的热点起爆机理。大量研究表明，孔洞塌缩是非均质固体炸药冲击起爆热点形成的主导机制[6，11-12]。近年来，段卓平团队建立并发展了PBX冲击起爆细观孔洞塌缩热点模型和DZK(Duan-Zhang-Kim)系列宏细观反应速率模型[12，15-19]，可较好地反映和预测加载压力、初始温度、炸药组分配比、炸药颗粒度、孔隙度及粘结剂含量和强度等多种因素对PBX冲击起爆爆轰成长过程的影响，具有较宽的适应性。DZK系列模型属于典型的3阶段式宏细观结合模型，解耦表征了冲击起爆孔洞塌缩热点点火阶段、低压慢反应阶段和高压快反应阶段的反应机理和贡献。需说明的是，DZK系列模型将炸药组分颗粒尺寸进行简单统计平均处理，采用平均颗粒度意义下的塌缩热点胞元代表炸药整体的热点点火过程，未考虑炸药颗粒尺寸分布、孔洞尺寸分布等，忽略了一些冲击起爆物理机制，因此模型还存在一定的局限性。此外，DZK系列模型的点火项描述细观尺度孔洞塌缩变形产生热点的过程，是通过子网格技术[20]的数值解实现的，推广到三维实际装置的工程计算时，工作量巨大。

目前，具有代表性的统计热点反应速率模型(如Grebenkin 3因素模型[21]、Cochran模型[22]和Nichols-Tarver模型[23]、Hill模型[24]和Hamate-Horie模型[25-26]等)可反映炸药初始孔洞尺寸分布对热点密度的影响，重点关注和描述了炸药冲击起爆过程中热点形成、点火、成长或消亡及快速转为爆轰的全过程。但是，上述统计热点反应速率模型均采用层流/表面燃烧机理描述热点点火后低压慢反应直至转为爆轰的过程，计算结果与实验结果差异较大。

本文构建非均质固体炸药冲击起爆统计热点反应速率模型时，采用点火-增长3阶段形式，热点项借鉴现有统计热点模型的思想，描述PBX冲击起爆过程中热点的形成、点火、成长或消亡；结合DZK反应速率模型后两项，描述点火后反应增长和快速转为爆轰的过程。将模型嵌入有限元程序DYNA-2D，数值模拟PBX冲击起爆和爆轰成长过程，与实验数据对比，验证模型的适应性。

1 PBX的孔洞尺寸分布

目前普遍认为，针对成熟定型炸药配方，炸药内部孔洞呈单/双/多峰分布特征主要与炸药颗粒度分布和级配有关。如图1(a)所示，Willey等[14]采用超小角X射线散射(USAXS)实验研究发现，三氨基三硝基苯(TATB)基LX-17炸药内孔洞体积随其尺寸满足双峰分布特征，并提出了孔洞体积-尺寸双峰/单峰分布函数表达式为

图1 LX-17和PBX9501炸药孔洞体积-尺寸呈双峰/单峰分布特征

(1)

式中：xv为孔洞尺寸；a1、a2、c1、c2、ω1和ω2均为常数，a1和a2代表双峰分布函数的两个峰值，c1和c2分别表示在对数坐标系下双峰分布函数的数学期望，ω1和ω2分别表示在对数坐标系下双峰分布函数的标准差。因此，在算术坐标系下，双峰分布函数的数学期望分别为10c1和10c2，标准差分别为10ω1和10ω2.LX-17炸药孔洞体积-尺寸双峰分布函数的参数值[14]如表1所示。

表1 LX-17和PBX9502炸药孔洞体积-尺寸双峰/单峰分布函数参数值

此外，Mang等[27]采用超小角中子散射技术测量了HMX基PBX9501炸药内孔洞尺寸呈单峰分布特征，如图1(b)所示，本文采用孔洞体积-尺寸单峰分布函数拟合实验数据，确定的参数值如表1所示。

2 统计热点DZK反应速率模型的构建

相比孔洞尺寸xv，热点尺寸xh非常小(几乎相差一个量级)，通常简单假设

xh=βxv，

(2)

式中：β≈0.1[21].联立(1)式和(2)式，获得潜在热点的体积-尺寸分布函数为

(3)

借鉴Cochran模型关于热点反应的统计处理方法[22]，假设冲击波作用下，在t时刻，宏观反应流场坐标系x位置处，单位体积炸药内产生M(x，r，t)Δr个半径为r的潜在球形热点，M(x，r，t)为局部球坐标系下潜在球形热点的数量密度，Δr为局部球坐标系下的无限小量长度；随后部分消亡，成为无效热点，而部分潜在热点发生反应，生成N(x，t，r)Δr个半径为r的球形反应核，N(x，r，t)为局部球坐标系下球形反应核的数量密度。于是，单位体积炸药的热点点火反应度为

(4)

式中：dc为临界热点尺寸(μm)，其大小取决于热点温度，同时与前导冲击波阵面压力pf[21]有关，

(5)

T0为炸药的初始温度；T1、η和γ为炸药的热力学相关常数；对于TATB，T1≈1 400 K，η=200 K/GPa，γ≈0.08[6]；对于HMX，T1≈1 000 K，η=150 K/GPa，γ≈0.06.

此外，潜在热点和反应核的统计动力学方程[22]为

(6)

(7)

为表征非均质固体炸药冲击起爆热点生成的主导机制，即孔隙塌缩，形核速率K采用以下表达式[23]：

(8)

式中：A为常数；p*为饱和压力(在压力峰值处，为避免出现不符合实际的较大的塌缩速率)；pc为热点燃烧率阈值压力，表示抗孔隙塌缩的能力；H表示阶跃函数。

(9)

采用矩量法求解热点反应控制方程(6)式和(7)式，获得单位体积炸药热点点火反应速率为

(10)

式中：λ表示在t时刻炸药的整体反应度(即已反应炸药的体积分数)。可见，热点点火项的待定热力学参数有T1、b、γ、v0、A、p*和pc.

热点点火反应结束后，炸药进入早期低压慢反应阶段，通常以表面燃烧反应机理描述；随着反应的进行，炸药随后进入高压快反应阶段，即转为爆轰。系列研究表明，DZK反应速率模型[12，16，20]即可较好地描述上述低压慢反应过程和高压快反应过程。于是，基于DZK模型[20]，构建统计热点DZK反应速率(SHS-DZK)模型：

(11)

式中：ro表示炸药颗粒度体积-尺寸分布函数的期望值，即统计平均粒度(μm)；a、b、G、m、n和s为常数，通过冲击起爆实验数据确定。(11)式中：第1项描述冲击载荷作用下热点成核反应过程(反应度为λi)；第2项描述热点形成后早期低压燃烧慢反应(反应度为λs；第3项为田占东等[28]提出的高压快反应速率方程(反应度为λf)，描述炸药整体爆轰反应。

统计热点DZK模型中点火项描述PBX炸药冲击起爆孔洞塌缩热点机制，但暂不关注具体的力学塌缩形式如粘塑性塌缩、射流等，热点项函数形式简单，嵌入现有限元程序或商用软件时，无需子网格技术，方便实际工程应用。

3 统计热点DZK反应速率模型的适应性

为对比验证统计热点DZK模型的适应性，分别将现有统计Cochran模型[22]和统计DZK模型(见(11)式)嵌入有限元程序DYNA2D，计算HMX基PBX9501(95% HMX，2.5% Estane，2.5% BDNPA-F nitroplasticizer)冲击起爆爆轰成长过程中不同Lagrange位置的压力成长历史，并与文献[25]中实验数据对比，如图2和图3所示。由图2和图3可看出：统计热点DZK模型计算结果与实验数据符合较好，波到达时间的计算值与实验值偏差不大于0.8%；而Cochran模型计算结果与实验结果差异较大，尤其是转爆轰后的前导冲击波阵面附近压力波形特征与实验现象明显不符。分析原因是：炸药在冲击起爆初期和接近爆轰阶段反应机制不同，冲击起爆初期表现为表面燃烧反应；而在接近爆轰阶段，高幅值冲击波作用下呈现整体均质快速反应，表明高压快反应阶段仍采用层流燃烧机理描述已不合适。

图2 Cochran模型计算的不同Lagrange位置压力历史与实验数据[25]对比(PBX9501)

图3 本文模型计算的不同Lagrange位置压力历史与实验数据[25]对比(PBX9501)

进一步计算TATB基LX-17(92.5%TATB，7.5%Kel-F)冲击起爆爆轰成长过程中不同Lagrange位置的压力成长历史，并与文献[2]中实验数据对比，如图4所示。由图4可看出，二者符合较好，波到达时间的计算值与实验值偏差不大于3.7%，验证了统计热点DZK模型的合理性和适应性。对于LX-17炸药的最后两个Lagrange位置，分析统计热点DZK模型计算的转爆轰时间比实验测试结果稍偏早的原因是，冲击起爆实验用LX-17炸药样品与小角散射实验测试孔洞尺寸分布的LX-17炸药样品是不同批次，虽然二者装药密度差异较小，但制备工艺、颗粒度级配等差异均会导致孔洞尺寸分布变化。

图4 LX-17炸药的不同Lagrange位置压力历史的数值模拟与实验数据[2]对比

对于HMX和TATB组分，统计热点DZK模型中热点点火项参数值[21]如表2所示，(11)式第2项、第3项参数[18]如表3所示。

表2 LX-17和PBX9501炸药的统计热点DZK模型中热点点火项参数[21]

表3 PBX9501和LX-17炸药的统计热点DZK模型中第2项、第3项参数[18]

此外，爆轰产物和未反应炸药状态方程均采用含温度形式的JWL状态方程[29]：

(12)

表4 常温下PBX9501和LX-17爆轰产物和未反应炸药JWL状态方程参数[1，29]

进一步提取PBX9501和LX-17冲击起爆爆轰成长过程中反应流场的不同Lagrange位置反应度-时间历史及统计热点DZK模型点火项、低压慢反应项和高压快反应项的反应度贡献，分别如图5和图6所示。图5和图6中，λf为高压快反应项，λs为低压慢反应项。由图5和图6中可以看出，由于流场中波的追赶及化学反应产生的压缩波，越靠后的Lagrange位置处反应速率越快，在形成爆轰的最后两个位置，反应在波阵面附近快速完成。此外，对于HMX基PBX9501炸药，低压慢反应阶段的燃烧反应贡献较多(λs占比较大)，表明HMX基PBX冲击起爆特性为点火后的加速反应特征；TATB基LX-17炸药的产物气体温度较低[21]，在低压慢反应阶段的燃烧反应贡献较少，表明TATB基PBX点火后呈稳定反应特征。

图5 PBX9501冲击起爆过程中反应度-时间历史

图6 LX-17冲击起爆过程中反应度-时间历史

如图7所示：载荷强度pf越高，PBX9501冲击起爆反应速率越快，波阵面附近的热点点火反应也明显越快，因为前导冲击波阵面压力越高，临界热点尺寸越小，能够发生反应的热点数量越多，表现为热点点火反应明显越快(见图7(b))；如果前导冲击波阵面压力足够强，波阵面附近的热点数量可达到饱和。此外，如图8(b)所示，初始温度越高，炸药内激发的反应热点数量也越多，热点点火反应也明显越快。由此可见，统计热点DZK模型可反映载荷强度和初始温度等变化对炸药热点点火的影响。值得解释的是，DZK系列模型中热点项定量描述PBX炸药冲击起爆孔洞塌缩热点机制，热点点火反应完成时，点火反应度约为1%[20]，因此，这里控制统计热点DZK模型的热点点火反应度达到1%时关闭热点项，启动后续燃烧反应项。

图7 在0 mm Lagrange位置处载荷强度对PBX9501炸药冲击起爆过程的影响(T0=25 ℃)

图8 在0 mm Lagrange位置处初始温度对PBX9501炸药冲击起爆过程的影响(pf=3.09 GPa)

为探索PBX炸药孔洞尺寸分布变化对其冲击起爆特性的影响，基于装药密度ρc=1.84 g/cm3的PBX9501炸药内孔洞体积-尺寸分布数据(见表1或表5的类型2)，通过调整分布函数(见(1)式)的期望值或标准差(参数值见表5)，可改变平均孔径或孔径分散度等特征，如图9(a)所示。积分后的孔洞体积分数即孔隙度保持一致，如图9(b)所示。相比类型2的孔洞分布，类型1集中较小尺寸孔洞，类型5集中较大尺寸孔洞，而类型3和类型4的孔洞尺寸更分散。

表5 具有相同孔隙度的PBX9501炸药多种孔洞体积-尺寸分布特征参数

图9 孔隙度相同的PBX9501炸药多种孔洞尺寸分布特征

固定PBX9501炸药的孔隙度和颗粒度等参数，在同一载荷条件下，炸药内孔洞尺寸分布变化对其冲击起爆爆轰成长过程中反应度-时间历史的影响如图10所示。上述5种分布特征中：集中分布较多小孔洞的炸药(类型1)或集中分布较多大孔洞的炸药(类型5)反应速率最慢，因为小孔洞能够形成的反应热点数量最少，热点点火反应速率最慢，而大孔洞容易形成反应热点，但同一孔隙度下，大孔洞的数量少，形成的热点数量少，热点点火反应速率也慢；孔洞尺寸分布越分散(类型4>类型3>类型2)，反应速率越慢，因为尺寸分散度越大，较小孔洞的占比越多，形成无效热点的比例越大，导致热点点火反应速率下降。由此可见，实际装药制备时，尽量控制使炸药内孔洞尺寸较小，炸药不易起爆，炸药的安全性更好，与基本认识一致。上述结果初步表明，目前统计热点DZK模型可反映炸药内部孔洞尺寸分布变化对其冲击起爆爆轰成长过程的影响。

图10 孔洞尺寸分布对PBX9501炸药冲击起爆过程的影响

4 结论

本文构建的统计热点反应速率模型，可反映PBX炸药内部孔洞尺寸分布、颗粒度、载荷强度、初始温度及炸药组分热力学参数等对炸药冲击起爆爆轰成长过程的影响。相比文献[21-24]的统计模型，该模型适应性更强，与实验结果吻合更好。得到以下主要结论：

1)HMX基PBX的临界起爆压力低，随着反应增长，波阵面压力升高，热点数量增多，冲击起爆过程受波阵面热点点火和燃烧反应共同作用，其燃烧反应速度较快，表现为加速反应特性。TATB基钝感PBX的临界起爆压力高，波阵面热点数量几乎饱和，冲击起爆主要受点火后的燃烧反应过程控制；由于TATB产物气体温度较低，点火后燃烧反应速度较慢，表现为稳定反应特性。进一步提高了对TATB/HMX混合基钝感高能炸药冲击起爆机理的认识。

2)载荷强度越高或初始温度越高，临界热点尺寸越小，满足点火条件的热点数量越多，炸药点火反应越快；对于相同孔隙率的炸药，小尺寸孔洞越多或孔洞尺寸分布越分散，满足点火条件的热点数量越少，炸药点火反应速率越小。实际装药制备时，尽量控制使得炸药孔洞尺寸较小，则炸药不易起爆，炸药安全性提高。

参考文献(References)

[1] GUSTAVSEN R L，SHEFFIELD S A，ALCON R R，et al.Embedded electromagnetic gauge measurements and modeling of shock initiation in the TATB based explosives LX-17 and PBX 9502[J].AIP Conference Proceedings，2002，620:1019-1022.

[2] URTIEW P A，TARVER C M.Shock initiation of energetic materials at different initial temperatures(review)[J].Combustion，Explosion，and Shock Waves，2005，41(6):766-776.

[3] 温丽晶，段卓平，张震宇，等.HMX基和TATB基PBX炸药爆轰成长差别的实验研究[J].爆炸与冲击，2013，33(增刊1): 135-139.

WEN L J，DUAN Z P，ZHANG Z Y，et al.Experimental research on differences of detonation growth process between HMX-based and TATB-based plastic bonded explosives[J].Explosion and Shock Waves，2013，33(S1):135-139.(in Chinese)

[4] FORBES J W，TARVER C M，URTIEW P A，et al.The effects of confinement and temperature on the shock sensitivity of solid explosives[C]∥Proceeding of the 11th International Detonation Symposium.Snowmass，CO US: USDOE Office of Defense Programs，1998.

[5] MULFORD R N，ALCON R R.Shock initiation of PBX-9502 at elevated temperatures[J].AIP Conference Proceedings，1996，370:855-858.

[6] PERRY W L，CLEMENTS B，MA X，et al.Relating microstructure，temperature，and chemistry to explosive ignition and shock sensitivity[J].Combustion and Flame，2018，190:171-176.

[7] URTIEW P A，TARVER C M，FORBES J W，et al.Shock sensitivity of LX-04 at elevated temperatures[C]∥Proceedings of 1997 APS Topical Conference on Shock Compression of Condensed Matter.Amherst，MA，US:American Physical Society，1997.

[8] SIMPSON R L，URTIEW P A，TARVER C M.Shock initiation of 1，3，3-trinitroazetidine(TNAZ)[J].AIP Conference Proceedings，1996，370:883.

[9] TARVER C M，SIMPSON R L，URTIEW P A.Shock initiation of an ε-CL-20-estane formulation[J].AIP Conference Proceedings，1996，370:891.

[10] URTIEW P A，TARVER C M，SIMPSON R L.Shock initiation of 2，4-dinitroimidazole(2，4-DNI)[J].AIP Conference Proceedings，1996，370:887.

[11] KHASAINOV B A，ERMOLAEV B S，PRESLES H N，et al.On the effect of grain size on shock sensitivity of heterogeneous high explosives[J].Shock Waves，1997，7(2): 89-105.

[12] WEN L J，DUAN Z P，ZHANG L S，et al.Effects of HMX particle size on the shock initiation of PBXC03 explosive[J].International Journal of Nonlinear Sciences and Numerical Simulation，2012，13(2):189-194.

[13] MOULARD H，DELCLOS A，KURY J.The effect of RDX particle size on the shock sensitivity of cast PBX formulations:2，Bimodal compositions[C]∥Proceedings of International Symposium on Pyrotechnics and Explosives.Beijing，China:China Ordnance Society，1987.

[14] WILLEY T M，BUUREN T V，LEE J R I，et al.Changes in pore size distribution upon thermal cycling of TATB-based explosives measured by ultra-small angle X-ray scattering[J].Propellants，Explosives，Pyrotechnics，2010，31(6):466-471.

[15] LIU Y R，DUAN Z P，ZHANG Z Y，et al.A mesoscopic reaction rate model for shock initiation of multi-component PBX explosives[J].Journal of Hazardous Materials，2016，317: 44-51.

[16] DUAN Z P，WEN L J，LIU Y R，et al.A pore collapse model for hot-spot ignition in shocked multi-component explosives[J].International Journal of Nonlinear Sciences and Numerical Simulation，2010，11:19-23.

[17] BAI Z L，DUAN Z P，WEN L J，et al.Comparative analysis of detonation growth characteristics between HMX-and TATB-based PBXs[J].Propellants，Explosives，Pyrotechnic，2019，44(7):858-869.

[18] BAI Z L，DUAN Z P，WEN L J，et al.Shock initiation of multi-component insensitive PBX explosives: experiments and MC-DZK mesoscopic reaction rate model[J].Journal of Hazardous Materials，2019，369:62-69.

[19] BAI Z L，DUAN Z P，WEN L J，et al.Embedded manganin gauge measurements and modeling of shock initiation in HMX-based PBX explosives with different particle sizes and porosities[J].Propellants，Explosives，Pyrotechnics，2020，45(6):908-920.

[20] 白志玲.PBX炸药冲击起爆机理及其系列反应速率模型研究[D].北京：北京理工大学，2019.

BAI Z L.Physical mechanism and series of chemical reaction rate models for detonation initiation in PBX explosives[D].Beijing:Beijing Institute of Technology，2019.(in Chinese)

[21] GREBENKIN K F.Comparative analysis of physical mechanisms of detonation initiation in HMX and in a low-sensitive explosive(TATB)[J].Combustion Explosion and Shock Waves，2009,45(1):78-87.

[22] COCHRAN S G.Statistical treatment of heterogeneous chemical reaction in shock-initiated explosives:UCID-18548[R].Livermore，CA，US:Lawrence Livermore National Laboratory，1980.

[23] NICHOLS A L，TARVER C M.A statistical hot spot reactive flow model for shock initiation and detonation of solid high explosives[C]∥Proceedings of the 12th International Detonation Sympo-sium.San Diego，CA，US:US Department of Energy，2002.

[24] HILL L.The Shock-triggered statistical hot spot model[J].AIP Conference Proceedings，2012，1426:307-310.

[25] HAMATE Y，HORIE Y.Ignition and detonation of solid explosives:a micromechanical burn model[J].Shock Waves，2006，16(2):125-147.

[26] HAMATE Y，HORIE Y.A statistical approach on mechanistic modeling of high-explosive ignition[J].AIP Conference Proceedings，2004，706:335-338.

[27] MANG J T，HJELM R P，FRANCOIS E G.Measurement of porosity in a composite high explosive as a function of pressing conditions by ultra-small-angle neutron scattering with contrast variation[J].Propellants，Explosives，Pyrotechnics，2010，35(1):7-14.

[28] 田占东，张震宇.PBX-9404炸药冲击起爆细观反应速率模型[J].含能材料，2007，15(5):464-467.

TIAN Z D，ZHANG Z Y.A mesomechanic model of shock initiation in PBX-9404 explosive[J].Chinese Journal of Energetic Materials，2007，15(5):464-467.(in Chinese)

[29] URTIEW P A，VANDERSALL K S，TARVER C M，et al.Initiation of heated PBX-9501 explosive when exposed to dynamic loading[C]∥Proceedings of AIP ZABABAKHIN SCIENTIFIC TALKS-2005:International Conference on High Energy Density Physics.Snezhinsk，Russia:AIP，2005.