高聚物黏结炸药的分子模拟进展

龙瑶; 陈军

doi:10.11858/gywlxb.20190755

高聚物黏结炸药的分子模拟进展

doi: 10.11858/gywlxb.20190755

龙瑶,
陈军^*, ,

北京应用物理与计算数学研究所，北京　100088

基金项目: 科学挑战专题（TZ2016001）

详细信息

作者简介:
龙　瑶（1980－），男，博士，研究员，主要从事含能材料研究. E-mail：long_yao@iapcm.ac.cn

通讯作者:
陈　军（1969－），男，博士，研究员，主要从事计算物理研究. E-mail：jun_chen@iapcm.ac.cn

中图分类号: O521.2
计量
- 文章访问数: 8454
- HTML全文浏览量: 2684
- PDF下载量: 68
出版历程
- 收稿日期: 2019-04-03
- 修回日期: 2019-04-25

Progress of Atomistic Simulations for Plastic Bonded Explosives

LONG Yao,
CHEN Jun^{*
, ,}

Beijing Institute of Applied Physics and Computation Mathematics, Beijing 100088, China

摘要

摘要: 回顾了近年来在高聚物黏结炸药(PBX)原子和分子尺度数值模拟方面取得的进展，主要研究领域包括以下6个方面：炸药分子力场、热力学参数计算、耗散/输运性能、相图/相变动力学、动力学响应行为和热点形成机制。针对当前研究现状，介绍了各领域的代表性工作和主要研究成果。目前对PBX炸药的结构和静力学性能已有较充分的认识，但对炸药的动力学响应行为和细观起爆机制尚缺少系统的科学认识，存在一系列挑战性问题，如结构缺陷在爆轰反应后期的形态和表征，以及初始缺陷对爆轰波波形畸变的影响机制。需要将理论计算与实验相结合，以解决爆轰物理领域中的难点问题。
- PBX炸药 /
- 分子动力学 /
- 热力学参数
Abstract: The recent atomistic simulations of plastic bonded explosive is reviewed in six aspects: the force-field, thermodynamic property, dissipation/transport property, phase transition, constitutive relation and ignition mechanism. In past decades, the structure and mechanical property of PBX are carefully investigated. However, the microscopic defect evolution and hot spot formation mechanisms are unclear. There are a set of challenging problems in detonation physics, such as the defect configuration at the chemical reaction zone, and the detonation wave deformation induced by defect. To investigate them, both atomistic simulation and experiment are required.
- plastic bonded explosive /
- molecular dynamics /
- thermodynamic properties

HTML全文

含能材料在现代国防工业、武器装备和民用工程中有广泛的应用，是先进武器战斗部、航天器推进剂和工程爆破装置的主要做功元件。目前，高聚物黏结炸药(PBX)已逐步取代传统炸药成为武器系统主装药。与梯恩梯(TNT)等传统单质炸药相比，PBX炸药由基底炸药、降感剂、黏结剂、耦联剂、增塑剂、热安定剂组成，各组分的形态和性能存在较大差异。性能迥异的组分结合在一起使得PBX炸药在外界刺激下表现出复杂的物理-力学耦合行为。

对PBX炸药物性、安全性和起爆特性的数值模拟研究涉及多个方面，包括热力学参数、缺陷演化、热点起爆等。PBX炸药由于具有复杂的组分和构型，其内部存在多种缺陷和微结构，如孔洞、位错、界面、裂纹等。缺陷在外界刺激作用下的演化过程使得PBX炸药表现出复杂的力学响应行为和热点起爆机制，因此缺陷演化是预测PBX炸药安全性和爆轰性能的重要依据。目前对PBX炸药的主要物性和起爆机制已开展了大量实验测量和数值模拟工作。如：通过X射线衍射、荧光光谱分析和中子散射技术对炸药的振动模式、缺陷形态和反应机理开展了系统的实验研究，积累了大量数据；通过第一性原理计算、分子动力学模拟、蒙特卡罗方法、物质点法、有限元法对炸药从微观到宏观的物理参数和动力学行为开展了大规模数值模拟，取得大量新成果。本文主要回顾近年来在微观数值模拟领域取得的成果，分6个方面阐述相关研究进展：（1）对炸药分子力场的研究工作，包括单质炸药力场、界面力场和化学反应力场；（2）对单晶、界面和复合炸药体系的热力学参数计算；（3）对PBX炸药耗散、输运性能的研究，包括热传导率和热耗散率；（4）对炸药相图和相变动力学的研究；（5）对动力学响应行为和本构关系的研究；（6）对炸药起爆初期的热点形成机制和降感机理研究。下面依次介绍各方面的研究工作。

1. PBX炸药分子力场研究

在微观数值模拟研究中，原子间相互作用势是材料模拟的关键。目前对金属材料、半导体材料、陶瓷材料已构建了较为成熟的嵌入原子势、Tersoff势等势模型。无机材料的势模型较容易构建，因为只需考虑原子间的相互作用。含能材料一般是分子晶体结构，需要考虑分子内部的相互作用和分子间的长程作用，力场构建较复杂。现在对于分子间相互作用机理的科学认识存在不足，只是简单地将该作用分为氢键、范德华力等。早期的炸药分子力场由于计算条件的限制，假定分子为刚体，只需要描述分子间的范德华力，由Sorescu等^[1]完成。该力场在预测奥克托金(HMX)炸药的晶格常数和等温线方面取得了一定成果^[2]；但由于其将分子设定为刚体，因此该力场无法正确描述在高压条件下炸药的分子形变过程，对高压区域的状态方程和物性预测不准确。

包含分子内和分子间相互作用的完整炸药力场是由Smith等^[3–5]开发的。该力场的总能可写为

${E_{{\rm{total}}}} = {E_{{\rm{pair}}}} + {E_{{\rm{Coul}}}} + {E_{{\rm{bond}}}} + {E_{{\rm{angle}}}} + {E_{{\rm{dihedral}}}} + {E_{{\rm{improper}}}}$

(1)

式中：E_pair表示范德华相互作用能，E_Coul表示库仑作用，E_bond表示键能，E_angle表示角度项，E_dihedral表示二面角项，E_improper表示奇异项。该力场能够正确预测晶体结构^[4–5]、状态方程^[6]、弹性常数^[6]和热导率^[7–9]，得到科研人员的广泛使用。Smith力场采用价键相互作用描述分子内部的力场，要求在分子动力学模拟过程中炸药分子不能发生断键和热分解行为，因此对描述炸药高温区域的性能有些不足。其他关于炸药体系的价键力场包括Gee开发的三氨基三硝基苯(TATB)分子力场^[10]，以及宋华杰开发的黑索金(RDX)分子力场^[11]，均得到广泛的应用。

在炸药添加剂的研究方面，孙淮等^[12–13]开发了适合聚合物添加剂的COMPASS力场，能够描述PBX炸药里常用的黏结剂，如石蜡、F2311、F2312、F2313、F2314、维通橡胶等。COMPASS力场的能量表达式为

$\begin{split} {E_{{\rm{total}}}} = & \sum\limits_b {[{k_2}{{(b - {b_0})}^2} + {k_3}{{(b - {b_0})}^3} + {k_4}{{(b - {b_0})}^4}]} + \sum\limits_\theta {[{k_2}{{(\theta - {\theta _0})}^2} + {k_3}{{(\theta - {\theta _0})}^3} + {k_4}{{(\theta - {\theta _0})}^4}]} + \\ & \sum\limits_{b,b'} {{k_{b,b'}}(b - {b_0})(b' - {b_0'})} + \sum\limits_{b,\theta } {{k_{b,\theta }}(b - {b_0})(\theta - {\theta _0})} + \sum\limits_{i,j} {\frac{{{q_i}{q_j}}}{{{r_{ij}}}}} + \sum\limits_{i,j} {{\varepsilon _{ij}}\left[ {2{{\left( {\frac{{r_{ij}^0}}{{{r_{ij}}}}} \right)}^9} - 3{{\left( {\frac{{r_{ij}^0}}{{{r_{ij}}}}} \right)}^6}} \right]} \end{split}$

(2)

式中：等号右边第1项表示键能，第2项表示键角能，第3项是化学键耦合能，第4项是键长-键角耦合能，第5项是库仑相互作用，第6项是范德华相互作用。COMPASS力场采用一种特殊的价键形式，考虑了键-键耦合等更为复杂的相互作用类型，能够比较准确地描述聚合物的链状分子结构。其他添加剂(如石墨)可以用传统的Tersoff力场^[14]描述。

陈军等^[15–19]基于大量由第一性原理计算的界面黏结能曲线，建立了炸药/添加剂界面原子势库，可以描述(RDX、HMX、TATB)/(石墨、石蜡、氟聚物)等三十余种界面。该参数库主要采用Morse势描述界面原子间的吸引力，采用指数衰减函数描述排斥力

