ε-CL-20在二元混合溶剂中晶体形貌的分子动力学模拟

赵盐鹏,卫宏远,党乐平

(天津大学化工学院,天津 300072)

自20世纪Nielsen首次合成以来,六硝基六氮杂异伍兹烷(hexanitrohexaazaisowurtzitane,简称CL-20)在航天和军事等领域一直备受关注[1]。CL-20作为一种新型含能材料,相较于奥克托今(HMX)和黑索金(RDX)等含能材料,CL-20具有更高的爆速,更大的爆压,更高的能量及密度。CL-20具有α、β、γ和ε 4种晶型,相较于其他晶型,ε-CL-20的能量密度最高,最具有实际应用价值[2]。然而,过高的感度使得ε-CL-20在使用过程中存在极大的安全隐患,容易引发爆炸等危险[3]。因此,如何降低ε-CL-20的感度成为当前含能材料领域的主要研究方向。

目前,降低ε-CL-20晶体感度的方法主要是通过与其他含能材料晶体形成共晶、制备聚合物黏性炸药(PBXs)以及改进工艺[4-7]。通过形成共晶、黏性炸药虽然可以在一定程度上降低晶体感度,但也会降低其能量密度。溶剂化物形成技术等改进的生产工艺由于受到实际生产条件的限制,往往无法在含能材料的实际生产中得到运用[8]。晶体性质由外部环境和晶体的内部结构共同决定。通过改变溶剂、添加剂或操作条件可以达到控制晶体形貌的目的。实验结果表明,晶体形貌对感度有着十分重要的影响。相较于椭球形晶体,针状或片状晶体具有更高的摩擦感度和撞击感度[9]。因此,改善晶体形貌,使其更趋近于椭球形是降低晶体感度的一种十分有效的方案。

目前,关于晶体球形化生产工艺的改进已经有了很多相关报道。朱明河等[10]通过探究结晶条件及晶体生长环境得到了制备球形氯化钠晶体的改进工艺;任晓婷等[11]采用溶剂-反溶剂法和超声波辅助技术重结晶,在乙酸乙酯-正庚烷体系中制备了细粒度椭球形ε-CL-20,显著降低了晶体的摩擦感度和撞击感度;Guo等[12]采用球磨法在乙醇-水的二元溶剂体系下得到了超细椭球形ε-CL-20晶体,晶体的粒度、热解温度、摩擦感度和撞击感度都得到了降低;荆肖凡等[13]以乙酸乙酯为溶剂,正庚烷为非溶剂,采用自制的喷雾辅助结晶装置制备了椭球形ε-CL-20,显著降低了ε-CL-20晶体的冲击感度、摩擦感度与热爆炸临界温度;徐洋等[14]通过采用喷雾和超声辅助重结晶装置制备了超细球形化ε-CL-20颗粒,与原料相比,超细ε-CL-20的感度明显降低,热爆炸临界温度略微提高。

分子动力学模拟是近年来飞速发展的一种分子模拟方法,以牛顿经典力学、量子力学、统计力学为基础,利用计算机数值求解分子体系经典力学运动方程的方法得到体系的相轨迹,并统计体系的结构特征与性质[15]。Song等[16]通过分子动力学方法,采用修正附着能(MAE)模型预测了3,4-二硝基-1H-吡唑(DNP)在不同溶剂影响下的晶体形貌,并通过实验结果验证,结果表明预测结果与实验结果一致。Lan等[17]通过分子动力学方法预测了不同温度下1,1-二氨基-2,2-二硝基乙烯(FOX-7)在二甲基亚砜(DMSO)溶剂影响下的晶体形貌,并与实验结果进行了比较,二者高度吻合。陈巍等[18]通过分子动力学方法预测了对二苯酚晶体在水溶液中的晶体形貌,探究了溶剂以及杂质对晶体溶解度、介稳区及晶习的影响。虽然通过分子动力学模拟方法探究晶体形貌已经在含能材料领域得到了广泛应用,但截至目前,分子动力学模拟方法对ε-CL-20的晶习研究还有待发掘。

本课题组通过先前的研究探究了不同二元溶剂体系对ε-CL-20晶体形貌的影响,得出了在乙酸乙酯/二溴甲烷混合溶剂中晶体形貌最接近椭球形。采用Material Studio 7.0软件,通过分子动力学方法,从分子层面上探究二元溶剂体系与ε-CL-20各个晶面之间的相互作用,对于加深对晶体生长机理的理解和晶体形貌的调控具有积极的意义,同时,为含能材料的形貌设计提供一个高效的方法。

