物理化学学报, 2017, 33(5): 949-959 doi: 10.3866/PKU.WHXB201702152

论文

冲击载荷下TATB晶体滑移和各向异性的分子动力学研究

周婷婷,1, 宋华杰1, 黄风雷,2

The Slip and Anisotropy of TATB Crystal under Shock Loading via Molecular Dynamics Simulation

ZHOU Ting-Ting,1, SONG Hua-Jie1, HUANG Feng-Lei,2

收稿日期: 2016-10-10  

基金资助: 国家自然科学基金.  11402031
国家自然科学基金.  11372053
国家自然科学基金.  11221202

Corresponding authors: ZHOU Ting-Ting, Email: zhou_tingting@iapcm.ac.cnHUANG Feng-Lei, Email: huangfl@bit.edu.cn; Tel: +86-10-59872646

Received: 2016-10-10  

Fund supported: theNationalNaturalScienceFoundationofChina.  11402031
the National Natural Science Foundation of China.  11372053
the National Natural Science Foundation of China.  11221202

摘要

采用ReaxFF反应力场和分子动力学方法，研究了1，3，5-三氨基-2，4，6-三硝基苯（TATB）炸药晶体在沿不同方向冲击载荷下的滑移和各向异性。冲击方向分别垂直于（101）、（111）、（011）、（110）、（010）、（100）和（001）晶面，冲击强度为10 GPa。研究结果表明，各冲击方向下可能被激发的滑移系均在{001}面，而其它滑移系均因很大的剪切阻力不容易被激发，这与TATB晶体沿c轴的层状结构和平面分子结构相符。预测了七个冲击方向下最容易被激发的滑移系，分别为（101）/{001}＜100＞、（111）/{001}＜010＞、（011）/{001}＜010＞、（110）/{001}＜010＞、（010）/{001}＜1$\bar 1$0＞、（100）/{001}＜120＞和（001）/{001}＜010＞。TATB晶体的冲击响应具有各向异性，动力学过程中体系的应力、能量、温度和化学反应都依赖于冲击方向。对垂直于（100）和（001）晶面的冲击，体系在滑移过程中遭遇的剪切阻力较高、持续时间较长，使得能量和温度较快升高，化学反应较容易发生；对垂直于（101）和（111）晶面的冲击，体系在滑移过程中遭遇的阻力较小且出现次数少，使得能量和温度缓慢升高，化学反应不易发生；对其余冲击方向，体系的响应居中。据此评价了7个冲击方向的相对敏感程度：（101）、（111）＜（011）、（110）、（010）＜（100）、（001）。本研究有助于在微观层次深入认识动载荷下TATB的响应机制、结构与性能的关系，为高能低感炸药的设计和研制提供理论参考。

关键词： TATB ; 冲击 ; 滑移 ; 各向异性 ; ReaxFF ; 分子动力学

Abstract

The slip and anisotropy of 2, 4, 6-triamino-1, 3, 5-trinitrobenzene (TATB) crystal under shock loading along various directions were investigated using molecular dynamics simulation combined with reactive force field (ReaxFF). The shock strength was approximately 10 GPa, and seven shock orientations normal to the (101), (111), (011), (110), (010), (100), and (001) crystal planes were considered. For these shock directions, the slip systems that are likely to be activated are predicted to be on the {001} plane, whereas others that could not be activated exhibit large shear stress barriers. These slip characteristics are consistent with the layered structure of TATB crystal along the c axis and the planar structure of TATB molecule. The most favorable slip systems are suggested to be (101)/{001}＜100＞, (111)/{001}＜010＞, (011)/{001}＜010＞, (110)/{001}＜010＞, (010)/{001}＜1$\bar 1$0＞, (100)/{001}＜120＞, and (001)/{001}＜010＞. TATB crystal exhibits anisotropic response to shock loading, that is, the shear stress, energy, temperature, and chemical reactivity during shear deformation depend on shock direction. For the (100) and (001) shock planes, the shear stress barrier is relatively high and lasts for a long time, leading to fast energy accumulation and temperature increment, which, in turn, increase the chemical reactivity. In contrast, for the (101) and (111) shock planes, the small shear stress barrier results in slow energy accumulation and temperature rise and, thus, low chemical reactivity. The (011), (110), and (010) shock planes exhibit intermediate responses. The sensitivity of the seven shock planes can be ranked as follows: (101), (111)＜(011), (110), (010)＜(100), (001). This study provides microscale insight into the response mechanisms and structure-property relationship of TATB crystal under dynamic loading and may facilitate designing explosives with high energy but low sensitivity.

Keywords： TATB ; Shock ; Slip ; Anisotropy ; ReaxFF ; Molecular dynamics