${\varphi _{{\rm{Morse}}}} = {D_0}\left[ {{{\rm{e}}^{ - 2\alpha (r - {r_0})}} - 2{{\rm{e}}^{ - \alpha (r - {r_0})}}} \right],\quad {\varphi _{{\rm{Exp}}}} = A{{\rm{e}}^{ - r/\rho }}$

(3)

式中：r表示原子间距离，D₀、r₀、α为Morse力场参数，A、ρ为指数力场参数。图1给出了计算得到的TATB/氟聚物界面原子势曲线。采用二体相互作用描述炸药/添加剂界面是因为该界面间没有化学成键现象，只有较弱的范德华相互作用，因此可以采用对势模型。经过检验，基于第一性原理黏结能曲线推导出的对势模型具有良好的移植性，能够正确描述样本库之外的界面构型能量曲面。

图 1 TATB/氟聚物界面原子势曲线

Figure 1. The potential curves for the TATB/fluoropolymer interface

下载: 全尺寸图片幻灯片

为解决传统分子力场在描述炸药热分解行为的缺陷，Goddard等^[20–21]基于键序理论和典型过渡态能垒曲线建立了化学反应力场 (REAXFF)，可以描述典型炸药的热分解反应过程及中间产物。键序(BO)的计算考虑了 ${\text{σ}}$ 键、 ${\text{π}}$ 键和 ${\text{ππ}}$ 键，计算公式如下

${\rm{B}}{{\rm{O}}_{ij}} = \exp \left[ {{p_{{\rm{bo}},1}}{{\left( {\frac{{{r_{ij}}}}{{r_0^{{\text{σ}}}}}} \right)}^{{p_{{\rm{bo}},2}}}}} \right] + \exp \left[ {{p_{{\rm{bo}},3}}{{\left( {\frac{{{r_{ij}}}}{{r_0^{{\text{π}}}}}} \right)}^{{p_{{\rm{bo}},4}}}}} \right] + \exp \left[ {{p_{{\rm{bo}},5}}{{\left( {\frac{{{r_{ij}}}}{{r_0^{{\text{ππ}}}}}} \right)}^{{p_{{\rm{bo}},6}}}}} \right]$

(4)

式中：p_bo,1～p_bo,6表示键序参数。基于键序建立的总能计算公式为

${E_{{\rm{total}}}} = {E_{{\rm{bond}}}} + {E_{{\rm{over}}}} + {E_{{\rm{under}}}} + {E_{{\rm{val}}}} + {E_{{\rm{pen}}}} + {E_{{\rm{tors}}}} + {E_{{\rm{conj}}}} + {E_{{\rm{vdW}}}} + {E_{{\rm{Coul}}}}$

(5)

REAXFF力场模拟得到的热分解终态产物和实验一致。但使用该力场模拟得到的HMX炸药常温热力学参数与实验有一定误差，如计算获得的弹性常数偏小、热膨胀系数偏大等。Goddard团队应用反应力场研究了典型炸药的热分解过程^[22]、热点形成过程^[23–24]、冲击起爆过程^[25]等。

2. 热力学参数计算

高能炸药的热力学参数(如弹性常数、热容、热膨胀系数等)研究对于炸药的宏观应用具有重要的意义，如状态方程构建、添加剂配方设计、弹药安全性评估等。利用第一性原理方法和分子动力学方法计算炸药的热力学参数是目前的主要研究手段。主要的算法介绍如下。

（1）DFT+vdW方法：计算炸药晶体结构、弹性常数、生成焓。

（2）简谐声子理论：计算自由能、热膨胀系数、热容、相变曲线。

（3）非谐声子理论：计算热导率。

（4）分子动力学模拟：计算弹性常数、热膨胀系数、状态方程。

（5）非平衡分子动力学：计算热导率、黏性系数。

（6）涨落理论：计算比热、弹性系数。

（7）线性响应理论：计算热导率、热耗散率、黏性系数、扩散系数。

PBX炸药的物性研究包括3部分：单晶热力学参数、界面热力学参数和复合炸药热力学参数。在算法方面，弹性常数通过对能量-应变曲线求二阶导数获得。热膨胀系数通过对晶格长度-温度的关系求一阶导数获得。比热通过NPT系综下求温度涨落获得，计算公式为

${C_V} = \frac{{k{T^2}}}{{\Delta {T^2}}}$

(6)

式中：C_V为定容比热容，k为玻尔兹曼常数，T为温度，ΔT为温度的涨落。雨贡纽关系通过求解雨贡纽方程获得

$E - {E_0} = \frac{1}{2}(P + {P_0})({V_0} - V)$

(7)

等熵线通过求解等熵方程获得

$\frac{{{\rm{d}}P}}{{{\rm{d}}T}} = {{{{\left( {\frac{{\partial H}}{{\partial T}}} \right)}_P}} \Bigg/ {\left[ {V - {{\left( {\frac{{\partial H}}{{\partial P}}} \right)}_T}} \right]}}$

(8)

式中：H表示焓，E表示能量，P表示压力，V表示体积。

PBX炸药由于结构的复杂性，通过分子动力学计算得到的能量-应变曲线不够光滑，需要先将它拟合成为二次曲线才能获得二阶导数，从而计算弹性常数。但用这种方法计算高压下PBX炸药弹性常数会带来较大误差，无法正确反映压力对弹性常数的影响。对PBX炸药高压弹性性能的研究可通过拟合等温线获得体模量随压力变化的关系。

另外，通过涨落理论计算炸药的等容比热容也存在一定误差，一般只能获得零压下的热容。采用(6)式计算得到的炸药高压热容因误差的关系看不出与零压下的差别。实际上由于压缩导致体积缩小会使声子的运动熵降低，从而降低热容。要获得热容与压力的关系，需要计算炸药在不同压缩率下的声子频谱，通过声子自由能公式推导热容。

目前的线性响应方法和声子理论在研究炸药高压热力学参数时都遇到一定困难。采用线性响应方法计算热力学参数需要对关联函数进行时间积分，由于分子动力学计算得到的关联函数存在较大涨落误差，计算结果难以反映物性随压力的变化趋势。采用声子理论能得到较准确的高压热力学参数，但目前计算声子谱的力常数理论用于高压条件下容易产生虚频，需要改进算法以降低虚频的数量。

在单晶热力学参数研究方面，肖鹤鸣等^[26-27]、肖继军等^[28–30]计算了大量新型炸药的分子结构、键能和生成焓等参数。Zerilli等^[31]采用密度泛函理论计算了HMX炸药的体模量、热膨胀系数、比热、格林爱森参数随温度和压力的变化规律。Valenzano等^[32]采用B3LYP密度泛函近似计算了TATB和太安(PETN)的弹性常数。Sewell等^[6]采用分子动力学方法计算了HMX的晶格常数、弹性常数和等温线。Smith等^[5]计算了TATB炸药在高压下的弹性常数。姬广富等^[33]采用COMPASS力场和分子动力学方法计算了HMX炸药 $\alpha$ 、 $\beta$ 、 $\delta$ 相的晶格常数和弹性常数随温度/压力的变化趋势。目前对炸药在常温、常压下的热力学参数已有较充分的研究，取得了与实验一致的计算结果；但是对高压下热力学参数的计算较为困难。现阶段能够计算的高压参数只有部分弹性常数，对高压区域的热容、热膨胀系数等参数的研究工作较少。

在炸药/添加剂界面力学性能研究方面，目前已报导了关于(HMX、TATB)/(石墨、石蜡、氟聚物)等界面拉伸力学性能的系统性研究工作^{[15–17, 34]}，获得了界面拉伸势垒和拉伸应力。研究方法是通过分子动力学模拟了界面的拉伸断裂过程，获得系统总能与拉伸距离的关系。拉伸能量曲线的最高点与最低点之差即为拉伸势垒，最大斜率即为拉伸应力。图2给出了HMX/F2312界面在拉伸断裂过程中的界面构型变化。肖鹤鸣等^[35–36]计算了各种聚合物分子在基底炸药上的吸附构型，获得了界面黏结能。通过上述研究发现，界面力学性能与添加剂成分、界面取向和缺陷分布都有关系。根据目前的研究结果，界面处的结构缺陷一般会降低拉伸应力，因为在拉伸过程中，缺陷处成为应力集中区域，容易首先发生断裂，然后裂纹再扩展到其他区域。

图 2 HMX/F2312界面的拉伸断裂过程

Figure 2. The tensile process for the HMX/F2312 interface

下载: 全尺寸图片幻灯片

在PBX炸药物性研究方面，肖鹤鸣等^[35–37]应用COMPASS力场计算了一系列PBX炸药的构型和弹性力学性能，提出了添加剂对PBX炸药力学性能的影响机制。Jaidann等^[38–39]采用分子动力学方法计算了以RDX和FOX-7为基的PBX炸药的界面吸附能和弹性系数。上述工作采用的模型是在基底炸药上吸附一个添加剂分子，然后计算整个体系的弹性系数。用这种模型可以定性比较不同添加剂对炸药力学性能的影响，从而选出最适合做黏接剂、增塑剂的添加剂类型。