1 计算理论与细节

1.1 计算理论

1980年,P Hartman和P Bennema提出了周期性键链理论(PBC),基于这一理论提出了Attachment Energy(AE)模型,能够预测晶体在真空条件下的形貌[19]。相较于Bravais-Friedel-Donnay-Harker(BFDH)模型,AE模型不仅考虑了晶体内部结构,也考虑了外部环境对晶体形貌的影响,因此对晶体形貌的预测更为准确。附着能Eatt是指一个生长单元附着在晶体表面上释放的能量,定义式为[20]:

式(1)中:Elatt是指晶体的能量,Eslice是指厚度为dhkl的生长单元的能量。真空下晶体各晶面的相对生长速率Rhkl,假定与附着能Eatt的绝对值呈正比:

然而,在实际条件下,很多晶体都是在溶剂环境中生长的。在溶剂环境中,溶剂分子会吸附在晶体表面,晶体若想生长则必须有生长单元吸附才行,为此必须要使溶剂分子从晶体表面脱离,而这个过程会消耗能量,也会降低晶体各晶面的生长速率。为了预测溶剂环境下生长晶体的形貌,需要对AE模型进行修正,引入能量修正项Es,其定义式如式(3):

式(3)中:Aacc为晶体晶面的溶剂可及面积,可通过Materials Studio软件的Atoms Volumes & Surfaces模块进行计算;Abox为双层模型(h k l)方向平面的的面积;Eint为溶剂分子与晶体晶面间的相互作用能,其计算如式(4):

式(4)中:Etot是双层模型中溶剂层和晶面层的总能量;Esur是晶面层的能量;Esol是溶剂层的能量。

1.2 模拟细节

CL-20晶体在自然界共有4种晶型,分别为α、β、γ和ε。其中,ε晶型最为稳定,能量最高,也最具有实际应用价值[3]。根据赵信岐等[21]解构的单晶数据,构建了ε-CL-20晶胞。其晶格参数为a=1.3696 nm,b=1.2554 nm,c=0.8833 nm,α=90°,β=111.18°,γ=90°。

在构建好ε-CL-20晶胞后,采用MS 软件中的Forcite模块进行能量最小化,计算精度为Fine,截断距离为1.55 nm,力场选用COMPASS力场。用Morphology模块中的Growth Morphology算法预测ε-CL-20晶体在真空下的形貌,如图1所示,得到真空中晶体的重要晶面,分别为(1 1 0)、(0 0 1)、(1 1 -1)、(0 1 1)、(2 0 0)和(2 0 -1),各晶面相关信息列于表1。

图1 AE模型得到的真空条件下的l(a)晶体形貌;(b)分子间相互作用Fig.1 (a) Crystal morphology and (b) molecular interactions of ε-CL-20 in vacuum predicted by AE model

表1 优化晶胞参数和实验晶胞参数Table 1 Optimized lattice parameters and experimental lattice parameters

依据Lan等[22]的研究,双层模型中溶剂层与晶体层的厚度不应低于截断距离,长和宽不应低于2倍的截断距离。将ε-CL-20晶胞以3层厚度按照重要晶面的方向进行切割,构建超晶胞,使其参数不低于3.1 nm。将乙酸乙酯分子与二溴甲烷按照体积比4∶1～1∶4构建溶剂层,使溶剂层各边参数不低于3.1 nm。构建双层模型,将溶剂层置于晶体层上方,2层之间留有0.5 nm的间隙,溶剂层上方留有5.0 nm真空层,以消除周期性结构附加自由边界的影响。

对构造好的双层模型进行能量最小化,约束晶体层第1层的晶体分子,然后进行1 000 ps分子动力学模拟。选用COMPASS力场,NVT系综,Anderson控温法控制温度为298 K。静电力由Ewald法计算,计算精度为4.186×10-4kJ·mol-1,范德华力由Atom based法计算,截断距离为1.55 nm,模拟步长为1.0 fs,每500步输出一帧结果。待体系达到平衡后,取后400 ps做动力学分析。用自己编写的Perl Script脚本计算平衡态下溶剂层与晶体层的相互作用能。

2 实验部分

2.1 实验材料

ε-CL-20,质量分数为99.5%,北京理工大学提供;二溴甲烷,纯度AR,上海迈瑞尔化学技术有限公司。

2.2 ε-CL-20晶体的制备