PDF (1840KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

周婷婷, 宋华杰, 黄风雷. 冲击载荷下TATB晶体滑移和各向异性的分子动力学研究. 物理化学学报[J], 2017, 33(5): 949-959 doi:10.3866/PKU.WHXB201702152

ZHOU Ting-Ting, SONG Hua-Jie, HUANG Feng-Lei. The Slip and Anisotropy of TATB Crystal under Shock Loading via Molecular Dynamics Simulation. Acta Physico-Chimica Sinica[J], 2017, 33(5): 949-959 doi:10.3866/PKU.WHXB201702152

1 引言

钝感弹药的研制是当今武器系统发展的重要方向之一，它要求弹药的主装药具有很好的热稳定性和低撞击/冲击感度，以便能避免因意外热力刺激而导致的爆炸甚至爆轰。因此，发展和使用低感高能炸药(insensitive high explosives，IHEs)既能显著提高武器系统的安全性、可靠性和使用寿命，也能避免意外爆轰所导致的人员伤亡和财产损失。

1, 3, 5-三氨基-2, 4, 6-三硝基苯(TATB)是著名的钝感高能炸药，于1887年合成1，从上世纪后半叶作为炸药使用。它是目前为止美国能源部炸药安全委员会批准的唯一的单质钝感炸药，具有极高的热稳定性和抵抗外界刺激的能力，如在冲击、撞击、振动、火烧、光、电磁等意外刺激下不容易发生点火和爆炸2-4。从爆轰性能来看，TATB的威力大于TNT (2, 4, 6-三硝基甲苯)，略小于RDX (三亚甲基三硝胺)。由于其特别优异的安全性能和良好的综合性能，TATB作为最常用的IHE广泛应用于武器、推进剂及民用爆破等领域。

为弄清TATB极低感度的原因，人们对TATB的分子/晶体结构、热力学性质及爆炸性能等进行了大量研究。在微观层次深入认识已知炸药的基本结构，并将基本结构与其已知性质联系起来，将有助于设计制造新的更加安全的炸药并使之得到实际应用。TATB的晶体及分子结构如图1所示。TATB分子(C6H6N6O6)的构型类似于苯环，是平面结构，其六个碳上交替连有三个硝基(―NO2)和三个氨基(―NH2)，具有高度的对称性和稳定性。TATB晶体属三斜晶系，空间群为PĪ，每个单胞中有2个TATB分子，沿c轴方向具有类似于石墨的层状结构，实验测得的晶胞参数为5：a = 0.901 nm，b = 0.903 nm，c = 0.681 nm，α = 108.58°，β = 91.82°，γ = 119.97°，V = 0.443 nm3。在分子内部和层内分子间都存在很强的氢键作用，因此在(001)面形成稳定的二维网络结构，而层与层之间的分子只有很弱的范德华相互作用6。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1   TATB的晶体及分子结构

Fig 1   Crystal and molecular structures of TATB

The different colors represent different atom types: gray represents carbon atom, blue nitrogen atom, red oxygen atom, and white hydrogen atom. color online


已有研究表明，TATB的极低感度和物理化学特性与它的分子及晶体结构紧密相关：TATB分子的高度对称性、稳定的平面结构和分子内的强氢键7-9；类似于石墨的层状晶体结构、层内分子间氢键作用以及分子π堆积形式3, 10-15；摩擦系数小、热传导系数大、比热大，使得热量不易积累16。肖鹤鸣8认为分子内的氢键以及氢键的迁移对TATB的热分解机理、撞击感度和耐热性都有着极其重要的影响。Manaa等9的研究表明，TATB分子内硝基上的氧原子与氨基上的氢原子之间的氢键使得硝基和氨基的二面角扭转存在很高的能垒。Ojeda和Cagin12在原子和电子水平对TATB晶体的研究发现：在静水压缩下分子之间的氢键和分子构型发生显著变化，进而影响到能量、振动光谱和弹性性质；高压下沿着c轴的层间分子也形成了氢键，使体系氢键数量增加，晶体结构更加稳定。Pravica13, 14、Manaa和Fried15等的研究也表明，压力显著地降低了分子间氢键的键长，增强了氢键作用，且高压促进了层间分子间氢键的形成。

TATB晶体中分子的堆积方式使其晶体结构具有明显的各向异性，进而导致了各向异性的热力学性质，如热膨胀17-21、热传导22, 23、熔化24以及静水和单轴压缩下的力学响应12, 15, 21, 25-28。Sun等19采用X射线粉末衍射技术研究了TATB晶体的热膨胀行为，结果表明TATB的体积膨胀主要源于c轴方向的晶格膨胀。Kroonblawd和Sewell22采用非平衡分子动力学(MD)方法研究了常温常压下TATB晶体沿a、b和c轴三个方向的热传导过程，计算得到的热传导系数表明c轴方向的热传导明显低于另外两个方向。Mathew等24采用MD方法模拟了TATB晶体不同晶面的表面熔化过程，结果表明熔化具有明显的各向异性，他们认为这是由于(001)面内分子之间的氢键以及在垂直于(100)和(010)面内堆垛层错的形成趋势所致。Bedrov等21采用MD法研究了TATB晶体在常压200-400 K温度范围内的热膨胀行为和在常温10 GPa内的静水压缩行为，结果表明TATB的弹性系数具有明显的各向异性。第一性原理的计算表明，TATB晶体在单轴压缩下的应力、带隙和能量都具有各向异性28。

凝聚相炸药的起爆被认为与热点形成紧密相关29，由于局部塑性变形而产生的热量是热点产生的机制之一，尤其是对于均质炸药30, 31。Dick等30, 32-36对PETN (季戊四醇四硝酸酯)、HMX (1, 3, 5, 7-四硝基-1, 3, 5, 7-四氮杂环辛烷)等炸药单晶的冲击起爆实验表明：炸药在冲击载荷下的塑性变形与滑移运动直接相关，热力学响应和感度具有明显的各向异性。他们认为30, 32, 33，晶体内滑移面两侧的空间位阻对各向异性的感度起主要作用，敏感晶面的空间位阻大，需要克服较大的能垒才能发生滑移，因此分子变形较大，使得分子被激发到高能量振动态，进而引发化学键的断裂。Sewell课题组37-41采用大规模MD方法从理论上研究了几种常用炸药单晶在冲击载荷下的动力学行为，结果表明炸药的塑性变形以沿着某些滑移系的位错运动为主。Conroy等42, 43认为剪切应力是激发塑性变形的驱动力以及冲击波阵面后化学反应的源泉，因此压缩下的剪切应力对于认识冲击响应的各向异性起关键作用。在此基础上，Zybin等44, 45采用ReaxFF反应力场和MD方法，先后研究了PETN和RDX单晶在压缩-剪切载荷下的物理化学响应，得到的热力学响应和冲击感度的各向异性与冲击实验结果一致。我们46, 47最近采用该方法研究了HMX的冲击各向异性行为，并细致分析了塑性变形与化学反应之间的关联。

Kuklja和Rashkeev48-50采用密度泛函理论(DFT)研究了FOX-7 (1, 1-二氨基-2, 2-二硝基乙烯)和TATB两种晶体在(010)面和(001)面的剪切作用下的化学反应，着重分析了几种初始分解反应的能垒随剪切应变的变化。结果表明：剪切变形使得FOX-7的分解能垒显著降低，但几乎不改变TATB的分解能垒，这可能是FOX-7比TATB更敏感的原因，据此他们认为剪切变形对含能晶体的分子化学键断裂起重要作用。张朝阳51采用DFT方法计算了TATB单胞内单层分子在(001)面滑移需要克服的能量，计算得到的最高能量远低于TATB及其它常用单质炸药分解所需的活化能，因此他认为TATB在(001)晶面的滑移是很可能发生的，且发生在分解反应之前，这可能是其致钝的主要原因。随后，张朝阳等11又比较了TATB、FOX-7及HMX单胞分别在滑移和单轴压缩下分子间作用的变化，结果表明：滑移作用下TATB的分子间作用变化最小，其次是FOX-7，HMX的变化要大得多，说明滑移在TATB晶体中很容易发生；压缩作用导致的内应力变化为TATB＜FOX-7＜HMX，这可能是TATB在机械作用下钝感的原因。但是，该研究中只考虑移动单层分子，且分子结构、层间距离及晶胞体积保持不变。Mathew和Sewell52研究了常压和5 GPa压力下TATB晶体在(001)面的滑移行为，两种压力下计算得到的堆垛层错能都很小，表明在(001)面的滑移很容易发生。他们认为这是由于沿c轴的层间分子间的弱范德华力所致，这将导致较低的滑移运动的临界应力；如果炸药冲击感度的空间位阻模型30, 32, 33对TATB也成立的话，可以推断TATB在(001)面的滑移将非常容易发生，这可能是其极低感度的原因之一。以上对TATB滑移行为的研究基本是人为移动一层分子，使得两层分子之间产生相对位移，然后计算移动前后体系的能量差，以此来判断滑移的可能性，即都是静态或平衡态计算，并不能真实反映动载荷下的非平衡过程。此外，除了(001)晶面，均没有给出在其它晶面的滑移不可能发生的直接证据。因此，目前对TATB在动载荷下的响应机制的研究和认识是不足的。

本文将采用ReaxFF反应力场和分子动力学方法，通过施加压缩-剪切载荷模拟TATB晶体在沿不同方向冲击作用下的动力学过程。与采用MD方法直接模拟冲击相比，该方法对计算资源的要求大幅降低，所需的研究体系更小，能够在很短的时间尺度内(10 ps)预测炸药热力学响应的各向异性，有效区分敏感晶面与不敏感晶面，并能考察大量滑移被激发的可能性44-47。通过本研究，拟获得TATB在冲击载荷下可能被激发的滑移系，各向异性的热力学响应，以及不同冲击方向下体系的相对敏感程度。这将有助于在微观层次深入认识动载荷下TATB的响应机制、结构与其性能的关系，为高能低感炸药的设计和研制提供理论参考。

2 模拟方法及说明

2.1 ReaxFF反应力场

ReaxFF反应力场是基于第一性原理的能够描述材料物理和化学性质的原子间作用势53。其核心是键级，所有与价键相关的势能项都依赖于键级，而键级由原子间距决定，在分子动力学的每一次循环进行计算。ReaxFF能够平滑地描述化学反应过程中化学键连接性的变化，因而能够描述化学键的断裂和生成，为化学反应过程提供丰富的信息。经过十多年的发展，ReaxFF反应力场已经可以正确地描述包括含能材料在内的多种体系在各种复杂条件(如冲击、剪切、高温高压)下的物理化学性质44, 54-57。

为更好地描述分子晶体的结构和密度，在ReaxFF的基础上增加了对伦敦色散力的描述，即修正了分子间长程相互作用力，形成了ReaxFF-lg反应力场58。这对含能材料至关重要，因为炸药的爆轰性能与密度紧密相关。ReaxFF-lg反应力场很好地预测了硝基甲烷(NM)、TATB、RDX及PETN的晶体结构和较高压力范围内的状态方程，对RDX相变的预测也有很大改善58，并且已用于研究HMX晶体在剪切和冲击载荷下的动力学响应46, 47, 59。本文采用的是ReaxFF-lg力场。

2.2 建模及分子动力学模拟细节

(1) ReaxFF-lg力场优化TATB晶体结构。初始晶体结构来源于实验5，首先沿着a、b和c轴将晶体扩大为3 × 3 × 4的超晶胞，含有1728个原子；然后采用共轭梯度法优化体系的原子坐标，得到零温下能量和应力较低的结构；接着采用正则系综(NVT)将体系缓慢升温到300 K，在该温度下对体系的原子坐标进行优化，时间为20 ps，获得300 K下的平衡结构；再采用等温等压系综(NPT)在常温常压(T = 300 K，p = 100000 Pa)下对体系的原子坐标和晶胞参数进行优化，时间为40 ps，获得常温、常压下的平衡结构；对最后10 ps的模拟结果进行统计分析，得到平衡结构的晶胞参数为：a = 0.893 nm，b = 0.891 nm，c = 0.650 nm，V = 0.413 nm3，与实验值5的相对误差分别为0.866%、3.414%、4.536%、6.662%符合得较好，表明ReaxFF-lg力场适用于TATB体系。

(2)确定不同冲击方向下的初始滑移系。选择7个冲击方向，分别垂直于TATB晶体的(101)、(111)、(011)、(110)、(010)、(100)和(001)晶面。基于TATB的晶体结构及分子排布特点，对其滑移系(滑移面和滑移方向)进行初步分析，寻找沿不同方向冲击载荷下的初始滑移系。根据位错理论59，这些滑移系的滑移面两侧的分子间接触应较少，滑移方向的分子排布应较密。

(3)沿冲击方向进行单轴压缩。沿上述七个冲击方向对TATB单胞(NPT获得的平衡结构)进行快速单轴压缩，将压缩后的单胞扩大为3 × 3 × 4的超晶胞，采用NVT系综在常温下(T = 300 K)对应力进行动力学弛豫，弛豫时间为2 ps (体系的能量和应力快速降低并趋于平衡)，统计应力，获得约为10 GPa的压力响应。由于TATB具有明显的各向异性，相同的冲击压力对应各压缩方向的压缩量会有所不同。表1给出了各冲击方向的压缩量和弛豫后体系的应力张量。

表1   各冲击方向的压缩量和NVT-MD优化后压缩体系的应力张量

Table 1  Compression rate for each shock direction and stress tensor for compressed systems after NVT-MD optimization

Shock planeCompression rate/%p/GPapxx/GPapyy/GPapzz/GPapxy/GPapxz/GPapyz/GPa
(101)159.13612.1167.4027.889-0.668-1.5070.748
(111)169.05612.4186.5718.1792.250-3.276-1.555
(011)1510.6819.24811.31311.4822.592-0.310-1.021
(110)1410.72013.1318.04210.9863.225-1.264-0.980
(010)139.8409.4529.77610.2923.192-1.164-1.465
(100)159.55812.1557.3239.197-0.689-0.0961.271
(001)189.9105.5166.85917.3541.3772.1992.137

新窗口打开| 下载CSV


(4)分解切应力计算。将体系的应力张量投影到第二步所确定的初始滑移系，得到每一组滑移系的分解切应力(RSS)。RSS是驱动晶体沿滑移系进行滑移的动力，故RSS越大体系越可能发生滑移。同时考察了冲击晶面与滑移面、冲击方向与滑移方向的夹角，若夹角在45°附近，则该滑移系应具有较大的RSS60。此外，参考TATB的晶体结构特点(沿c轴方向分子按层状排列)以及已有研究结果(在(001)晶面的滑移较容易发生)11, 51, 52，在(001)晶面选取了几组滑移系，虽然它们的RSS比其它滑移系要小得多。表2给出了沿垂直于(001)晶面的冲击方向下具有较大RSS的滑移系和在(001)晶面的滑移系的RSS值。其余冲击方向下的滑移系见表S1 (Supporting Information)。本文将这些滑移系称为潜在滑移系。

表2   沿垂直于(001)晶面的冲击方向下潜在滑移系的RSS值

Table 2  The RSS values for the potential slip systems along the shock direction normal to (001) plane

Shock plane Slip system RSS/GPa
(001) {201}＜10$\bar 2$＞ 6.423
{301}＜10$\bar 3$＞ 6.071
{$\bar 1$01}＜101＞ 5.913
{401}＜10$\bar 4$＞ 5.561
{021}＜01$\bar 2$＞ 5.447
{031}＜01$\bar 3$＞ 5.239
{011}＜01$\bar 1$＞ 5.062
{041}＜01$\bar 4$＞ 4.968
{001}＜120＞ 2.061
{001}＜010＞ 1.826
{001}＜110＞ 1.742

新窗口打开| 下载CSV


(5) MD模拟剪切动力学过程。压缩后的单胞根据第四步确定的潜在滑移系被重新定义，使x轴沿滑移方向，x-z面为滑移面，这主要是为了便于在同一个方向(即x方向)施加剪切载荷。将重新定义后的单胞扩大为较大的超晶胞：以(001)晶面、{001}＜010＞滑移系为例，扩大后的晶胞尺寸为10.696 nm × 9.297 nm × 3.092 nm，超晶胞含1728个分子，41472个原子。采用NVT-MD对超晶胞在常温下进行2 ps的结构优化，再采用微正则分子动力学(NVE-MD)方法模拟体系在恒速率剪切作用下的动力学过程。剪切速率等于剪切速度与沿该剪切方向的滑移距离的比值，大小为0.5 ps-1，剪切面在x-z平面，剪切方向沿x方向。由于TATB很钝感，选用较高的剪切速率是为了能够在较短的模拟时间内快速地分辨出不同冲击方向下体系的差异。NVE-MD的模拟时间为8 ps，时间步长为0.1 fs，采用周期性边界条件。文中MD模拟程序为GRASP (General Reactive Atomistic Simulation Package)61。

2.3 化学反应分析

为了获得剪切动力学过程中的化学反应机制及产物信息，对MD模拟中产生的轨迹文件进行分析。基于每一步MD模拟给出的原子坐标、原子间连接性及键级等信息，获得化学反应生成的产物及其数量随时间的变化；对产物进行编号并跟踪其在动力学过程中的变化，获得产物的反应路径46, 47, 59。由于每组滑移系所对应的模拟体系所含的TATB分子数不一定相等，为了比较不同体系的反应情况，将体系生成的产物数量平均到每个TATB分子。

3 结果与讨论

3.1 各冲击方向下可能被激发的滑移系

表3给出了在垂直于(001)晶面的冲击下体系在剪切过程中的关键物理量。τbarrier为滑移过程中的剪切阻力，即最大剪切应力与初始剪切应力的差值(τbarrier = τmax-τ0)，T为4 ps时体系的温度，treaction表示化学反应开始的时间，NTATB表示模拟结束时体系所剩下的未反应的TATB分子数占初始TATB分子数的百分比。之所以选择测量4 ps时的温度，是因为体系遭遇的剪切阻力集中在前4 ps内，在该过程中滑移面两侧的分子因相互作用会挤压碰撞而变形，体系由晶体结构变为无序结构，温度变化比较显著，因此4 ps时的温度可以较好地体现各体系的差别。从表中可以看出，不同滑移系的响应具有明显差异。对发生概率较高的滑移系，除了具有较大的RSS，还应在滑移过程中遇到较小的剪切阻力，以便更容易克服该阻力而发生滑移。对滑移面不在{001}面的滑移系，体系在滑移过程中遭受到的阻力非常大，有的甚至达到近16 GPa，因此这些滑移系几乎是不可能被激发的。而对滑移面在{001}面的滑移系，体系在滑移过程中遭遇的阻力要小得多，因此比较容易被激发。根据较小剪切阻力原则，在垂直于(001)晶面的冲击下，{001}＜120＞、{001}＜010＞和{001}＜110＞三组滑移系是可能被激发的滑移系。

表3   垂直于(001)晶面的冲击下潜在滑移系的MD模拟结果

Table 3  MD Simulation results for the potential slip systems along the shock direction normal to (001) plane

Shock planeSlip systemτ0/GPaτbarrier/GPaT/Ktreaction/psNTATB/%
(001){201}＜10$\bar 2$＞6.47615.85918600.850
{301}＜10$\bar 3$＞5.96911.38419660.950
{$\bar 1$01}＜101＞5.9597.58014130.903.61
{401}＜10$\bar 4$＞5.5689.72718010.950
{021}＜01$\bar 2$＞5.29210.58816190.800
{031}＜01$\bar 3$＞4.98213.64615910.850
{011}＜01$\bar 1$＞4.7238.16616461.050.26
{041}＜01$\bar 4$＞4.94913.29018210.800
{001}＜120＞2.0652.5717845.4597.45
{001}＜010＞1.9001.6936925.6098.73
{001}＜110＞1.8512.7947274.9097.22

τbarrier is the shear stress barrier, τ0 is the initial shear stress, T is measured at 4 ps, treaction is the time at which chemical reaction begins, NTATB is the undecomposed TATB ratio, measured at 8 ps.

新窗口打开| 下载CSV


在垂直于(001)晶面的冲击下，根据剪切阻力的大小，选取了五组滑移系(表3中黑体文字所示)，分别具有较大({201}＜10$\bar 2$＞滑移系)、居中({01$\bar 1$}＜011＞滑移系)和较小({001}＜120＞、{001}＜010＞和{001}＜110＞三组滑移系)的阻力，比较了它们在剪切过程中的剪切应力、能量、温度和未反应TATB的百分比随时间的变化，如图2所示。对于剪切阻力较大和居中的体系在克服阻力的过程中，体系的能量快速增加，温度急剧升高，化学反应很快发生。相反，剪切阻力较小的体系在克服阻力的过程中，体系的能量和温度都缓慢增加，化学反应进展很慢。在t = 4 ps时，五个体系的温度分别升高至1860、1646、784、692和727 K，剪切阻力小的体系的温升比其它体系低了1000 K左右。在t = 8 ps时，{201}＜10$\bar 2$＞滑移系中已没有未反应的TATB分子，{011}＜01$\bar 1$＞滑移系中只剩下0.26%的TATB分子未反应，而{001}＜120＞、{001}＜010＞和{001}＜110＞三组滑移系中绝大多数的TATB分子都没有发生反应。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2   沿着垂直于(001)晶面的冲击方向下，五组潜在滑移系的剪切阻力、能量、温度和未反应TATB的百分比随时间的变化

Fig 2   Time evolutions of shear stress, potential energy, temperature, and the ratio of undecomposed TATB for the five representative potential slip systems along the shock direction normal to (001) plane

color online


对其余六个冲击晶面，也观察到与(001)晶面类似的现象，在{001}面的滑移系均因较小的剪切阻力而可能被激发，其它滑移系均因较大的剪切阻力而不容易被激发(见表S2(Supporting Information))。即无论是沿着垂直于哪个晶面的冲击，TATB晶体内的滑移都是在{001}面最容易发生，这可以从TATB的晶体内分子的排布特征得到解释(见图1)。TATB晶体内分子沿c轴方向按层状分布，层内分子排列紧密且以氢键相连，而层间分子相距较远且只存在弱分子间相互作用。因此，在(001)面的滑移只需克服弱的分子间作用力而发生，而在其它晶面的滑移都或多或少需要克服层内分子间氢键作用而不容易进行。

为从微观结构更清晰地揭示造成剪切阻力显著差异的原因，图3显示了冲击晶面为(001)，滑移系为{201}＜10$\bar 2$＞、{011}＜01$\bar 1$＞和{001}＜010＞的三个体系的晶体结构。对{201}＜10$\bar 2$＞和{011}＜01$\bar 1$＞滑移系，在滑移过程中都需要克服层内分子间的强相互作用，且会发生分子间的碰撞和挤压，因此剪切阻力大，滑移不容易发生，且容易导致能量和温度的升高。{201}＜10$\bar 2$＞与{011}＜01$\bar 1$＞相比较而言，沿着滑移方向前者的分子间距离更近，更容易发生碰撞和挤压，因此剪切阻力更大，滑移更难进行，能量和温度升高得更快。对{001}＜010＞滑移系，在滑移过程中只需克服层间分子间弱相互作用，且不会发生分子间碰撞和挤压，因此剪切阻力小，滑移容易发生，不易造成能量和温度的升高。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3   冲击晶面为(001)，TATB晶体内不同滑移系在剪切过程中的分子间接触情况

Fig 3   Schematic illustration of molecule contacts for various slip systems in TATB crystal during shear process for (001) shock plane

(a) slip system {201}＜10$\bar 2$＞; (b) slip system {011}＜01$\bar 1$＞, (c) slip system {001}＜010＞. The crystal structures are the ones after redefinition. The different colors represent different atom types: gray represents carbon atom, blue nitrogen atom, red oxygen atom, and white hydrogen atom. The orange arrows indicate the slip systems. color online


表4总结了基于较小剪切阻力原则确定的7个冲击方向下可能被激发的滑移系。这些滑移系往往能量和温度升高较慢，化学反应发生较晚且反应程度较低，这可能是TATB在动载荷下具有较低感度的原因之一。

表4   沿着垂直于(101)、(111)、(011)、(110)、(010)、(100)和(001)晶面的七个冲击方向下可能被激发的滑移系的MD模拟结果

Table 4  MD simulation results for the possible slip systems along the shock directions normal to (101), (111), (011), (110), (010), (100), and (001) planes

Shock planeSlip systemτ0/GPaτbarrier/GPaT/Ktreaction/psNTATB/%Sensitivity
(101) {001}＜120＞0.1790.753(0.813/1.033)459 100
{001}＜110＞0.1540.454(0.908/1.015)482 100
{001}＜100＞0.1190.406(1.020)470 100I
(111){001}＜010＞0.3450.850502 100I
{001}＜120＞0.360.905494899.92
{001}＜110＞0.2940.92(0.797)532899.94
(011){001}＜110＞0.2162.017(1.052/1.171)576899.94
{001}＜120＞0.2341.637(1.617/0.875)540 100
{001}＜010＞0.1340.681(0.538/0.712)552 100M
(110){001}＜100＞0.8350.423(0.769/0.775)6174.198.76
{001}＜010＞0.3090.601(0.235/0.482)5657.299.54M
{001}＜110＞0.3981.797(0.899)6991.8594.51
{001}＜120＞0.1981.544(0.482)6152.3095.19
(010){001}＜010＞0.3101.5656732.2595.70
{001}＜1$\bar 2$0＞0.3291.704623698.96
{001}＜1$\bar 1$0＞0.3561.4245987.699.65M
(100){001}＜110＞0.3842.3666844.696.39
{001}＜120＞0.3071.5236556.7599.17S
(001){001}＜120＞2.0652.5717845.4597.45
{001}＜010＞1.91.693(0.384)6925.6098.73S
{001}＜110＞1.8512.7947274.9097.22

sensitivity: (I) insensitive, (M) intermediate, (S) sensitive. The shear stresses shown in the parenthesis are the second and third shear stress barriers.

新窗口打开| 下载CSV


3.2 各向异性的冲击响应

为比较TATB晶体在不同冲击方向下的响应，根据最小剪切阻力原则，确定了各冲击方向下最容易被激发的滑移系(表4中加粗文字所示)。图4比较了这7组滑移系在剪切过程中的应力、能量、温度及化学反应随时间的变化。对垂直于(101)晶面的冲击，体系先克服了第一个很小的剪切阻力(0.406 GPa)，随后又克服了第二个较大的剪切阻力(1.02 GPa)，总的持续时间约为1 ps，之后不再遭遇明显的阻力。对垂直于(111)晶面的冲击，体系在1.5 ps内克服了0.85 GPa的剪切阻力，之后不再遭遇明显的阻力。因此，在垂直于(101)和(111)晶面的冲击作用下，体系在初始滑移阶段所遇到的阻力较小，且该阻力在短时间内就能被克服，之后滑移很容易进行。对垂直于(011)和(110)晶面的冲击，体系在滑移过程中所遇到的阻力虽然较小(0.6 GPa左右)，但次数较多，因此其滑移要比(101)和(111)晶面更困难一些。对垂直于(100)和(010)晶面的冲击，体系都需要在初始滑移阶段克服一个较大的阻力(分别为1.523和1.424 GPa)，且(100)晶面阻力所持续的时间较长约为2 ps，说明两者的滑移都较难进行，且(100)晶面遇到的阻力更难克服。对垂直于(001)晶面的冲击，体系在初始滑移阶段先后克服了一个较大的阻力(1.693 GPa)和一个较小的阻力(0.384 GPa)，之后又陆续遭遇多次大小不一的阻力，说明在该冲击方向下，体系在滑移过程中遇到的阻力较多且较难克服。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4   沿着垂直于(101)、(111)、(011)、(110)、(010)、(100)和(001)晶面的七个冲击方向下7组最易被激发的滑移系在剪切过程中的应力、能量、温度以及化学反应随时间的变化

Fig 4   Time evolutions of shear stress, potential energy, temperature, and the ratio of undecomposed TATB for the seven most favorable slip systems along the shock directions normal to (101), (111), (011), (110), (010), (100), and (001) planes during shear process

color online


不同冲击方向下体系剪切应力的差别导致体系的能量变化有所不同，如图4(b)所示。由于体系在滑移过程中需要克服分子间作用力以及滑移面两侧的分子发生接触、挤压导致分子变形，使得体系的势能发生变化。在克服剪切阻力的过程中，体系的势能升高较快，阻力被克服后，势能的变化更加缓慢。对(101)和(111)晶面，由于体系在滑移过程中遭遇的阻力小且次数少，势能的升高比较缓慢；(111)晶面相比(101)晶面由于阻力更大且持续时间更长，势能的初始上升速率更大。对(011)和(110)晶面，体系在滑移过程中遭遇的阻力虽小但次数较多，因此势能的升高比(101)和(111)晶面稍快。对(100)和(010)晶面，体系在克服阻力的过程中势能升高得较快，势能曲线在早期出现明显的局域峰值，这与较大的初始剪切阻力相符；与(010)晶面相比，(100)晶面因更长的阻力持续时间使得势能保持了更长时间的快速升高。对(001)晶面，因体系在滑移过程中克服的阻力最大且次数较多，势能的升高也最快，势能曲线在早期出现明显的局域峰值。

体系温度的变化与势能对应，在时间上稍微有所滞后，如图4(c)所示。能量升高得越快的体系，温度也越高。对垂直于(101)、(111)、(011)、(110)、(010)、(100)和(001)晶面的冲击，体系的温度在t = 4 ps时分别升高到470、502、552、565、598、655和692 K。温度是控制化学反应的重要因素，不同冲击方向下体系的温升不同必然会导致化学反应的差异。图4(d)比较了不同冲击方向下体系在滑移过程中未反应的TATB分子数占初始TATB分子数的百分比随时间的变化。对(101)、(111)和(011)晶面，体系在模拟时间内没有发生反应；对(110)、(010)、(100)和(001)晶面，虽然都只有少数TATB分子发生了反应，但(100)和(001)晶面的化学反应发生得更早，消耗的反应物更多。

综上，在不同方向的冲击作用下，TATB单晶在滑移过程中的应力、能量、温度及化学反应都表现出差异，即各向异性的冲击响应。(100)和(001)晶面因相对较高的剪切阻力和较长的持续时间使得体系出现较快的能量积聚和温度升高，导致较为快速的化学反应，因此比较敏感；(101)和(111)晶面因较小的剪切阻力和较少的阻力次数使得体系的能量和温度升高缓慢，化学反应不易发生，因此较为钝感；其余晶面居中。因此，7个冲击晶面的相对敏感程度为：(101)、(111)＜(011)、(110)、(010)＜(100)、(001)，(001)晶面最敏感。

为进一步分析化学反应的差别，选取了三组敏感程度不同的冲击晶面(101)、(110)和(001)进行比较。对反应产物的分析发现，对钝感的(101)晶面，体系没有发生化学反应，无产物生成。对敏感度居中的(110)晶面，主要发生了C―H键断裂生成游离态H原子和TATB分子间反应生成少量C12H12O12N12大分子。对敏感的(001)晶面，反应生成了较多C12H12O12N12和少量C18H18O18N18大分子，以及少数HONO、HNO、NH2和NH3等小分子。因此，对敏感晶面的冲击，会诱发更多的反应通道，生成更多的产物。

此外，为考察体系的敏感程度对剪切速率的依赖性，研究了上述三个体系在剪切速率为1.0和1.5 ps-1时的动力学响应。结果表明，体系的敏感程度与剪切速率无关，仍然是(101)＜(110)＜(001)。

4 结论

本文采用ReaxFF反应力场和分子动力学方法，通过施加压缩-剪切载荷模拟了TATB晶体在沿不同方向冲击作用下的动力学过程。冲击方向分别垂直于(101)、(111)、(011)、(110)、(010)、(100)和(001)晶面，冲击强度为10 GPa。根据分解切应力和剪切阻力的大小，预测了各冲击方向下可能被激发的滑移系。这些滑移系都在{001}晶面，剪切阻力小，其它滑移系均因很大的阻力而不容易被激发，这与TATB晶体沿c轴的层状结构和TATB分子的平面结构相符。对在{001}面的滑移系，体系在剪切过程中的能量和温度升高较慢，化学反应发生较晚且反应程度较低，这可能是TATB在动载荷下具有较低感度的原因之一。基于最小阻力原则，给出了七个冲击方向下最可能被激发的滑移系，分别为(101)/{001}＜100＞、(111)/{001}＜010＞、(011)/{001}＜010＞、(110)/{001}＜010＞、(010)/{001}＜1$\bar 1$0＞、(100)/{001}＜120＞和(001)/{001}＜010＞。TATB晶体在剪切过程中的应力、能量、温度和化学反应与冲击方向相关，即具有各向异性的冲击响应。对(100)和(001)冲击晶面，滑移过程中体系遭遇的剪切阻力相对较高、持续时间较长，使得能量和温度升高较快，化学反应较早发生，因此体系比较敏感。对(101)和(111)冲击晶面、滑移过程中遭遇的剪切阻力较小，且出现次数少，使得能量和温度缓慢升高，化学反应不易发生，因此体系较为钝感；对其余冲击晶面，体系的响应处于中等程度，敏感程度居中。因此，7个冲击晶面的相对敏感程度为：(101)、(111)＜(011)、(110)、(010)＜(100)、(001)，(001)晶面最敏感。本研究为动载荷下TATB晶体的滑移和各向异性提供了微观水平的认识，有助于揭示炸药晶体结构与其性能的关系，为高能低感炸药的设计和研制提供理论参考。

Supporting Information:  The RSS values for the potential slip systems along the shock directions normal to (101), (111), (011), (110), (010), and (100) planes are summarized in Table S1. MD simulation results for the potential slip systems along the shock directions normal to (101), (111), (011), (110), (010), and (100) planes are summarized in Table S2. This information is available free of charge via the internet at http://www.whxb.pku.edu.cn.

参考文献

Jackson C. L. ; Wing J. F. J. Am. Chem. Soc. 1887, 9, 354.

[本文引用: 1]

Cady H. H. ; Larson A. C. Acta Cryst. 1965, 18, 485.

[本文引用: 3]

Agrawal J. P. Progr. Energy Combust. Sci. 1998, 24, 1.

[本文引用: 1]

Wu C. J. ; Fried L. E. J. Phys. Chem. A 2000, 104, 6447.

[本文引用: 1]

Xiao H. M. The Theory of the Molecular Orbits for Nitrocompound Beijing: National Defence Industry Press, 1993.

[本文引用: 1]

肖鹤鸣. 硝基化合物的分子轨道理论, 北京: 国防工业出版社, 1993.

[本文引用: 1]

Manaa M. R. ; Gee R. H. ; Fried L. E. J. Phys. Chem. A 2002, 106, 8806.

[本文引用: 2]

Roszak S. ; Gee R. H. ; Balasubramanian K. ; Fried L. E. Chem. Phys. Lett. 2003, 374, 286.

[本文引用: 1]

Zhang C. Y. ; Wang X. ; Huang H. J. Am. Chem. Soc. 2008, 130, 8359.

[本文引用: 2]

Ojeda O. U. ; Cagin T. J. Phys. Chem. B 2011, 115, 12085.

[本文引用: 2]

Pravica M. ; Yulga B. ; Liu Z. X. ; Tschauner O. Phys. Rev. B 2007, 76, 064102.

[本文引用: 1]

Pravica M. ; Yulga B. ; Tkachev S. ; Liu Z. X. J. Phys. Chem. A 2009, 113, 9133.

[本文引用: 1]

Manaa M. R. ; Fried L. E. J. Phys. Chem. C 2012, 116, 2116.

[本文引用: 3]

Dong H. S. Chin. J. Energy Mater. 2004, 12, 1.

[本文引用: 1]

Kolb J. R. ; Rizzo H. F. Propellants Explos. Pyrotech. 1979, 4, 10.

[本文引用: 1]

Gee R. H. ; Roszak S. ; Balasubramanian K. ; Fried L. E. J. Chem. Phys. 2004, 120, 7059.

Sun J. ; Kang B. ; Xue C. ; Liu Y. ; Xia Y. X. ; Liu X. F. ; Zhang W. Journal of Energetic Materials 2010, 28, 189.

[本文引用: 1]

Taylor D. E. J. Phys. Chem. A 2013, 117, 3507.

Bedrov D. ; Borodin O. ; Smith G. D. ; Sewell T. D. ; Dattelbaum D. M. ; Stevens L. L. J. Chem. Phys. 2009, 131, 224703.

[本文引用: 3]

Kroonblawd M. P. ; Sewell T. D. J. Chem. Phys. 2013, 139, 074503.

[本文引用: 2]

Kroonblawd M. P. ; Sewell T. D. J. Chem. Phys. 2014, 141, 184501.

[本文引用: 1]

Mathew N. ; Sewell T. D. ; Thompson D. L. J. Chem. Phys. 2015, 143, 094706.

[本文引用: 2]

Liu H. ; Zhao J. ; Du J. ; Gong Z. ; Ji G. F. ; Wei D. Q. Phys. Lett. A 2007, 367, 383.

[本文引用: 1]

Stevens L. L. ; Velisavljevic N. ; Hooks D. E. ; Dattelbaum D. M. ; Propellants Explos. Pyrotech. 2008, 33, 286.

Valenzano L. ; Slough W. J. ; Perger W. AIP Conf. Proc. 2012, 1426, 1191.

Budzevich M. M. ; Landerville A. C. ; Conroy M. W. ; Lin Y. ; Oleynik I. I. ; White C. T. J. Appl. Phys. 2010, 107, 113524.

[本文引用: 2]

Bowden F. P. ; Yoffe A. D. Initiation and Growth of Explosion in Liquids and Solids 1st ed Cambridge: Cambridge University Press, 1985.

[本文引用: 1]

Dick J. J. ; Mulford R. N. ; Spencer W. J. ; Pettit D. R. ; Garcia E. ; Shaw D. C. J. Appl. Phys. 1991, 70, 3572.

[本文引用: 4]

Armstrong R.W. ; Ammon H. L. ; Elban W. L. ; Tsai D. H. Thermochim. Acta 2002, 384, 303.

[本文引用: 1]

Dick J. J. ; Ritchie J. P. J. Appl. Phys. 1994, 76, 2726.

[本文引用: 3]

Dick J. J. J. Appl. Phys. 1997, 81, 601.

[本文引用: 2]

Yoo C. S. ; Holmes N. C. ; Souers P. C. ; Wu C. J. ; Ree F. H. ; Dick J. J. J. Appl. Phys. 2000, 88, 70.

Dick J. J. ; Hooks D. E. ; Menikoff R. ; Martinez A. R. J. Appl. Phys. 2004, 96, 374.

Menikoff R. ; Dick J. J. ; Hooks D. E. J. Appl. Phys. 2005, 97, 023529.

[本文引用: 1]

Jaramillo E. ; Sewell T. D. ; Strachan A. Phys. Rev. B 2007, 76, 064112.

[本文引用: 1]

Ramos K. J. ; Hooks D. E. ; Sewell T. D. ; Cawkwell M. J. J. Appl. Phys. 2010, 108, 066105.

Cawkwell M. J. ; Ramos K. J. ; Hooks D. E. ; Sewell T. D. J. Appl. Phys. 2010, 107, 063512.

Bedrov D. ; Hooper J. B. ; Smith G. D. ; Sewell T. D. J. Chem. Phys. 2009, 131, 034712.

Eason R. M. ; Sewell T. D. J. Phys. Chem. C 2012, 116, 2226.

[本文引用: 1]

Conroy M. W. ; Oleynik I. I. ; Zybin S. V. ; White C. T. Phys. Rev. B 2008, 77, 094107.

[本文引用: 1]

Conroy M. W. ; Oleynik I. I. ; Zybin S. V. ; White C. T. J. Appl. Phys. 2008, 104, 053506.

[本文引用: 1]

Zybin S. V. ; GoddardI W. A. II. ; Xu P. ; van Duin A. C. T. ; Thompson A. P. Appl. Phys. Lett. 2010, 96, 081918.

[本文引用: 3]

An Q. ; Liu Y. ; Zybin S. V. ; Kim H. ; Goddard III W. A. J. Phys. Chem. C 2012, 116, 10198.

[本文引用: 1]

Zhou T. T. ; Zybin S. V. ; Liu Y. ; Huang F. L. ; Goddard W. A. III. J. Appl. Phys. 2012, 111, 124904.

[本文引用: 3]

Song H. J. ; Zhou T. T. ; Huang F. L. ; Hong T. Acta Phys.-Chim. Sin. 2014, 30, 2024.

URL     [本文引用: 4]

宋华杰; 周婷婷; 黄风雷; 洪滔. 物理化学学报, 2014, 30, 2024.

URL     [本文引用: 4]

Kuklja M. M. ; Rashkeev S. N. Appl. Phys. Lett. 2007, 90, 151913.

[本文引用: 1]

Kuklja M. M. ; Rashkeev S. N. J. Phys. Chem. Lett. 2010, 1, 363.

Kuklja M. M. ; Rashkeev S. N. Journal of Energetic Materials 2010, 28, 66.

[本文引用: 1]

Zhang C. Y. J. Phys. Chem. B 2007, 111, 14295.

[本文引用: 2]

Mathew N. ; Sewell T. D. Philosophical Magazine 2015, 95, 424.

[本文引用: 2]

van Duin A. C. T. ; Dasgupta S. ; Lorant F. ; Goddard III W. A. J. Phys. Chem. A 2001, 105, 9396.

[本文引用: 1]

Zhou T. T. ; Shi Y. D. ; Huang F. L. Acta Phys.-Chim. Sin. 2012, 28, 2605.

URL     [本文引用: 1]

周婷婷; 石一丁; 黄风雷. 物理化学学报, 2012, 28, 2605.

URL     [本文引用: 1]

Strachan A. ; van Duin A. C. T. ; Dasgupta S. ; Chakraborty D. ; Goddard III W. A. Phys. Rev. Lett. 2003, 91, 098301.

Nomura K. ; Kalia R. K. ; Nakano A. ; Vashishta P. Appl. Phys. Lett. 2007, 91, 183109.

An Q. ; Zybin S. V. ; Goddard III W. A. ; Botero A. J. ; Blanco M. ; Luo S. N. Phys. Rev. B 2011, 84, 220101.

[本文引用: 1]

Liu L. C. ; Liu Y. ; Zybin S. V. ; Goddard III W. A. J. Phys. Chem. A 2011, 115, 11016.

[本文引用: 2]

Zhou T. T. ; Lou J. F. ; Zhang Y. G. ; Song H. J. ; Huang F. L. Phys. Chem. Chem. Phys. 2016, 18, 17627.

[本文引用: 3]

Wang Y. N. ; Chen S. J. ; Dong X. C. Dislocation Theory and Its Application Beijing: Metallurgical Industry Press, 2007.

[本文引用: 1]

王亚男; 陈树江; 董希淳. 位错理论及其应用, 北京: 冶金工业出版社, 2007.

[本文引用: 1]

[本文引用: 1]

/

〈 〉