高压下立方BC<sub>3</sub>的力学和热力学性质

摘要: 采用基于密度泛函理论的第一性原理赝势方法，系统地研究了立方BC₃在常压和高压下的晶格常数和力学性质，包括弹性常数、弹性模量和力学各向异性。利用准简谐近似下的德拜模型研究了高温高压条件下的热力学性质。研究结果表明：常压下立方BC₃具有较大的弹性模量和力学各向异性；高压下，立方BC₃的晶格常数、弹性常数和弹性模量显著增加。热力学性质的计算结果表明，立方BC₃具有较高的德拜温度，其摩尔定容热容和摩尔定压热容在高温高压条件下呈现明显的变化。立方BC₃的德拜温度随着压力的增大而增加，但随着温度的增大而明显减小。
- 第一性原理 /
- 立方BC₃ /
- 高压 /
- 力学性质 /
- 热力学性质
Abstract: The lattice constant and mechanical properties of cubic BC₃ under ambient and high pressure, including the elastic constants, the elastic modulus, and the mechanical anisotropy, were investigated using the first principle method in the framework of the density functional theory.The thermodynamic properties under high temperature and high pressure were calculated in terms of the quasi-harmonic Debye model.The results obtained show that the cubic BC₃ possesses a large elastic modulus and a high degree of anisotropy under ambient pressure.Under high pressure, the lattice constant, elastic constants, and elastic modulus of cubic BC₃ increase significantly.The results obtained from the thermodynamic calculations suggest that the cubic BC₃ has a large Debye temperature, and the molar heat capacity at constant volume and pressure exhibits obvious variation under high temperature and high pressure.Meanwhile, The Debye temperature of cubic BC₃ increases with the increase of pressure, but decreases with the increase of temperature.
- first principle /
- cubic BC₃ /
- high pressure /
- mechanical properties /
- 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 立方BC₃的晶体结构

Figure 1. Crystal structure of cubic BC₃

下载: 全尺寸图片幻灯片

图 2 立方BC₃归一化晶格常数a/a₀和密度ρ/ρ₀随压力的变化关系(a)以及能量-体积物态方程(b)

Figure 2. Pressure dependence of normalized lattice constants a/a₀ and density ρ/ρ₀ (a), and the equation of state for cubic BC₃ (b)

下载: 全尺寸图片幻灯片

图 3 立方BC₃弹性常数随压力的变化关系

Figure 3. Pressure dependence of elastic constants for cubic BC₃

下载: 全尺寸图片幻灯片

图 4 杨氏模量各向异性的三维图(a)和平面投影图(b)

Figure 4. Three-dimensional plot (a) and the corresponding projection (b) of anisotropy of Young's modulus

下载: 全尺寸图片幻灯片

图 5 立方BC₃的摩尔定容热容c_V、摩尔定压热容c_p和德拜温度随压力和温度的变化关系

Figure 5. Molar heat capacity c_V, c_p and Debye temperature vs. pressure and temperature for cubic BC₃

下载: 全尺寸图片幻灯片

表 1 不同压力下立方BC₃的体弹性模量、剪切模量和杨氏模量

Table 1. Bulk modulus, shear modulus, and Young's modulus for cubic BC₃ under pressure

Material	Pressure/GPa	B/GPa	G/GPa	Y/GPa
d-BC₃	0	34 5	318	730
	10	381	332	772
	20	416	343	807
	30	450	354	84 1
	4 0	483	363	870
	50	515	371	898
	60	54 7	379	923
	70	578	386	94 8
	80	610	393	971
	90	64 0	399	992
	100	671	4 05	1 012
c-BN^[28]	0	376	390
Diamond^[29-30]	0	432	517	1109

下载: 导出CSV

参考文献(33)