通过溶析结晶法制备ε-CL-20晶体。将过量ε-CL-20原料溶解于乙酸乙酯溶剂中,静置12 h,取一定量上清液加入100 mL三颈烧瓶中,水浴加热控制温度在298 K,以0.2 mL·min-1的速率缓慢加入一定量二溴甲烷溶剂,溶剂与反溶剂体积比为4∶1～1∶4,搅拌速率为150 r·min-1,溶剂滴加完成后停止搅拌,静置12 h。得到单晶后将其收集并进行进一步表征。实验装置如图2所示。

图2 实验装置设备图Fig.2 Flow chart of experiment device

2.3 测试条件与表征

1)扫描电子显微镜:采用德国蔡司EVO MA 15。测样前,用导电碳带将试样固定在试板上,在真空中用高速电子轰击试样,采集信号进行形貌分析。

2)X-射线粉末衍射:采用X-射线粉末衍射仪(D/max-2500,Rigaku,Japan)收集CL-20 溶剂化物的PXRD 谱图。测试条件如下:射线源Cu_Kα,波长λ=0.1541 nm,电压40 kV,电流40 mA,扫描角度(2θ)范围5°～40°,扫描速率5(°)·min-1,扫描步长0.02°,测试温度为室温。

3 结果与讨论

3.1 主要晶面相互作用分析

采用AE模型预测了ε-CL-20晶体在真空下的形貌,如图1(a)所示。在真空条件下,ε-CL-20形貌近似椭球形,长径比为1.436,有6个主要晶面,分别为(0 0 1)、(0 1 1)、(1 1 0)、(1 1 -1)、(2 0 0)和(2 0 -1),各晶面的具体信息列于表2。图1(b)展示了单位晶胞内分子间的相互作用,其中颜色越接近红色表示相互作用越弱,颜色越接近蓝色表示相互作用越强。根据周期性键链理论,相互作用越强的方向晶面生长速度越快,Zhao等[23]的研究证明了这一点。生长速度越快的晶面在晶体最终形貌中占比越小,甚至容易消失,生长速率越慢的晶面在最终形貌中越重要。在ε-CL-20晶体的6个晶面中,(1 1 0) 晶面占比39.62%,最为重要,这意味着该晶面在真空条件下生长速度最慢。

表2 真空条件下ε-CL-20的晶习参数Table 2 Crystal habit values of ε-CL-20 in vacuum conditions

由(3)式得知,晶体修正附着能与溶剂分子的可接触面积有关。溶剂分子的可接触面积由不同晶面晶体的堆叠结构相关。为了进一步探究溶剂分子与晶体分子之间相互作用的机理,各晶面的堆叠结构如图3所示。在所有晶面中,(0 1 1)晶面最不平整,Gang等[24]的研究表明,晶面越不平整,越有利于在溶剂分子与晶体分子形成强非键相互作用,这有利于溶剂分子吸附于该晶面上。晶面的平整程度可以通过粗糙度S表示,S的定义式如式(7):

图3 ε-CL-20主要晶面上的Connolly面Fig.3 The geometric structures of different ε-CL-20 planes represented by the Connolly surface model

式(7)中:Aacc为单位晶胞中晶面的溶剂可及面积,Ahkl为单位晶胞中晶面的表面积,S的值越大,则表明该晶面粗糙度越大,越不平整,越有利于在溶剂分子与晶体分子形成强非键相互作用,生长速率越慢,在晶体的最终形貌中越重要。

表3 ε-CL-20晶习表面的参数Table 3 The parameter values for the crystal habit planes of ε-CL-20

此外,晶体表面的极性也会影响溶剂分子与晶体分子间的相互作用。晶体晶面暴露的基团决定了该晶面的极性,极性越大,与溶剂分子的相互作用越强,溶剂分子越容易吸附在晶面上。为了系统性的探究晶面极性与溶剂分子极性对二者相互作用的影响,采用DMol3模块PBE的GGA梯度修正泛函计算了各晶面晶体分子的表面静电势。

如图4所示,对于各晶面的硝基基团,静电势为负值,主要为蓝色区域,这表明硝基具有强负极性,而笼状结构中的氢原子静电势为正值,主要为红色区域,表明笼状结构中的氢原子具有强正极性。正负电荷区交替出现在周期结构中,更有利于溶剂分子附着在平面上。由于(0 1 1)面具有较强的静电势点和较大的粗糙度,(0 1 1)面上的相互作用能最强。

图4 ε-CL-20主要晶面的静电势Fig.4 Electrostatic potentials of ε-CL-20 crystal planes

3.2 溶剂层与晶体层相互作用分析