更复杂的PBX炸药模型是添加剂包覆在多个基底炸药颗粒外形成的混合模型^[15–19]。基于多晶包覆模型计算了RDX、HMX、TATB基复合炸药的热力学参数^[15–19]，包括二元、三元和四元体系。图3(a)～图3(c)分别给出了RDX多晶、石蜡包覆RDX、F2311包覆RDX 3种PBX炸药的初始结构。这类多晶包覆模型和实际的PBX炸药微观构型比较接近，计算结果具有合理性。所计算的热力学参数包括体模量、弹性常数、热容、热膨胀系数、格林爱森参数、等温线、等熵线和雨贡纽线。通过大量的数值计算和理论分析发现，PBX炸药物性预测的难点在于炸药的总体物性不仅与组分有关，还与各组分之间的界面有关。例如加入的添加剂的热膨胀系数小于基底炸药时，在升温过程中炸药颗粒膨胀速度大于包覆层，在包覆层上会产生热应力，抑制颗粒的膨胀。因此，界面热应力显著地影响了PBX炸药的总体热膨胀性能。除了界面热应力之外，界面声阻抗匹配和界面折射效应对PBX炸药的冲击、耗散性能也有显著的影响。

图 3 (a) RDX多晶，(b) 石蜡包覆RDX，(c) F2311包覆RDX

Figure 3. (a) RDX polycrystal, (b) paraffin coated RDX, (c) F2311 coated RDX

下载: 全尺寸图片幻灯片

3. 耗散和输运性能研究

PBX炸药的耗散与输运性能包括热传导、热耗散、黏性、扩散等性能，是分析炸药在外界刺激作用下热耗散、热累积过程的关键参数，对评估炸药的安全性和起爆机制有重要意义。目前对炸药耗散与输运性能的理论计算工作难度较大，并且方法也有限。传统计算方法主要有3种：Green-Kubo方法^[40]、非简谐声子理论^[41–43]和非平衡分子动力学^[44]。Green-Kubo方法基于线性响应理论，通过对关联函数积分获得热导率，可计算晶体材料的各向异性热导张量。非简谐声子理论基于量子多体理论，需要第一性原理方法计算三阶力常数，通过声子输运模型计算热导率，同样适用于晶体结构。对炸药来说，由于单晶原子数较多，获得三阶力常数需要完成数万个构型的电荷自洽计算，计算量很大。非平衡分子动力学按一定频率交换原子动量获得热流，通过分子动力学模拟获得温度梯度，两者相除即为热导，计算量较小。并且非平衡分子动力学方法可获得单晶和界面系统的热导率。

Smith等^{[7, 45]}采用非平衡分子动力学方法计算了HMX炸药的热导率。Sewell等^{[8–9, 46]}采用非平衡分子动力学方法计算了TATB炸药的各向异性热导率。非平衡分子动力学方法可计算不同晶格取向的热传导率，但仅限于热导张量的对角项。非对角项只能通过声子输运模型或Green-Kubo方法得到。

陈军等针对PBX炸药的热导性能开展了较为系统的研究工作，发展了多种方法计算单晶、界面和PBX炸药的热导率。针对含能材料分子晶体的结构特征，发展了利用二阶力常数计算声子寿命的方法，其思路是声子热应力改变了材料的声速，声速影响德拜频率，从而反过来改变声子频率，由此可推导出声子-声子相互作用模型，获得声子寿命。基于该方法获得的声子寿命公式为

$\tau = \frac{1}{{\omega \left| {\lambda + \mu E\left( {\dfrac{{h{\omega _0}}}{{2{\text{π}}kT}}} \right)} \right|}},\quad E(x) = \frac{9}{{{x^4}}}\int_0^x {\dfrac{{{y^3}}}{{{{\rm{e}}^y} - 1}}{\rm{d}}y}$

(9)

式中： $\tau$ 为声子寿命， $\omega$ 为声子频率，k为玻尔兹曼常数，h为普朗克常数，T为温度， $\lambda$ 和 $\mu$ 为材料的特征物理参数。采用该方法计算了HMX、TATB炸药的各向异性热导张量^[47–50]。并且声子输运的方法可解析出每种振动模式对热导的贡献，从而获得对热导起关键作用的分子振动模式。对于TATB炸药，总共发现二十余种关键分子振动模式，图4给出了其中两种模式的振动图谱。

图 4 对TATB热导起关键作用的部分分子振动模式：(a) TATB分子结构，(b) 振动模式1，(c) 振动模式2

Figure 4. Some key vibrational modes for thermal conduction in TATB: (a) TATB molecular structure, (b) the first vibrational mode, (c) the second vibrational mode

下载: 全尺寸图片幻灯片

针对炸药/添加剂界面体系，基于线性响应理论和力常数理论发展了两种计算界面声子折射率的方法^[50–54]，并基于声子折射率推导了界面热导率的计算公式

$\lambda = \frac{1}{{k{T^2}V}}\sum\limits_{{{ c}_i} \cdot {{n}} > 0} {{g_i}{\hbar ^2}\omega _i^2{P_i}{{ c}_i} \cdot {{n}}\frac{{{{\rm{e}}^{{{\hbar {\omega _i}}/ {kT}}}}}}{{{{\left( {{{\rm{e}}^{{{\hbar {\omega _i}}/ {kT}}}} - 1} \right)}^2}}}}$

(10)

式中： $\lambda$ 表示界面热导率，V表示晶胞体积，n表示界面法向， ${\omega _i}$ 表示声子频率，g_i表示声子简并度，c_i表示声子速度，P_i表示声子的界面折射率。界面热导率研究的关键是界面声子折射率的计算，这一步通过力常数理论或线性响应理论均可完成。后续的声子输运公式是通用的。应用上述方法研究了(HMX、TATB)/(石墨、石蜡、氟聚物)界面，计算了界面热导率与温度的关系，如图5所示。基于均匀热流模型，在单晶热导率和界面热导率的研究基础上推导了PBX炸药的总热导，发现PBX炸药热导与颗粒粒度有关。总热导率( ${\kappa _{{\rm{PBX}}}}$ )的计算公式为

图 5 TATB/添加剂界面热导率的计算结果

Figure 5. The interfacial thermal conductivity vs. temperature for TATB/additive interfaces

下载: 全尺寸图片幻灯片

${\kappa _{{\rm{PBX}}}} = \frac{{{\kappa _{\rm{I}}}{\kappa _{\rm{C}}}d}}{{{\kappa _{\rm{I}}}d + 2{\kappa _{\rm{C}}}}}\quad\quad\quad\quad\quad\quad\quad(11)\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\;$

(11)

式中： ${\kappa _{\rm{C}}}$ 为基底炸药的热导率， ${\kappa _{\rm{I}}}$ 为炸药/添加剂界面热导率，d为颗粒直径。

在热耗散率研究方面，Smith等^[52]采用线性响应方法计算了HMX炸药单晶黏性系数与温度的关系。单晶黏性系数( $\eta$ )的计算公式为

$\eta = \frac{V}{{10kT}}\int_0^\infty {\sum\limits_{i,j} {{q_{ij}}\left\langle {{\sigma _{ij}}(t){\sigma _{ij}}(0)} \right\rangle } {\rm{d}}t}$

(12)

式中： ${{\rm{\sigma }}_{ij}}$ 为应力分量。根据计算结果，获得了含能材料黏性系数与温度的关系

$\eta = {\eta _0}{{\rm e}^{ - \Delta E/kT}}$

(13)

界面热耗散效应的研究主要基于线性响应理论，发现不同加载强度下PBX炸药具有不同的热耗散机制：低强度加载下，界面热耗散主要通过弹性波散射机制完成；中等强度加载下，界面热耗散通过界面摩擦造成；强冲击作用下，热耗散通过界面塑性功完成。基于上述机制推导了界面黏性系数、界面摩擦系数和弹性波散射率的计算公式^[53]，并应用于HMX/(TATB、石墨、石蜡、氟聚物)界面体系。界面摩擦系数(η)的计算公式如下

$\eta=\frac{k T}{S \int_{0}^{\infty}\langle{\text{δ}}v(0) {\text{δ}}v(\tau)\rangle \mathrm{d} \tau}$

(14)

式中：S表示界面面积， ${{\text{δ}}v}$ 表示界面两侧的相对质心速度。通过分析弹性波耗散率与频率的关系发现，高频条件下热耗散系数与分子振动频率有关，低频条件下热耗散系数与频率的平方成正比。

4. 相图和相变动力学研究