[1]	HAINES J, LÉGER J M, BOCQUILLON G.Synthesis and design of superhard materials[J]. Annual Review of Materials Research, 2001, 31(1):1-23. doi: 10.1146/annurev.matsci.31.1.1
[2]	TIAN Y, XU B, ZHAO Z.Microscopic theory of hardness and design of novel superhard crystals[J]. International Journal of Refractory Metals and Hard Materials, 2012, 33:93-106. doi: 10.1016/j.ijrmhm.2012.02.021
[3]	NOVIKOV N V.Synthesis of superhard materials[J]. Journal of Materials Processing Technology, 2005, 161(1):169-172. https://www.sciencedirect.com/science/article/pii/S0924013604009094
[4]	SOLOZHENKO V L, GREGORYANZ E.Synthesis of superhard materials[J]. Materials Today, 2005, 8(11):44-51. doi: 10.1016/S1369-7021(05)71159-7
[5]	NOVIKOV N V, DUB S N.Fracture toughness of diamond single crystals[J]. Journal of Hard Materials, 1991, 2(1):3-11. https://www.researchgate.net/publication/284697397_Fracture_toughness_of_diamond_single_crystals
[6]	SOLOZHENKO V L, DUB S N, NOVIKOV N V.Mechanical properties of cubic BC₂N, a new superhard phase[J]. Diamond and Related Materials, 2001, 10(12):2228-2231. doi: 10.1016/S0925-9635(01)00513-1
[7]	ZININ P V, MING L C, ISHⅡ H A, et al.Phase transition in BC_x system under high-pressure and high-temperature:synthesis of cubic dense BC₃ nanostructured phase[J]. Journal of Applied Physics, 2012, 111(11):114905. doi: 10.1063/1.4723275
[8]	SOLOZHENKO V L, KURAKEVYCH O O, ANDRAULT D, et al.Ultimate metastable solubility of boron in diamond:synthesis of superhard diamondlike BC₅[J]. Physical Review Letters, 2009, 102(1):015506. doi: 10.1103/PhysRevLett.102.015506
[9]	SOLOZHENKO V L, ANDRAULT D, FIQUET G, et al.Synthesis of superhard cubic BC₂N[J]. Applied Physics Letters, 2001, 78(10):1385-1387. doi: 10.1063/1.1337623
[10]	BADZIAN A R.Superhard material comparable in hardness to diamond[J]. Applied Physics Letters, 1988, 53(25):2495-2497. doi: 10.1063/1.100528
[11]	CHUNG H Y, WEINBERGER M B, LEVINE J B, et al.Synthesis of ultra-incompressible superhard rhenium diboride at ambient pressure[J]. Science, 2007, 316(5823):436-439. doi: 10.1126/science.1139322
[12]	GOU H, DUBROVINSKAIA N, BYKOVA E, et al.Discovery of a superhard iron tetraboride superconductor[J]. Physical Review Letters, 2013, 111(15):157002. doi: 10.1103/PhysRevLett.111.157002
[13]	CROWHURST J C, GONCHAROV A F, SADIGH B, et al.Synthesis and characterization of the nitrides of platinum and iridium[J]. Science, 2006, 311(5765):1275-1278. doi: 10.1126/science.1121813
[14]	KUMAR N R S, CHANDRA S, AMIRTHAPANDIAN S, et al.Investigations of the high pressure synthesized osmium carbide by experimental and computational techniques[J]. Materials Research Express, 2015, 2(1):016503. doi: 10.1088/2053-1591/2/1/016503
[15]	DOMNICH V, REYNAUD S, HABER R A, et al.Boron carbide:structure, properties, and stability under stress[J]. Journal of the American Ceramic Society, 2011, 94(11):3605-3628. doi: 10.1111/jace.2011.94.issue-11
[16]	GLASS C W, OGANOV A R, HANSEN N.Uspex-evolutionary crystal structure prediction[J]. Computer Physics Communications, 2006, 175(11):713-720.
[17]	LYAKHOV A O, OGANOV A R, STOKES H T, et al.New developments in evolutionary structure prediction algorithm uspex[J]. Computer Physics Communications, 2013, 184(4):1172-1182. doi: 10.1016/j.cpc.2012.12.009
[18]	WANG Y, LV J, ZHU L, et al.CALYPSO:a method for crystal structure prediction[J]. Computer Physics Communications, 2012, 183(10):2063-2070. doi: 10.1016/j.cpc.2012.05.008
[19]	ZHANG M, LIU H, LI Q, et al.Superhard BC₃ in cubic diamond structure[J]. Physical Review Letters, 2015, 114(1):015502. doi: 10.1103/PhysRevLett.114.015502
[20]	KRESSE G, FURTHMVLLER J.Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[J]. Physical Review B, 1996, 54(16):11169-11186. doi: 10.1103/PhysRevB.54.11169
[21]	PERDEW J P, BURKE K, ERNZERHOF M.Generalized gradient approximation made simple[J]. Physical Review Letters, 1996, 77(18):3865-3868. doi: 10.1103/PhysRevLett.77.3865
[22]	KRESSE G, JOUBERT D.From ultrasoft pseudopotentials to the projector augmented-wave method[J]. Physical Review B, 1999, 59(3):1758-1775. doi: 10.1103/PhysRevB.59.1758
[23]	HILL R.The elastic behaviour of a crystalline aggregate[J]. Proceedings of the Physical Society Section A, 1952, 65(5):349-354. doi: 10.1088/0370-1298/65/5/307
[24]	OTERO-DE-LA-ROZA A, LUAÑA V.Gibbs 2:a new version of the quasi-harmonic model code.Ⅰ.robust treatment of the static data[J]. Computer Physics Communications, 2011, 182(8):1708-1720. doi: 10.1016/j.cpc.2011.04.016
[25]	BIRCH F.Finite elastic strain of cubic crystals[J]. Physical Review, 1947, 71(11):809-824. doi: 10.1103/PhysRev.71.809
[26]	MOUHAT F, COUDERT F X.Necessary and sufficient elastic stability conditions in various crystal systems[J]. Physical Review B, 2014, 90(22):224104. doi: 10.1103/PhysRevB.90.224104
[27]	WU Z, ZHAO E, XIANG H, et al.Crystal structures and elastic properties of superhard IrN₂ and IrN₃ from first principles[J]. Physical Review B, 2007, 76(5):054115. doi: 10.1103/PhysRevB.76.054115
[28]	ZHANG R F, VEPREK S, ARGON A S.Anisotropic ideal strengths and chemical bonding of wurtzite BN in comparison to zincblende BN[J]. Physical Review B, 2008, 77(17):998-1002.
[29]	WANG Y J, WANG C Y.Mechanical properties and electronic structure of superhard diamondlike BC₅:a first-principles study[J]. Journal of Applied Physics, 2009, 106(4):043513. doi: 10.1063/1.3195082
[30]	ZHANG R F, LIN Z J, VEPREK S.Anisotropic ideal strengths of superhard monoclinic and tetragonal carbon and their electronic origin[J]. Physical Review B, 2011, 83(15):4400-4408.
[31]	CAZZANI A, ROVATI M.Extrema of Young's modulus for cubic and transversely isotropic solids[J]. International Journal of Solids and Structures, 2003, 40(7):1713-1744. doi: 10.1016/S0020-7683(02)00668-6
[32]	KLEIN C A.Anisotropy of Young's modulus and Poisson's ratio in diamond[J]. Materials Research Bulletin, 1992, 27(12):1407-1414. doi: 10.1016/0025-5408(92)90005-K
[33]	ZHENG B, ZHANG M, LUO H G.Pressure effect on structural, elastic, and thermodynamic properties of tetragonal B₄C₄[J]. AIP Advances, 2015, 5(3):436-439. http://www.osti.gov/scitech/biblio/22454476