溶剂-晶体相互作用能可以描述溶剂分子与晶体之间的相互作用。溶液中溶剂与晶体之间存在相互作用,阻碍了溶质分子的附着,抑制了晶面的生长。每个晶面的生长速率发生变化,导致晶体最终形貌发生变化。为了研究反溶剂分子对相互作用能的影响,表4和图5显示了不同二元溶剂体系中的相互作用能。

表4 ε-CL-20在不同V(乙酸乙酯)∶V(二溴甲烷)混合溶剂中的相互作用能和晶习aTable 4 Interaction energy and crystal habit of ε-CL-20 in ethyl acetate/dibromomethane mixed solvents with different volume ratiosa

不同体系的相互作用能是不同的。由图5可知,二溴甲烷做反溶剂时,随着V(乙酸乙酯)∶V(二溴甲烷)发生变化,溶剂层与晶体层的相互作用也随之变化。二溴甲烷体积分数越大,溶剂层与晶体层的相互作用越弱,反之则越强。这表明溶剂层与晶体层间的相互作用主要由乙酸乙酯与CL-20提供,而二溴甲烷与CL-20晶体间的相互作用较弱,即二溴甲烷作反溶剂减弱了乙酸乙酯与CL-20晶体的相互作用。

图5 不同V(乙酸乙酯)∶V(二溴甲烷)下混合溶剂中溶剂层与晶体层相互作用能Fig.5 Interaction energy between solvent layer and crystal layer with different volume ratio of ethyl acetate to dibromomethane

3.3 径向分布函数(RDF)分析

径向分布函数(RDF)可用于分析溶剂层与晶体层之间相互作用的组成,通过量化指定原子为中心的周围原子的分布情况来探究原子间相互作用力的性质。通常,氢键作用力的范围为0～0.31 nm,范德华力(vdW)作用范围0.31～0.50 nm,库伦力作用范围大于0.5 nm。为了探究晶体层与溶剂层之间相互作用能的组成,对双层模型进行RDF分析。以溶剂层V(二溴甲烷)∶V(乙酸乙酯)为1∶1时的双层模型为例,图6展示了(0 0 1)晶面的溶剂分子与晶体分子之间的径向分布函数。

图6 乙酸乙酯与二溴甲烷体积比为1∶1时CL-20(0 0 1)晶面的RDFFig.6 RDF of (0 0 1) plane of CL-20 while the volume ratio of ethyl acetate to dibromomethane is 1∶1

由图6可以看出,在(0 0 1)晶面上,在0～0.31和0.50 nm以外范围内均可以发现尖峰,这表明溶剂层与晶体层之间的相互作用主要由氢键作用力与库仑力组成,而范德华力相较于氢键作用力与库仑力较小,在晶体层与溶剂层相互作用的组成中贡献较小。晶体层与溶剂层之间的氢键作用力和库仑力主要由晶体层的氢原子与溶剂层的氧原子贡献,这表明影响晶体层与溶剂层分子相互作用的主要是CL-20笼状结构上的氢原子和乙酸乙酯溶剂中的氧原子。

3.4 均方位移(MSD)分析

溶剂分子在各晶面的扩散系数表示了该分子在各晶面的扩散能力,这对晶体各晶面的生长有着至关重要的作用。在该研究中,通过爱因斯坦的扩散方程计算了乙酸乙酯溶剂分子在各晶面的扩散系数D,计算方程如式(8)所示。

式(8)中:ri表示第i个粒子的位置向量,角括号表示系综平均值。因此,可以通过均方位移随时间变化的斜率计算三维随机布朗运动粒子的扩散系数D。采用Forcite模块计算了298 K,1 atm下各体系中乙酸乙酯分子在(0 1 1)晶面上0～300 ps的MSD,每0.5 ps输出一帧结果,如图7所示。

图7 不同V(乙酸乙酯)∶V(二溴甲烷)体系下(0 1 1)晶面的MSDFig.7 MSD of (0 1 1) plane in the system of ethyl acetate and dibromomethane with different volume ratio

由图7可知,在(0 1 1)晶面上,V(乙酸乙酯)∶V(二溴甲烷)为4∶1,2∶1,1∶1,1∶2,1∶4的体系中MSD分别为1.219×10-8m2·s-1,1.148×10-8m2·s-1,1.108×10-8m2·s-1,0.766×10-8m2·s-1,0.739×10-8m2·s-1,混合溶剂中乙酸乙酯的体积分数越大,扩散系数越高,乙酸乙酯更容易扩散到晶面上。由于溶剂层与晶体层的相互作用主要由乙酸乙酯中的氧原子与晶体层中的氢原子提供,因此扩散系数越高,乙酸乙酯在二元溶剂体系中的扩散速率越大,同时溶剂层与晶体层的相互作用也越强。因此,当混合溶剂中V(乙酸乙酯)∶V(二溴甲烷)为4∶1时,乙酸乙酯的扩散速率最高,此时溶剂层与晶体层的相互作用也最强,结果与图5一致。