含能材料是分子晶体结构，在加压、加温条件下分子会发生形变从而产生大量新相，如HMX具有4种常见晶型，CL-20具有6种常见晶型，TATB、RDX在加压至5～7 GPa时产生高压相结构。相变发生后，颗粒体积的变化会在PBX炸药内部造成孔洞、裂纹等缺陷；同时冲击加载条件下的高压相变会在冲击波波前形成相变波。由于炸药不同晶相之间的感度差异较大，并且相变可能导致相变潜热、裂纹和相变波，因此有关炸药相结构和相图的研究对评估炸药的安全性和武器系统的可靠性具有重要意义。

在相结构研究方面，Kholod等^[54]采用密度泛函理论计算了CL-20各晶相的拉曼谱和红外谱。肖鹤鸣等^[55]计算了CL-20的4种晶相的声子能带和声子态密度分布。Brand等^[56]采用第一性原理方法计算了HMX $\alpha$ 、 $\beta$ 、 $\delta$ 相的特征振动频谱。Munday等^[57]采用分子动力学方法研究了RDX的高压相结构。魏冬青等^[58]采用分子动力学模拟得到了HMX的高压相变点。目前实验上对RDX、HMX等典型炸药高压相的研究主要基于拉曼谱，通过拉曼谱的峰值变化确定新相是否产生。对RDX在单轴加载和动加载条件下找到了很多新相，但无法确定新相的分子结构。需要大量理论计算工作的配合，来预测实验获得的相结构。

在相变动力学研究方面，Smilowitz等^[59–63]建立了HMX炸药的相变速率唯象模型，获得了 $\beta$ 相到 $\delta$ 相的相变曲线。Levitas等^[61]基于熔化成核理论提出了描述炸药固-固相变过程的物理模型，根据该模型获得的相界面扩散速度(v)为

$v = {v_0}({{\rm{e}}^{ - \Delta {E_1}/RT}} - {{\rm{e}}^{ - \Delta {E_2}/RT}})$

(15)

式中：T为相变点温度， $\Delta {E_1}$ 为正向相变活化能， $\Delta {E_2}$ 为反向相变活化能。理论预测的相变速率与实验一致。Brill等^[62]提出了CL-20固-固相变的速率模型，并计算了相变活化能。Sewell等^[63]采用分子动力学方法计算了TATB炸药的固液相变曲线，同时计算了液相的黏性系数和热导率。对炸药相变动力学的研究基本是基于阿雷尼乌斯公式，通过拟合实验数据获得相变活化能。目前通过过渡态理论直接计算相变活化能的工作较少。由于相变过程和化学反应过程的相似性，实际可以将化学反应研究中的过渡态方法应用于相变研究中。

在冲击相图和相变波研究方面，目前针对含能材料体系发展了较为精确的电子-声子自由能模型^[64–68]，应用该模型计算了HMX炸药在高温、高压区域的准静态相变和冲击相变曲线。炸药的电子-声子自由能计算公式如下

$F(V,T) = {E_{\rm{c}}}(V) - kT\ln {Z_{\rm{m}}}(T) - kT\ln {Z_{\rm{c}}}\left[ {T{{\left( {\frac{V}{{{V_0}}}} \right)}^{\frac{1}{3}}}} \right] - kT\ln {Z_\chi }\left[ {T{{\left( {\frac{V}{{{V_0}}}} \right)}^{\frac{n}{3}}}} \right] - kT\ln {Z_{\rm{e}}}(T,\mu ) + {N_{\rm{e}}}\mu$

(16)

式中：F为自由能，V为体积，T为温度，E_c为冷能，Z_m为分子振动配分函数，Z_c为晶格振动配分函数， ${Z_\chi }$ 为高阶声子配分函数，Z_e为电子配分函数， $\mu$ 为费米能级，N_e为电子数。准静态相变点通过对比两相的吉布斯自由能获得，冲击相变点通过对比两相在相同飞片撞击速度下的熵获得。图6给出了基于德拜声子理论计算得到的HMX炸药相图。实际计算表明炸药的准静态相变曲线和冲击相变曲线不重合，这是分析相变对炸药冲击起爆的影响机制时需要注意的。

图 6 通过德拜声子理论计算得到的HMX的三相相图（