3.5 粉末X射线衍射(PXRD)分析

图8显示了分别在不同体系中的CL-20晶体与ε-CL-20晶体的PXRD。结果表明,实验得到的晶体仍为ε晶型,晶型未发生变化。

图8 不同体系下ε-CL-20的PXRDFig.8 PXRD of ε-CL-20 under different systems

3.6 溶剂效应对晶体形貌的影响

表4中的模拟结果表明,在乙酸乙酯与二溴甲烷不同体积比混合溶剂中生长的ε-CL-20晶体具有不同的长径比和主要晶面。在所有体系中,由于dhkl的面间距较大,(1 1 0)和(0 0 1)晶面总是存在的。同时,由于(2 0 0)晶面面粗糙度低,极性弱,不易被溶剂分子吸附,使得溶剂分子对(2 0 0)面生长的抑制作用最弱,任何体系中都不可能存在(2 0 0)面。从图3可以看出,(0 1 1)面和(2 0 -1)面上暴露了大量具有强负极性的硝基,有利于溶剂分子的吸附,形成强相互作用。

基于以上讨论,溶剂效应和晶体内部结构是影响晶体形貌的重要因素。在溶剂极性、官能团、分子结构、扩散速率等因素的影响下,不同的溶剂分子可能以氢键、范德华力和库仑力的形式选择性地与不同粗糙度的主要晶面相互作用。二元溶剂体系中晶体形貌的模拟和实验结果如图9所示。

图9列出了不同V(乙酸乙酯)∶V(二溴甲烷)混合溶剂中CL-20的模拟形貌与SEM得到的实验结果。当V(乙酸乙酯)∶V(二溴甲烷)为1∶2时,模拟得到晶体形貌的长径比为1.817,且最终形貌中存在的主要晶面最多,此时CL-20晶体形貌最接近椭球形。在5种组成的二元溶剂中,预测结果与实际晶习拟合度较高,个别晶面的比例与实际晶体有所偏差,这是由于分子模拟结果处于理想的平衡状态。在实际实验条件下,受限于实验时间、方法、温度、搅拌速度等实验条件的影响,很难达到理想的平衡状态,这是造成模拟结果与实验结果差异的主要原因,但总的来说预测出晶体形貌的趋势可以为晶体设计提供重要指导。

图9 CL-20在不同体系下的模拟形貌与实验结果Fig.9 Simulation crystal morphology and experimental results of CL-20 in different systems

综上所述,通过研究不同溶剂分子在平面上的相互作用,分子动力学模拟是预测ε-CL-20在不同溶剂体系中的形貌和晶体表面积变化趋势的一种可行方法。

4 结论

通过溶析结晶法制备了不同体系下的ε-CL-20晶体。通过分子动力学模拟的方法对不同体系下的晶体形貌进行了模拟,模拟结果与实验结果拟合良好。此外,通过RDF和MSD分析进一步探究了相互作用能以及溶剂分子在不同表面的扩散系数,并比较分析了不同反溶剂分子影响下相互作用能、氢键作用力以及扩散系数的变化规律,该研究的主要结论总结如下。

1)随着V(乙酸乙酯)∶V(二溴甲烷)的变化,溶剂层与晶体层相互作用及溶剂对晶面的抑制作用发生改变,影响了各晶面的生长速率,使得ε-CL-20的晶体形貌也随之改变。当混合溶剂中V(乙酸乙酯)∶V(二溴甲烷)为1∶2时,长径比为1.817,ε-CL-20晶体的形貌更接近椭球形。

2)RDF分析表明溶剂层与晶体层间的相互作用主要为ε-CL-20(H)与乙酸乙酯(O)提供的氢键相互作用与库仑力相互作用。

3)MSD分析表明,V(乙酸乙酯)∶V(二溴甲烷)对溶剂的扩散速率有显著影响。混合溶剂中乙酸乙酯的体积分数越大,扩散系数越高,溶剂层与晶体层间的相互作用越强。

4)PXRD分析表明,乙酸乙酯/二溴甲烷混合溶剂体系下得到的CL-20晶体仍为ε晶型,这表明该体系下通过溶剂-反溶剂重结晶法得到的晶体爆轰性能未受到影响。

本研究对CL-20晶体的结晶具有一定的指导意义,有助于控制合适的晶体生长环境得到想要的晶体。