${\beta '}$ 表示

${\beta }$ 相通过声子软化形成的新相）

Figure 6. The phase diagram of HMX, obtained by Debye theory (The

${\beta '}$ phase is the phonon soften state of

${\beta }$ phase.)

下载: 全尺寸图片幻灯片

另外，陈军等^[65]应用爆轰物理领域中CJ理论的物理思想研究相变波。由于不同固体相之间的界面天然地不连续，因此由相变传递造成的相变波是间断波，符合定常流体力学方程的间断解，可用雨贡纽方程表征。以 $\beta$ 相作为起始相，通过以下公式求解三相的雨贡纽线

$\left\{ \begin{aligned} & {E_\alpha } - {E_0} = \frac{1}{2}({P_\alpha } + {P_0})({V_0} - {V_\alpha })\\ & {E_\beta } - {E_0} = \frac{1}{2}({P_\beta } + {P_0})({V_0} - {V_\beta })\\ & {E_\delta } - {E_0} = \frac{1}{2}({P_\delta } + {P_0})({V_0} - {V_\delta }) \end{aligned} \right.\quad\quad\quad(17)\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\;\;$

(17)

式中：E₀、P₀、V₀表示初始状态下 $\beta$ -HMX的能量、压力和体积。根据CJ理论，从初始状态(O点)出发的瑞利线与 $\alpha$ 、 $\delta$ 相雨贡纽线的切点代表相变波的稳定传播状态。图7给出了HMX三相的雨贡纽线和CJ点(A、B、C、D点)。根据理论分析，原点(O)左边的CJ点(A、B)表示冲击作用下的 $\beta$ → $\alpha$ 和 $\beta$ → $\delta$ 相变波，对应的相变波波速(D_p)分别为3796 m/s和4625 m/s。原点右边的CJ点(C、D)表示自发相变的 $\beta$ → $\alpha$ 和 $\beta$ → $\delta$ 相变波，对应的相变波波速分别为1800 m/s和1531 m/s。相变波波速通过瑞利线的斜率计算得到。计算结果表明，冲击作用下形成的相变波波速略高于炸药声速，非冲击条件下形成的相变波波速低于炸药声速。

图 7 HMX三相的雨贡纽线和瑞利线

Figure 7. The Hugoniot curves and Rayleigh curve for the three phases of HMX

下载: 全尺寸图片幻灯片

5. 动力学响应行为和本构关系研究

炸药的动力学响应行为指外界刺激作用下炸药表现出的应力弛豫现象。对PBX炸药来说，由于压制成型过程中产生了大量内应力，应力弛豫现象可能导致局部内应力的释放，通过热-功转换形成热点，产生安全性问题。同时，应力弛豫关系可用来构造炸药的黏弹性本构方程，为宏观有限元模拟提供模型和参数。

在炸药的动力学响应行为研究方面，目前根据声子数重分布机制和分子旋转机制发展了两种物理模型描述典型炸药的动力学响应行为^[66]。两种机制分别介绍如下：在正应力作用下，炸药晶体的声子谱发生变化，故声子数需重分布以适应新的态密度，声子重分布所需时间造成了声子热应力响应的滞后，即表现为应力-应变响应效应；在剪切作用下，分子需发生旋转，由于转动惯量和分子间摩擦力的存在，该旋转相对于应变有一定滞后，使剪切应力表现出应力-应变响应行为。因此上述两种机制代表了正应变和剪切应变的动力学响应关系。基于上述物理机制，正应力的应力弛豫关系为

$P(t) = {\lambda _1}\int_0^t {{{\rm{e}}^{ - (t - \tau )/{\tau _1}}}\varepsilon (\tau )} {\rm{d}}\tau + {\lambda _2}\int_0^t {{{\rm{e}}^{ - (t - \tau )/{\tau _2}}}\varepsilon (\tau )} {\rm{d}}\tau + {\lambda _3}\int_0^t {{{\rm{e}}^{ - (t - \tau )/{\tau _3}}}\varepsilon (\tau )} {\rm{d}}\tau + \sum\limits_{n = 0}^\infty {{\mu _n}\frac{{{{\rm{d}}^n}\varepsilon }}{{{\rm{d}}{t^n}}}}$

(18)

式中：P为正应力，t为时间， $\varepsilon$ 为正应变， ${\lambda _i}$ （i=1, 2, 3）为应力弛豫系数， ${\tau _i}$ 为应力弛豫时间， ${\mu _n}$ 为高阶黏性系数。剪切应力的动力学响应关系为

$\sigma (t) = \varepsilon \gamma + {C_1}\int_0^t {{{\rm{e}}^{ - (t - \tau )/{\eta _1}}}\varepsilon (\tau )} {\rm{d}}\tau + {C_2}\int_0^t {{{\rm{e}}^{ - (t - \tau )/{\eta _2}}}\varepsilon (\tau )} {\rm{d}}\tau + {C_3}\int_0^t {{{\rm{e}}^{ - (t - \tau )/{\eta _3}}}\varepsilon (\tau )} {\rm{d}}\tau$

(19)

式中： $\sigma$ 表示剪切应力， $\gamma$ 表示剪切形变弹性系数，C_i表示应力弛豫系数， ${\eta _i}$ 表示剪切应力的弛豫时间。计算得到的HMX炸药应力弛豫时间随温度的升高而降低。

在本构关系研究方面，Jones、Wilkins和Lee建立的JWL状态方程^[67]被广泛应用于描述固体炸药和爆轰产物。该方程的等熵形式为

$\left\{\begin{aligned} & P = A{{\rm{e}}^{ - {R_1}V/{V_0}}} + B{{\rm{e}}^{ - {R_2}V/{V_0}}} + C{\left( {\frac{V}{{{V_0}}}} \right)^{ - \omega - 1}} \\ & E = \frac{A}{{{R_1}}}{{\rm e}^{ - {R_1}V/{V_0}}} + \frac{B}{{{R_2}}}{{\rm e}^{ - {R_2}V/{V_0}}} + \frac{C}{\omega }{\left( {\frac{V}{{{V_0}}}} \right)^{ - \omega }} \end{aligned}\right.$

(20)

式中：A、B、C、R₁、R₂和 $\omega$ 为状态方程参数，P为等熵条件下压力，E为等熵条件下能量，V为体积。朱建士等建立的朱-王-唐本构模型^[68–72]被用于研究炸药的冲击损伤过程。该模型的本构方程为

$\sigma = {C_0}\varepsilon + \alpha {\varepsilon ^2} + \beta {\varepsilon ^3} + {C_1}\int_0^t {{{\rm{e}}^{ - (t - \tau )/{\theta _1}}}\varepsilon (\tau )} {\rm{d}}\tau + {C_2}\int_0^t {{{\rm{e}}^{ - (t - \tau )/{\theta _2}}}\varepsilon (\tau )} {\rm{d}}\tau$

(21)

式中： $\sigma$ 为应力， $\varepsilon$ 为应变，C₀、C₁、C₂、 $\alpha$ 、 $\beta$ 为应力弛豫系数， ${\theta _1}$ 、 ${\theta _2}$ 为应力弛豫时间。赵艳红等^[70–71]研究了炸药爆轰产物的状态方程，提出了产物状态方程的混合规则，能够描述多种产物组成的混合体系。杨明理等^[72]研究了爆轰产物的相分离机制，计算了多相产物体系的热力学参数和输运系数，为状态方程建模提供了参数。目前对未反应炸药和爆轰产物的状态方程研究较多，但对爆轰的中间产物、碳聚物的状态方程研究较少。实际上，现在对爆轰中间产物的种类、分子结构和碳聚物的相图缺少系统的科学认识，需要加强这方面的研究工作。

6. 热点形成机制和降感机理研究

炸药的缺陷演化和热点形成过程是理解冲击起爆微观机制的重要途径，受到国内外科研人员的广泛关注。PBX炸药内部的典型缺陷包括孔洞、位错、剪切带、裂纹等。外界刺激作用下不同缺陷产生耦合演化行为，表现出复杂的力学响应行为。如：强冲击作用下，炸药内部孔洞塌缩，并在塌缩点周围产生微裂纹；孔洞塌缩过程和裂纹扩展过程相互耦合，形成复合缺陷。近年来对PBX炸药缺陷演化和热点形成过程的研究已成为含能材料领域的前沿方向。

Menikoff^[73–74]应用欧拉-拉格朗日方法模拟了冲击条件下PBX炸药的孔洞塌缩和热点形成过程，提出了热点质量这一概念。Austin等^[75–76]模拟了HMX炸药的孔洞塌缩过程，发现孔洞附近产生了微裂纹。Springer等^[77]模拟了强冲击作用下的孔洞塌缩过程，获得了由孔洞塌缩造成的热点形貌和温度场分布。Goddar等^[23–24]模拟了复合炸药的冲击起爆过程，发现界面处容易形成热点。周婷婷等^[78]采用分子动力学方法研究了典型炸药的孔洞塌缩过程，获得了热点温度随时间演化的规律。

陈军等^[79–81]研究了PBX炸药的缺陷演化和热点形成机制，提出冲击条件下PBX炸药的结构缺陷附近容易产生塑性功，使得局部快速升温形成热点。主要的塑性功机制包括3种：界面摩擦、包覆层扭曲和孔洞塌缩。其中孔洞塌缩导致的热点升温速率最快。炸药塑性功( $\Delta E$ )的计算思路为总能减去冲击绝热压缩的能量，剩下部分即为塑性功

$\Delta E({ x}) = U({ x}) - E(u({ x}))$

(22)

式中：x表示位置；U表示总能分布函数；u表示流体速度分布；E为冲击绝热压缩的能量，通过雨贡纽关系获得。基于塑性功研究了热点的形成和演化过程，发现热点温度场的径向分布符合Bessel函数，随时间的演化符合指数关系。温度场方程可写为

$T(r,t) = {T_0} + T{{\rm e}^{ - t/\eta }}{J_0}\left( {r\sqrt {\frac{{\rho C}}{{\kappa \eta }}} } \right)$

(23)

式中： $\eta$ 为热点弛豫时间， $\rho$ 为材料密度，C为热点热容， $\kappa$ 为热点热导，J₀为零阶Bessel函数。通过(23)式获得的热点半径(R)与热点弛豫时间( $\eta$ )的关系为

$R = {x_0}\sqrt {\frac{{\kappa \eta }}{{\rho {C_{\rm{W}}}}}}$

(24)

式中：x₀为Bessel函数的零点。因此热点弛豫时间与热点半径的平方成正比。

另外，通过研究典型炸药缺陷演化的细观机制发现孔洞在强冲击作用下会演变为微裂纹^[79–81]，如图8所示。微裂纹的分布比较有规律，呈现“X”形，与孔洞的初始形状和冲击速度无关。通过分析孔洞塌缩初期的稀疏波提出了微裂纹的产生机制：孔洞塌缩过程中水平方向产生稀疏波，垂直方向产生冲击波，稀疏波和冲击波交叠的区域产生强剪切应力，导致该方向出现微裂纹。孔洞塌缩和裂纹扩展造成的塑性功是热点产生的重要机制。

图 8 孔洞塌缩后的场分布：(a) 温度场，(b) 密度场

Figure 8. The temperature field (a) and density field (b) after pore collapsing

下载: 全尺寸图片幻灯片

在炸药降感机制研究方面，张朝阳^[82–83]研究了添加剂对炸药的降感机制，提出石墨、石蜡等添加剂对HMX炸药颗粒的表面有润滑作用，从而降低了机械感度。陈军等^[84]分析了冲击后PBX炸药内部的能量场分布，基于能量禁锢机制提出了层状、链状添加剂对基底炸药的降感机制。基本思路是层状添加剂在PBX炸药中由于界面张力的原因具有较高的声速，炸药受到冲击后，添加剂层中的冲击波能量扩散速度较快，结合炸药/添加剂界面的折射效应分流了炸药颗粒中的冲击能量，因此降低了感度。基于上述机制，添加剂的降感能力可通过界面能禁锢率表征。该禁锢率的物理意义是冲击条件下被添加剂包覆层分流的冲击波能量。界面能量禁锢率( $\lambda$ )的计算公式为

$\lambda = \int_0^{\textstyle\min \left\{ {\dfrac{{\rm{\pi }}}{2},\arcsin \sqrt {\dfrac{{Yd{\rho _2}}}{{\sigma {\rho _1}}}} } \right\}} {\frac{{\sin \theta \sin 2\theta \sin \left( {\arctan \dfrac{{\sigma \sin \theta }}{{\sqrt {{C_{33}}d} \sqrt {\dfrac{{Yd{\rho _2}}}{{{\rho _1}}} - \sigma {{\sin }^2}\theta } }}} \right)}}{{\sin \left( {\theta + \arctan \dfrac{{\sigma \sin {\theta}}}{{\sqrt {{C_{33}}d} \sqrt {\dfrac{{Yd{\rho _2}}}{{{\rho _1}}} - \sigma {{\sin }^2}\theta } }}} \right)}}{\rm{d}}\theta }$

(25)

式中： $\theta$ 为入射角， ${\rho _1}$ 、 ${\rho _2}$ 为界面两侧材料密度，Y为基底炸药杨氏模量，C₃₃为添加剂剪切模量， $\sigma$ 为界面张力。图9给出了对一系列HMX/(石墨、石蜡、氟聚物)界面的计算结果。

图 9 HMX/添加剂界面能量禁锢率和界面张力的关系。

Figure 9. The energy constraint rate vs. interfacial tension curves for HMX/additive interfaces

下载: 全尺寸图片幻灯片

7. 结　论

目前，对PBX炸药的结构和热力学性能已有较充分的研究，系统计算了复合炸药的弹性常数、状态方程和雨贡纽关系；但对PBX炸药的耗散、输运性能研究较少，对炸药热导率、热耗散率的大部分理论计算工作还局限于单晶或熔融态。需要对界面和复合炸药系统的耗散和输运性能开展进一步的研究工作。

在相图和相变动力学研究方面：对炸药的高压相变研究较多，因为高压相变通过焓判定，较容易计算；对高温区域的相结构变化由于涉及精确的吉布斯自由能计算，理论工作较少。需要针对炸药分子晶体体系开发高置信度的自由能算法。在相变动力学研究方面：实验工作比较多，通过差式热计量仪能够测量相变速率，但理论建模工作较少。可将化学反应理论中的过渡态方法用于研究炸药的相变过程，通过理论计算获得相变速率模型，并与实验进行比较。

在损伤演化和热点形成机制方面，目前已有大量数值模拟工作研究冲击作用下缺陷的演化过程和热点的产生；但对结构缺陷在炸药爆轰反应后期的形态和表征研究较少。在炸药化学反应阶段，初始的缺陷可能会转换形态，变成爆轰波的局部畸变区，对PBX炸药的爆轰延迟和非理想爆轰现象构成影响。结构缺陷对炸药爆轰波形态的影响机制是爆轰物理学研究中的难点问题，需要建立新的理论框架和物理模型。

总的来说，对PBX炸药的结构和静力学行为已开展大量工作，但对炸药的耗散/输运性能、相变动力学、动力学响应行为、损伤演化机制等方面研究尚不充分，需要在理论算法和物理模型上进一步开展工作。

图 1 TATB/氟聚物界面原子势曲线

Figure 1. The potential curves for the TATB/fluoropolymer interface

下载: 全尺寸图片幻灯片

图 2 HMX/F2312界面的拉伸断裂过程

Figure 2. The tensile process for the HMX/F2312 interface

下载: 全尺寸图片幻灯片

图 3 (a) RDX多晶，(b) 石蜡包覆RDX，(c) F2311包覆RDX

Figure 3. (a) RDX polycrystal, (b) paraffin coated RDX, (c) F2311 coated RDX

下载: 全尺寸图片幻灯片

图 4 对TATB热导起关键作用的部分分子振动模式：(a) TATB分子结构，(b) 振动模式1，(c) 振动模式2

Figure 4. Some key vibrational modes for thermal conduction in TATB: (a) TATB molecular structure, (b) the first vibrational mode, (c) the second vibrational mode

下载: 全尺寸图片幻灯片

图 5 TATB/添加剂界面热导率的计算结果

Figure 5. The interfacial thermal conductivity vs. temperature for TATB/additive interfaces

下载: 全尺寸图片幻灯片

图 6 通过德拜声子理论计算得到的HMX的三相相图（ ${\beta '}$ 表示 ${\beta }$ 相通过声子软化形成的新相）

Figure 6. The phase diagram of HMX, obtained by Debye theory (The ${\beta '}$ phase is the phonon soften state of ${\beta }$ phase.)

下载: 全尺寸图片幻灯片

图 7 HMX三相的雨贡纽线和瑞利线

Figure 7. The Hugoniot curves and Rayleigh curve for the three phases of HMX

下载: 全尺寸图片幻灯片

图 8 孔洞塌缩后的场分布：(a) 温度场，(b) 密度场

Figure 8. The temperature field (a) and density field (b) after pore collapsing

下载: 全尺寸图片幻灯片

图 9 HMX/添加剂界面能量禁锢率和界面张力的关系。

Figure 9. The energy constraint rate vs. interfacial tension curves for HMX/additive interfaces

下载: 全尺寸图片幻灯片

参考文献(84)

[1]	SORESCU D C, RICE B M, THOMPSON D L. A transferable intermolecular potential for nitramine crystals [J]. The Journal of Physical Chemistry A, 1998, 102(43): 8386–8392. doi: 10.1021/jp9820525
[2]	SORESCU D C, RICE B M, THOMPSON D L. Isothermal-isobaric molecular dynamics simulations of 1, 3, 5, 7-tetranitro-1, 3, 5, 7-tetraazacyclooctane (HMX) crystals [J]. The Journal of Physical Chemistry A, 1998, 102(35): 6692–6695. doi: 10.1021/jp981661+
[3]	SMITH G D, BHARADWAJ R K. Quantum chemistry based force field for simulations of HMX [J]. The Journal of Physical Chemistry B, 1999, 103(4): 3570–3575.
[4]	BEDROV D, AYYAGARI C, SMITH G, et al. Molecular dynamics simulations of HMX crystal polymorphs using a flexible molecule force field [J]. Journal of Computer-Aided Materials Design, 2001, 8: 77–85. doi: 10.1023/A:1020046817543
[5]	BEDROV D, BORODIN O, SMITH G, et al. A molecular dynamics simulation study of crystalline 1, 3, 5-triamino-2, 4, 6-trinitrobenzene as a function of pressure and temperature [J]. The Journal of Chemical Physics, 2009, 131: 224703. doi: 10.1063/1.3264972
[6]	SEWELL T D, MENIKOFF R, BEDROW D, et al. A molecular dynamics simulation study of elastic properties of HMX [J]. The Journal of Chemical Physics, 2003, 119(14): 7417–7426. doi: 10.1063/1.1599273
[7]	BEDROV D, SMITH D. Thermal conductivity of molecular fluids from molecular dynamics simulations: application of a new imposed-flux method [J]. The Journal of Chemical Physics, 2000, 113(18): 8080–8084. doi: 10.1063/1.1312309
[8]	KROONBLAWD M P, SEWELL T D. Theoretical determination of anisotropic thermal conductivity for crystalline 1, 3, 5-triamino-2, 4, 6- trinitrobenzene (TATB) [J]. The Journal of Chemical Physics, 2013, 139: 074503. doi: 10.1063/1.4816667
[9]	KROONBLAWD M P, SEWELL T D. Theoretical determination of anisotropic thermal conductivity for initially defect-free and defective TATB single crystals [J]. The Journal of Chemical Physics, 2014, 141: 184501. doi: 10.1063/1.4901206
[10]	GEE R, ROSZAK S, BALASUBRAMANIAN L, et al. Ab initio based force field and molecular dynamics simulations of crystalline TATB [J]. The Journal of Chemical Physics, 2004, 120(15): 7059–7066. doi: 10.1063/1.1676120
[11]	SONG H J, ZHANG Y G, LI H, et al. All-atom, non-empirical, and tailor-made force field for α-RDX from first principles [J]. RSC Advances, 2014, 4(76): 40518–40533. doi: 10.1039/C4RA07195F
[12]	SUN H. COMPASS: An ab initio force-field optimized for condensed-phase applications-overview with details on alkane and benzene compounds [J]. The Journal of Physical Chemistry B, 1998, 102: 7338–7364. doi: 10.1021/jp980939v
[13]	BUNTE S, SUN H. Molecular modeling of energetic materials: the parameterization and validation of nitrate esters in the COMPASS force field [J]. The Journal of Physical Chemistry B, 2000, 104(11): 2477–2489. doi: 10.1021/jp991786u
[14]	TERSOFF J. Empirical interatomic potential for carbon, with applications to amorphous carbon [J]. Physical Review Letters, 1988, 61(25): 2879–2882. doi: 10.1103/PhysRevLett.61.2879
[15]	LONG Y, LIU Y G, NIE F D, et al. Force-field derivation and atomistic simulation of HMX/graphite interface and polycrystal systems [J]. Communications in Theoretical Physics, 2012, 57(1): 102–114. doi: 10.1088/0253-6102/57/1/16
[16]	LONG Y, LIU Y G, NIE F D, et al. The force-field derivation and atomistic simulation of HMX-fluoropolymer mixture explosives [J]. Colloid & Polymer Science, 2012, 290(18): 1855–1866.
[17]	LONG Y, LIU Y G, NIE F D, et al. Force-field derivation and atomistic simulation of HMX-TATB-graphite mixture explosives [J]. Modelling and Simulation in Materials Science and Engineering, 2012, 20(6): 065010. doi: 10.1088/0965-0393/20/6/065010
[18]	LONG Y, CHEN J. The force-field derivation and application of explosive/additive interfaces [J]. Modelling and Simulation in Materials Science and Engineering, 2016, 24(7): 075013. doi: 10.1088/0965-0393/24/7/075013
[19]	LONG Y, CHEN J. Theoretical study of the interfacial force-field, thermodynamic property, and heat stress for plastic bonded explosives [J]. The Journal of Physical Chemistry C, 2017, 121(5): 2778–2788. doi: 10.1021/acs.jpcc.6b11203
[20]	STRACHAN A, DUIN A, GODDARD W, et al. Shock waves in high-energy materials: the initial chemical events in nitramine RDX [J]. Physical Review Letters, 2003, 91: 098301. doi: 10.1103/PhysRevLett.91.098301
[21]	LIU L, LIU Y, GODDARD W, et al. ReaxFF-lg: correction of the ReaxFF reactive force field for London dispersion, with applications to the equations of state for energetic materials [J]. The Journal of Physical Chemistry A, 2011, 115(40): 11016–11022. doi: 10.1021/jp201599t
[22]	ZHANG L, ZYBIN S V, DUIN A, et al. Carbon cluster formation during thermal decomposition of octahydro-1, 3, 5, 7-tetranitro-1, 3, 5, 7- tetrazocine and 1, 3, 5-triamino-2, 4, 6-trinitrobenzene high explosives from ReaxFF reactive molecular dynamics simulations [J]. The Journal of Physical Chemistry A, 2009, 113(40): 10619–10640. doi: 10.1021/jp901353a
[23]	AN Q, ZYBIN S V, GORRDAR W A, et al. Elucidation of the dynamics for hot-spot initiation at nonuniform interfaces of highly shocked materials [J]. Physical Review B: Condensed Matter, 2011, 84: 220101. doi: 10.1103/PhysRevB.84.220101
[24]	AN Q, GORRDAR W A, ZYBIN S V, et al. Highly shocked polymer bonded explosives at a nonplanar interface: hot-spot formation leading to detonation [J]. The Journal of Physical Chemistry C, 2013, 117(50): 26551–26561. doi: 10.1021/jp404753v
[25]	CHERUKARA M J, WOOD M A, KOBER E M, et al. Ultra-fast chemistry under non-equilibrium conditions and the shock to deflagration transition at the nanoscale [J]. The Journal of Physical Chemistry C, 2015, 119: 22008–22015. doi: 10.1021/acs.jpcc.5b05362
[26]	肖鹤鸣, 居学海. 高能体系中的分子间相互作用 [M]. 北京: 科学出版社, 2004.
[27]	肖鹤鸣, 许晓娟, 邱玲. 高能量密度材料的理论设计 [M]. 北京: 科学出版社, 2008.
[28]	肖继军, 朱卫华, 朱伟, 等. 高能材料分子动力学 [M]. 北京: 科学出版社, 2013.
[29]	XIAO J J, WANG W R, CHEN J, et al. Study on the relations of sensitivity with energy properties for HMX and HMX-based PBXs by molecular dynamics simulation [J]. Physica B: Condensed Matter, 2012, 407(17): 3504–3509. doi: 10.1016/j.physb.2012.05.010
[30]	CAO Q, XIAO J J, GAO P, et al. Molecular dynamics simulations for CL-20/TNT co-crystal based polymer-bonded explosives [J]. Journal of Theoretical and Computational Chemistry, 2017, 16: 1750072. doi: 10.1142/S0219633617500729
[31]	ZERILLI F J, KUKLJA M. Ab initio equation of state of the organic molecular crystal: β-Octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine [J]. The Journal of Physical Chemistry A, 2010, 114(16): 5372–5376. doi: 10.1021/jp911767q
[32]	VALENZANO L, SLOUGH W J, PERGER W. Accurate prediction of second-order elastic constants from first principles: PETN and TATB [J]. AIP Conference Proceedings, 2012, 1426: 1191–1194.
[33]	CUI H, JI G F, CHEN X, et al. Phase transitions and mechanical properties of octahydro-1, 3, 5, 7-tetranitro-1, 3, 5, 7- tetrazocine in different crystal phases by molecular dynamics simulation [J]. Journal of Chemical & Engineering Data, 2010, 55(9): 3121–3129.
[34]	LONG Y, LIU Y G, NIE F D, et al. Theoretical study of breaking and slipping processes for HMX/graphite interface [J]. Applied Surface Science, 2012, 258(7): 2384–2392. doi: 10.1016/j.apsusc.2011.10.052
[35]	QIU L, XIAO H M. Molecular dynamics study of binding energies, mechanical properties, and detonation performances of bicyclo-HMX-based PBXs [J]. Journal of Hazardous Material, 2009, 164(1): 329–336. doi: 10.1016/j.jhazmat.2008.08.030
[36]	XIAO J, HUANG H, LI J S, et al. A molecular dynamics study of interface interactions and mechanical properties of HMX-based PBXs with PEG and HTPB [J]. Journal of Molecular Structure: THEOCHEM, 2008, 851(1): 242–248.
[37]	XIAO J, ZHANG H, HUANG H, et al. NPT ensemble MD simulation investigation on the mechanical properties of HMX/F2311 polymer-bonded explosive [J]. Chinese Journal of Chemistry, 2008, 26(11): 1969–1972. doi: 10.1002/cjoc.v26:11
[38]	JAIDANN M, LUSSIER L S, BOUAMOUL A, et al. Effects of interface interactions on mechanical properties in RDX-based PBXs HTPB-DOA: molecular dynamics simulations [C]//Computational Science-ICCS 2009. Heidelberg: Springer-Berlin, 2009: 131-140.
[39]	JAIDANN M, LUSSIER L S, BOUAMOUL A, et al. Atomistic studies of RDX and FOX-7-based plastic-bonded explosives: molecular dynamics simulation [J]. Procedia Computer Science, 2011, 4: 1177–1185. doi: 10.1016/j.procs.2011.04.126
[40]	KUBO R, TODA M, HASHITSUME N. Statistical physics II: nonequilibrium statistical mechanics [M]. Springer-Verlag, 1997.
[41]	PARLINSKI K, LI Z Q, KAWAZOE Y. First-principles determination of the soft mode in cubic ZrO₂ [J]. Physical Review Letters, 1997, 78(21): 4063–4066. doi: 10.1103/PhysRevLett.78.4063
[42]	MARADUDIN A, FEIN A. Scattering of neutrons by an anharmonic crystal [J]. Physical Review, 1962, 128(6): 2589–2608. doi: 10.1103/PhysRev.128.2589
[43]	TOGO A, CHAPUT L, TANAKA I. Distributions of phonon lifetimes in Brillouin zones [J]. Physical Review B: Condensed Matter, 2015, 91: 094306. doi: 10.1103/PhysRevB.91.094306
[44]	MÜLLER-PLATHE F. Reversing the perturbation in nonequilibrium molecular dynamics: an easy way to calculate the shear viscosity of fluids [J]. Physical Review E, 1999, 59: 4894–4898.
[45]	BEDROV D, SMITH G D, SEWELL T D. Thermal conductivity of liquid octahydro-1, 3, 5, 7-tetranitro-1, 3, 5, 7-tetrazocine (HMX) from molecular dynamics simulations [J]. Chemical Physics Letters, 2000, 324: 64–68. doi: 10.1016/S0009-2614(00)00559-5
[46]	KROONBLAWD M P, SEWELL T D. Predicted anisotropic thermal conductivity for crystalline 1, 3, 5-triamino-2, 4, 6-trinitobenzene (TATB): temperature and pressure dependence and sensitivity to intramolecular force field terms [J]. Propellants, Explosives, Pyrotechnics, 2016, 41(3): 502–513. doi: 10.1002/prep.v41.3
[47]	LONG Y, CHEN J, LIU Y G, et al. A direct method to calculate thermal conductivity and its application in solid HMX [J]. Journal of Physics: Condensed Matter, 2010, 22: 185404. doi: 10.1088/0953-8984/22/18/185404
[48]	LONG Y, LIU Y G, NIE F D, et al. A method to calculate the thermal conductivity of HMX under high pressure [J]. Philosophical Magazine, 2012, 92(8): 1023–1045. doi: 10.1080/14786435.2011.637981
[49]	LONG Y, CHEN J. Theoretical study of the phonon-phonon scattering mechanism and the thermal conductive coefficients for energetic material [J]. Philosophical Magazine, 2017, 97: 2575–2595. doi: 10.1080/14786435.2017.1343962
[50]	LONG Y, CHEN J. A theoretical study of wave dispersion and thermal conduction for HMX/additive interfaces [J]. Modelling and Simulation in Materials Science and Engineering, 2014, 22: 035013. doi: 10.1088/0965-0393/22/3/035013
[51]	LONG Y, CHEN J. Theoretical study of the phonon spectrum, phonon refraction and thermodynamic properties for explosive/additive interfaces [J]. Modelling and Simulation in Materials Science and Engineering, 2018, 26: 015002. doi: 10.1088/1361-651X/aa944d
[52]	BEDROV D, SMITH G D, SEWELL T D. Temperature-dependent shear viscosity coefficient of octahydro-1, 3, 5, 7-tetranitro-1, 3, 5, 7-tetrazocine (HMX): a molecular dynamics simulation study [J]. The Journal of Chemical Physics, 2000, 112(16): 7203–7208. doi: 10.1063/1.481285
[53]	LONG Y, CHEN J. The heat dissipation model and desensitizing mechanism of the HMX/additive interfaces: a theoretical investigation based on linear response theory [J]. Modelling and Simulation in Materials Science and Engineering, 2013, 21: 055025. doi: 10.1088/0965-0393/21/5/055025
[54]	KHOLOD Y, OKOVYTYY S, KURAMSHINA G, et al. An analysis of stable forms of CL-20: a DFT study of conformational transitions, infrared and Raman spectra [J]. Journal of Molecular Structure, 2007, 843: 14–25. doi: 10.1016/j.molstruc.2006.12.031
[55]	XU X J, ZHU W H, XIAO H M. DFT Studies on the four polymorphs of crystalline CL-20 and the influences of hydrostatic pressure on ε-CL-20 crystal [J]. The Journal of Physical Chemistry B, 2007, 111(8): 2090–2097. doi: 10.1021/jp066833e
[56]	BRAND H V, RABIE R L, FUNK D J, et al. Theoretical and experimental study of the vibrational spectra of the α, β and δ phases of octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX) [J]. The Journal of Physical Chemistry B, 2002, 106(41): 10594–10604. doi: 10.1021/jp020909z
[57]	MUNDAY L B, CHUNG P W, RICE B M, et al. Simulations of high-pressure phases in RDX [J]. The Journal of Physical Chemistry B, 2011, 115(15): 4378–4386. doi: 10.1021/jp112042a
[58]	LU L Y, WEI D Q, CHEN X R. The pressure-induced phase transition of the solid β-HMX [J]. Molecular Physics, 2009, 107(22): 2373–2385. doi: 10.1080/00268970903313642
[59]	SMILOWITZ L, HENSON B F, ASAY B W, et al. The β-δ phase transition in the energetic nitramine-octahydro-1, 3, 5, 7-tetranitro-1, 3, 5, 7-tetrazocine: kinetics [J]. The Journal of Chemical Physics, 2002, 117(8): 3789–3798. doi: 10.1063/1.1495399
[60]	SMILOWITZ L, HENSON B F, ASAY B W, et al. The β-δ phase transition in the energetic nitramine-octahydro-1, 3, 5, 7-tetranitro-1, 3, 5, 7-tetrazocine: thermodynamics [J]. The Journal of Chemical Physics, 2002, 117(8): 3780–3788. doi: 10.1063/1.1495398
[61]	LEVITAS V I, HENSON B F, SMILOWITZ L B, et al. Solid-solid phase transformation via internal stress-induced virtual melting, significantly below the melting temperature. application to HMX energetic crystal [J]. The Journal of Physical Chemistry B, 2006, 110(20): 10105–10119. doi: 10.1021/jp057438b
[62]	LI J, BRILL T B. Kinetics of solid polymorphic phase transitions of CL-20 [J]. Propellants, Explosives, Pyrotechnics, 2007, 32(4): 326–330. doi: 10.1002/(ISSN)1521-4087
[63]	MATHEW N, KROONBLAWD M, SEWELL T D, et al. Predicted melt curve and liquid-state transport properties of TATB from molecular dynamics simulations [J]. Molecular Simulation, 2018, 44: 613–622. doi: 10.1080/08927022.2017.1418084
[64]	LONG Y, CHEN J. Theoretical study of phonon density of states, thermodynamic properties and phase transitions for HMX [J]. Philosophical Magazine, 2014, 94: 2656–2677. doi: 10.1080/14786435.2014.927598
[65]	LONG Y, CHEN J. Theoretical study of the thermodynamic properties, phase transition wave, and phase transition velocity for octahydro-1, 3, 5, 7- tetranitro-1, 3, 5, 7-tetrazocine [J]. Journal of Applied Physics, 2015, 118: 115901. doi: 10.1063/1.4930812
[66]	LONG Y, CHEN J. A theoretical study of the stress relaxation in HMX on the picosecond time scale [J]. Modelling and Simulation in Materials Science and Engineering, 2015, 23: 085001. doi: 10.1088/0965-0393/23/8/085001
[67]	KURY J W, HORNIG H C, LEE E L, et al. Metal accelaration by chemical explosives [C]//Proceedings of the 4th International Symposium on Detonation. Maryland: White Oak, 1966: 3–13.
[68]	WANG L L, ZHU X X, SHI S Q. An impact dynamics investigation on some problems in bird strike on windshield of high speed aircrafts [J]. Chinese Journal of Aeronautics, 1991, 12(3): 27–33.
[69]	周风华, 王礼立, 胡时胜. 有机玻璃在高应变率下的损伤型非线性粘弹性本构关系及破坏准则 [J]. 爆炸与冲击, 1992, 12(4): 333–341. ZHOU F H, WANG L L, HU S S. A damage-modified nonlinear visco-elastic constitutive relation and its failure criterion of PMMA at high strain rates [J]. Explosion and Shock Waves, 1992, 12(4): 333–341.
[70]	赵艳红, 刘海风, 张弓木. PETN炸药爆轰产物状态方程的理论研究 [J]. 高压物理学报, 2009, 23(2): 143–149. doi: 10.3969/j.issn.1000-5773.2009.02.011 ZHAO Y H, LIU H F, ZHANG G M. Equation of state of detonation products for PETN explosive [J]. Chinese Journal of High Pressure Physics, 2009, 23(2): 143–149. doi: 10.3969/j.issn.1000-5773.2009.02.011
[71]	赵艳红, 刘海风, 张广财. PBX9502炸药爆轰产物的状态方程 [J]. 爆炸与冲击, 2010, 30(6): 647–651. ZHAO Y H, LIU H F, ZHANG G C. Equation of state of detonation products for PBX9502 explosive [J]. Explosion and Shock Waves, 2010, 30(6): 647–651.
[72]	LÜ L, ZHANG L, YANG M L. Understanding the phase separation of N₂/H₂O and CO₂/H₂O binary systems through reactive force fields-based molecular dynamics simulations [J]. Journal of Applied Physics, 2018, 124: 235901. doi: 10.1063/1.5066585
[73]	MENIKOFF R. Compaction wave profiles in granular HMX [J]. AIP Conference Proceedings, 2002, 620: 979–982. doi: 10.1063/1.1483701
[74]	MENIKOFF R. Pore collapse and hot spots in HMX [J]. AIP Conference Proceedings, 2004, 706: 393–396. doi: 10.1063/1.1780261
[75]	AUSTIN R A, BARTON N R, HOWARD W M, et al. Modeling pore collapse and chemical reactions in shock-loaded HMX crystals [J]. Journal of Physics: Conference Series, 2014, 500: 052002. doi: 10.1088/1742-6596/500/5/052002
[76]	AUSTIN R A, BARTON N R, REAUGH J E, et al. Direct numerical simulation of shear localization and decomposition reactions in shock-loaded HMX crystal [J]. Journal of Applied Physics, 2015, 117: 185902. doi: 10.1063/1.4918538
[77]	SPRINGER H K, TARVER C M, BASTEA S. Effects of high shock pressures and pore morphology on hot spot mechanisms in HMX [J]. AIP Conference Proceedings, 2017, 1793: 080002. doi: 10.1063/1.4971608
[78]	ZHOU T T, LOU J F, ZHANG Y G, et al. Hot spot formation and chemical reaction initiation in shocked HMX crystals with nanovoids: a large-scale reactive molecular dynamics study [J]. Physical Chemistry Chemical Physics, 2016, 18(26): 17627–17645. doi: 10.1039/C6CP02015A
[79]	LONG Y, CHEN J. An investigation of the hot spot formation mechanism for energetic material [J]. Journal of Applied Physics, 2017, 122: 175105. doi: 10.1063/1.4996385
[80]	LONG Y, CHEN J. A molecular dynamics study of the early-time mechanical heating in shock-loaded octahydro-1, 3, 5, 7-tetranitro-1, 3, 5, 7- tetrazocine-based explosives [J]. Journal of Applied Physics, 2014, 116: 033516. doi: 10.1063/1.4890715
[81]	LONG Y, CHEN J. Theoretical study of the defect evolution for molecular crystal under shock loading [J]. Journal of Applied Physics, 2019, 125: 065107. doi: 10.1063/1.5067284
[82]	ZHANG C Y. Computational investigation on the desensitizing mechanism of graphite in explosives versus mechanical stimuli: compression and glide [J]. The Journal of Physical Chemistry B, 2007, 111(22): 6208–6213. doi: 10.1021/jp070918d
[83]	ZHANG C Y. Understanding the desensitizing mechanism of olefin in explosives versus external mechanical stimuli [J]. The Journal of Physical Chemistry C, 2010, 114(11): 5068–5072. doi: 10.1021/jp910883x
[84]	LONG Y, LIU Y G, NIE F D, et al. Theoretical study of impacting and desensitizing for HMX-graphite mixture explosive [J]. Shock Waves, 2012, 22(6): 605–614. doi: 10.1007/s00193-012-0394-7