1 引言
高聚物黏结炸药(PBXs)主要由高分子聚合物黏结剂与炸药颗粒间实现包覆粘连形成，其具有许多特殊性能，由于黏结剂的加入，对混合体系的感度及力学性能均有所改善[1]。1,1’-二羟基-5,5’-联四唑二羟胺盐(TKX-50)作为一种多氮元素的含能离子盐，具有爆速高，机械感度与冲击波感度较低的特性[2]，苯并三氧化呋咱(BTF)具有较高的能量与较好的传爆性能[3-5]。因此，将固相高能炸药TKX-50与BTF混合以期综合两者的优良性能。
一直以来，构建合适的多组分炸药体系以解决单一炸药在应用中存在的缺陷被众多学者广泛研究。张蕾等[6]根据第一性原理对CL-20/BTF共晶结构进行计算，结果表明相对无氢BTF晶体，CL-20/BTF共晶的热力学稳定性得到了较大提升。王伟等[7]运用最小自由能法计算了含TKX-50/CL-20的聚叠氨缩水甘油醚(GAP)基固体推进剂能量性能，预估了TKX-50/CL-20复配推进剂的应用效能。刘佳辉等[8]制备了TKX-50/F2314混合炸药，与常见的HMX、RDX、TATB基混合炸药相比其具有较高的爆速，撞击感度显著降低，且对冲击波刺激较为钝感。Yu Yuehai等[9]通过分子动力学(MD)模拟计算了4种黏结剂F2311、F2641、聚乙二醇(PEG)和乙烯-醋酸乙烯酯共聚物(EVA)与TKX-50基炸药的力学性能，表明PEG对TKX-50力学性能改善效果最佳。Ma song等[10]计算出TKX-50/三氟氯乙烯(PCTFE)的结合能一般低于TKX-50/PVDF，且黏结剂的加入对TKX-50的力学性能得到改善，感度有所降低。
本研究针对TKX-50/BTF炸药体系，搭建了F2311、F2602与直链式酚醛树脂(PF)3种高聚物黏结剂结构，进行MD模拟计算炸药体系与不同黏结剂体系间结合能、界面作用方式及力学性能，为TKX-50/BTF双组分炸药筛选物理相容性较好、力学性能改善效果较优的热塑性黏结剂提供理论支撑，并对改善混合炸药体系压制成型性有一定的实际应用价值。
2 计算模型与计算方法
2.1 高聚物黏结剂分子结构搭建
本文选用3种热塑性高分子黏结剂：F2311、F2602与PF，搭建相应PBXs模型。F2311为偏氟乙烯(PVDF)与PCTFE以 1∶1摩尔比形成的无规共聚物，在本次模拟中端基加H或F以饱和，取链接数n=10；F2602与F2311同属二元系氟橡胶，为PVDF与六氟丙烯(HFP)形成的聚合物，尾端加H以饱和；PF是以酚类与醛类在酸或碱的催化作用下相互缩合而成的高分子化合物，其因催化条件的不同分为热塑性与热固性两类[11-12]，本文搭建的热塑性酚醛树脂是一种线性预聚物，其分子结构为线性结构，端基以苯酚加以饱和。后续计算均在Materials Studio 2020中完成。
利用Build Polymer工具搭建黏结剂分子模型，经结构优化后在Forcite模块中执行MD计算，在NVT系综下，温度设定为298 K，计算时间为2 500 ps，待体系能量收敛后，取计算结果最后一帧为黏结剂稳定结构，3种黏结剂结构如图1所示。
图1 3种高聚物黏结剂分子模型
Fig.1 Three molecular models of polymer binders
2.2 炸药分子表面静电势计算
模拟搭建的两种炸药TKX-50[13]与BTF[14]晶体结构均取自单晶衍射数据，TKX-50属单斜晶系，空间群为P21/c，BTF属斜方晶系，空间群为Pna21。TKX-50作为一种含能离子盐[15]，在模拟过程中为保证计算的可靠性，需利用Dmol3模块对TKX-50及BTF的电荷进行计算分配。密度泛函选择GGA-PBE，基组选择DNP，自洽收敛标准(SCF tolerance)定为Fine，最大自洽循环次数(Max.SCF cycles)定为200步，同时进行电子密度(electron density)、表面静电势(electrostatics)、布居分析(population analysis)的计算，计算将重新分配电荷给每一个原子，结果以炸药分子的电子等密度轮廓线形式表达。
如图2所示，根据密度泛函理论(DFT)，分子的局部正静电势反映键的强弱，局部正静电势越大，其相应键的键能反而越低[16-17]，TKX-50的强正电区域集中于羟胺阳离子(NH3OH+)处，BTF具有很高吸电子成分的氮杂原子，在分子的中心部分消耗电子电荷，因此在此区域留下较强的正电势，强正电区域集中于其缺电性六元环处。在后续计算中，均以Mulliken电荷分布作为原子电荷，TKX-50的羟胺阳离子局部正电势较高，在MD模拟过程中易游离分散于体系模型中
图2 TKX-50与BTF分子表面的静电势
Fig.2 Electrostatic potential on the surface of TKX-50 and BTF
2.3 PBXs构型搭建及MD计算细节
为保证PBXs模型与实际情况相符，实现黏结剂与炸药层的充分接触，采用Build Layers工具搭建PBXs层状模型，即TKX-50层/黏结剂/BTF层构型。使用Morphology模块的BFDH方法分别对TKX-50及BTF的主要晶面进行预测，发现TKX-50的主要晶面为(0 1 1)、(0 2 0)、(1 0 0)、(1 1 0)及(1 1 -1)，其所占面积分别为42.47%、19.88%、18.42%、15.68%、3.55%，BTF的主要晶面为(0 2 0)、(1 1 0)、(0 1 1)、(0 1 -1)、(1 1 1)及(1 1 -1)，其所占面积分别为34.98%、32.96%、15.36%、15.36%、0.67%、0.67%。因此，本文选取TKX-50中的(0 1 1)晶面与BTF的(0 2 0)晶面搭建PBXs层状模型进行后续研究。
为比较不同黏结剂与炸药体系的界面作用，此次模拟选择搭建的炸药体系TKX-50与BTF质量比均为7∶3，黏结剂体系相对总混合体系的质量占比为5%左右。首先，对TKX-50单胞进行5×3×4的超晶胞扩建，BTF单胞扩建为4×3×4的超晶胞，分别切割TKX-50超晶胞的(0 1 1)晶面与BTF的(0 2 0)晶面，使用Build Layers工具，根据黏结剂质量占比，将5条F2311链放入炸药层间，完成TKX-50/BTF/F2311体系PBXs初始构型的搭建。对F2602及PF体系的PBXs同样进行扩建超晶胞、切割晶面、构建层状结构的操作，将3种PBXs分别记作TKX-50/BTF/F2311、TKX-50/BTF/F2602和TKX-50/BTF/PF，相应的黏结剂含量分别为 5.4%、5.1%、4.7%。
选用PCFF力场，使用Forcite模块对3种体系进行结构优化。对优化后结构进行MD计算时，选用NPT系综，采用Nose恒温器和Berendsen恒压器，在常压下，温度设置为298 K，时间步长1 fs，总模拟计算时长400 ps，其中前200 ps用于体系结构的平衡，后200 ps用于数据统计分析。计算完成后得到的TKX-50/BTF/F2311、TKX-50/BTF/F2602和TKX-50/BTF/PF平衡构型如图3所示。通过对体系平衡结构进行单点能的计算分析，得到不同黏结剂与炸药体系的结合能，利用径向分布函数进行体系界面作用进行非键力的作用分析，最后任务选择Mechanical Properties，方式选择恒定应变，对3种PBXs体系进行力学性能的预测。
图3 3种PBXs体系平衡构型
Fig.3 Equilibrium structures of the three kinds of PBXs models
3 结果与讨论
3.1 体系平衡的判定
在MD模拟中，针对不同类型的体系，平衡的判断方式也不相同。本次模拟的PBXs体系包含炸药层与黏结剂层，在NPT系综下进行控温与控压，因此需要体系温度、能量、密度均达到收敛条件，才能判定体系达到平衡。
以TKX-50/BTF/PF体系为例，其温度、能量及密度曲线如图4所示，从密度变化的角度分析，体系初始构型的搭建需确保黏结剂与炸药分子均处于模型盒子中，因此初始结构较为稀疏，密度偏低，MD模拟过程中在设定压力条件下对盒子进行压缩，体系密度逐渐上升，最终体系温度、能量及密度的瞬时值均随模拟进行在可控范围内存在小幅波动，其平均值已趋于稳定，证明此时体系已经达到平衡。
图4 TKX-50/BTF/PF体系平衡收敛曲线
Fig.4 TKX-50/BTF/PF convergence curve
3.2 黏结剂与炸药体系的结合能
对各体系进行炸药与黏结剂间结合能(Ebind)的计算，结合能定义为分子间相互作用能的负值。针对PBXs体系，结合能反映的是黏结剂与炸药体系混合的能力，结合能越大，表明炸药体系与黏结剂间的相互作用越强，进而反映体系物理相容性越好[18]。炸药体系与黏结剂的结合能可由式(1)计算得到：
Ebind=-(Etotal-Epolymer-Ebase)
(1)
式中：Etotal为PBXs平衡构型下的平均总能量；Epolymer为体系去掉炸药部分黏结剂的单点能；Ebase为在平衡构型下TKX-50与BTF去掉黏结剂链体系的单点能；
定义为单位黏结剂(质量分数%)计算结合能的归一值，结合能计算结果如表1所示。
表1 不同黏结剂与TKX-50/BTF炸药体系的结合能(kJ·mol-1)
Table 1 Binding energies of different binders and the TKX-50/BTF explosive system
BinderEtotalEbaseEpolymerEbindE′bindF2311-63465.9-53542.1-2234.17689.71424.0F2602-28432.4-16208.0-5484.06740.41321.6PF-43277.1-41013.84370.66633.91411.5
结合能具有容量性质，相互比较时需换算成单位质量才有意义。由表1数据可知，3种不同黏结剂与炸药体系的结合能大小顺序为F2311>PF>F2602。这3种热塑性高聚物黏结剂在PBXs中作为惰性组分表现较为稳定，故结合能的大小反映着体系的物理相容性，结合能越大，表明黏结剂与炸药体系的相互作用越强，形成的PBXs更为稳定。因此，F2311与炸药体系TKX-50/BTF的相容性优于F2602与PF。
3.3 径向分布函数
通过对径向分布函数(RDF)的计算分析可判断高聚物黏结剂与炸药层间界面相互作用，通常分子间的相互作用包括氢键作用与范德华力(vdW)，在RDF分析中，氢键的距离一般在2.6～3.1 Å，强范德华作用通常在3.1～5.0 Å[19]，在g(r)曲线中上述距离范围内有较明显峰的出现，可以说明存在相对稳定的相互吸引作用。因此，选取与高聚物黏结剂距离较近的炸药层进行分析，将炸药层分子中的氢、氧、氮分别标记为H、O、N，F2311中氟与氢分别标记为F1、H1，F2602中氟与氢分别标记为F2、H2，PF中氧与氢分别标记为O1、H3，3种PBXs的g(r)曲线如图5所示。
图5 3种PBXs的径向分布函数
Fig.5 Radial distribution function of the three kinds of PBXs
由图5(a)F1-H曲线可知，F2311中的氟(F1)与炸药层中的H同时出现在2.6 Å附近的概率较大，其g(r)值最大达到6.5左右，表示这两种原子间存在较强的氢键作用，且在3.2 Å附近再次出现峰值，表明该原子对间还存在一定的vdW作用，图5(a)O-H1与N-H1曲线可见，炸药层中的O原子和N原子与F2311中的H原子(H1)主要存在较弱的氢键作用。因此，TKX-50/BTF与F2311间存在氢键作用与vdW作用，以氟氢原子对产生的氢键作用为主，这也是他们之间相互作用较强的主要根源。由图5(b)F2-H曲线可见，g(r)曲线在2.6 Å与3.4 Å附近均出现较高峰值，表明在F2602中的氟(F2)与炸药层中的H之间同时存在一定的氢键作用与vdW作用，图5(b)N-H2曲线在3.6 Å附近出现峰值约为1.4，说明炸药层中N原子与F2602中H原子(H2)存在较弱的vdW作用，图5(b)O-H2曲线峰值出现在2.1?附近，表示炸药层中O原子与F2602中的H原子(H2)也存在一定的氢键作用。图5(c)O1-H曲线表明PF与炸药层中O1-H原子对也存在一定的氢键和vdW作用，图5(c)N-H3曲线呈现逐渐上升趋势，说明PF中H(H3)与炸药层中N存在着vdW作用。
3.4 PBXs的力学性能
优异的力学性能对压装型PBXs炸药显得尤为重要，本文分别对3种加入了不同黏结剂的PBXs层状模型进行力学性能计算。将TKX-50/BTF/F2311、TKX-50/BTF/F2602和TKX-50/BTF/PF分别命名为PBXⅠ、PBXⅡ和PBXⅢ，计算得到各体系的弹性系数矩阵Cij下表2中。
表2 TKX-50/BTF混合体系与3种PBXs的弹性系数(GPa)
Table 2 Elastic constants of three PBXs and TKX-50/BTF
ElasticConstantTKX-50/BTFPBXⅠPBXⅡPBXⅢC1123.6119.3117.5719.52C2224.6318.4721.0319.22C3322.6817.8019.3019.99C446.014.824.754.07C555.504.464.694.40C665.334.214.004.37C1212.1310.149.8310.61C1312.3410.009.3410.90C2312.988.979.9710.24
依据Voigt理论与Reuss理论方法[20]，分别给出晶体的体积模量KV、KR，剪切模量GV、GR的计算公式
9KV=(C11+C22+C33)+2(C12+C23+C13)
(2)
15GV=(C11+C22+C33)-(C12+C23+C13)+
3(C44+C55+C66)
(3)
1/KR=(S11+S22+S33)+2(S12+S23+S13)
(4)
15/GR=4(S11+S22+S33)-4(S12+S23+S13)+
3(S44+S55+S66)
(5)
将2种理论方法计算的模量进行算术平均得到平均体积模量K与平均剪切模量G，并根据K、G、杨氏模量E、泊松比ν之间的关系得到杨氏模量E与泊松比ν的计算公式：
(6)
(7)
杨氏模量E与材料抵抗弹性形变的能力有关，其值越小，表明材料的弹性越好；剪切模量G越小，表明材料的塑性越好；体积模量K用于关联材料的断裂强度，其值越大，表明材料的断裂强度越大；K/G为材料体积模量与剪切模量的比值，用来表示材料受到较强冲击作用而不致破坏的性能，其值越大材料的韧性越强；柯西压(C12-C44)可用来衡量材料的延展性，其值为负数，则材料显脆性，若为正，则表明材料的延展性较好[21-22]。据此计算出各体系的力学性能，如表3所示。
表3 TKX-50/BTF混合体系与3种PBXs的力学性能
Table 3 Mechanical properties of three kinds of PBXs and TKX-50/BTF
ElasticConstantTKX-50/BTFPBXⅠPBXⅡPBXⅢTensilemodulus(E/GPa)7.536.276.205.96Poisson’sratio(ν)0.350.330.340.35Bulkmodulus(K/GPa)8.116.316.466.81Shearmodulus(G/GPa)2.802.352.302.20K/G2.892.742.813.10C12-C446.125.325.086.54
由表3数据分析可知，相对TKX-50/BTF混合体系，加入了黏结剂的3种PBXs的E、K与G值均有所下降，表明黏结剂的加入使材料的刚性减弱，柔性增强；4种体系的泊松比均在0.3左右，与可压材料实际情况较为符合，具有一定的塑性；各体系的C12-C44值均为正数，表明TKX-50/BTF混合体系自身具有一定的延展性，而随着黏结剂PF的加入，TKX-50/BTF/PF体系的C12-C44值与K/G值相较TKX-50/BTF有所增加且数值最大，而F2311与F2602的加入对体系的力学性能改善效果并不明显。因此预测随着PF的加入，材料的韧性与延展性有所提升，力学性能的改善效果优于其他两种氟橡胶F2311与F2602。
4 结论
1) 以结合能判断体系物理相容性来分析，预示与TKX-50/BTF混合体系的相容性顺序：F2311>PF>F2602，其中F2311与混合体系结合能最大，表明其相较F2602与PF更适合加入TKX-50/BTF混合体系中，形成的PBXs稳定性更好。
2) 径向分布函数g(r)分析表明，黏结剂F2311、F2602和PF与炸药层间存在一定强度的氢键作用与vdW作用，其中以F2311中的F与炸药层中的H存在较强的氢键作用，同时也从侧面说明了其结合能较大、相容性较好的原因。
3) TKX-50/BTF与3种添加了不同黏结剂的PBXs力学性能计算结果表明，黏结剂的加入均使体系刚性减弱，柔性增强，力学性能得以改善；其中PF对体系力学性能的改善效果最为突出，使材料韧性与延展性得到增强，更利于PBXs的加工、运输与使用。
[1] Niu H,Chen S S,Jin S H,et al.Preparation,nonisothermal decomposition kinetics,heat capacity,and safety parameters of TKX-50-based PBX[J].Journal of Thermal Analysis And Calorimetry,2018,131(03):3193-3199.
[2] 张坤,王健,王晓峰,等.1,1′-二羟基-5,5′-联四唑二羟胺盐在火炸药中的应用研究进展[J].火炸药学报,2022,45(01):20-29.
Zhang K,Wang J,Wang X F,et al.Research progress on the application of 1,1′-dihydroxy-5,5′-bitetrazole dihydroxyamine salt in explosives[J].Journal of Explosives,2022,45(01):20-29.
[3] 陈捷,汪佩兰,张孝仪,等.微量杂质对苯并三氧化呋咱 (BTF) 感度影响研究[J].化学研究应用,2004(01):120-122.
Chen J,Wang P L,Zhang X Y,et al.Study on the effect of trace impurities on the sensitivity of benzotrioxyfurazan (BTF)[J].Application in Chemical Research,2004(01):120-122.
[4] Ju X H,Xiao H M,Chen L T.Periodic DFT Approach to Benzotrifuroxan Crystal[J].International Journal of Quantum Chemistry,2005,102(07):224-229.
[5] Zhang H B,Guo C Y,Wang X C,et al.Five Energetic Cocrystals of BTF by Intermolecular Hydrogen Bond and π-Stacking Interactions[J].Crystal Growth &Design,2013(13):679-687.
[6] 张蕾,赵艳红,姜胜利,等.CL-20及其共晶炸药热力学稳定性与爆轰性能的理论研究[J].含能材料,2018,26(06):464-470.
Zhang L,Zhao Y H,Jiang S L,et al.Theoretical study on thermodynamic stability and detonation performance of CL-20 and its eutectic explosives[J].Energetic Materials,2018,26(06):464-470.
[7] 王伟,王健,付晓梦,等.TKX-50/CL-20复配在固体推进剂中应用效能的预估[J].含能材料,2021,29(09):827-832.
Wang W,Wang J,Fu X M,et al.Prediction of the application performance of TKX-50/CL-20 compound in solid propellant[J].Energetic Materials,2021,29(09):827-832.
[8] 刘佳辉,范桂娟,卢校军,等.TKX-50基混合炸药的爆轰及安全性能[J].含能材料,2019,27(11):902-907.
Liu J H,Fan G J,Lu X J,et al.Detonation and safety performance of TKX-50-based mixed explosives[J].Energetic Materials,2019,27(11):902-907.
[9] Yu Y,Chen S,Li X,et al.Molecular Dynamics Simulations for 5,5′-bistetrazole-1,1′-diolate (TKX-50) and its PBXs[J].RSC Adv,2016(01):112-118.
[10] Ma S,Li Y J,Li Y,et al.Research on structures,mechanical properties,and mechanical responses of TKX-50 and TKX-50 based PBX with molecular dynamics[J].Mol Model,2016,22:43-53.
[11] 马子强.酚醛树脂的合成及应用研究进展[J].合成树脂及塑料,2021,38(06):71-74.
Ma Z Q.Research progress on the synthesis and application of phenolic resins[J].Synthetic Resins and Plastics,2021,38(06):71-74.
[12] 孙鑫,王刚,凌芳.酚醛树脂的研究进展及发展趋势[J].上海化工,2022,47(01):59-63.
Sun X,Wang G,Ling F.Research progress and development trend of phenolic resin[J].Shanghai Chemical Industry,2022,47(01):59-63.
[13] Fischer N,Fischer D,Klap?tke T M,et al.Pushing the limits of energetic materials-the synthesis and characterization of di-hydroxylammonium 5,5’-bistetrazole-1,1’-diolate[J].Journal of Materials Chemistry,2012,22(38):20418-20422.
[14] Golovina N I,Titkov A N,Raevskii A V,et al.Kinetics and Mechanism of Phase Transitions in the Crystals of 2,4,6-Trinitrotoluene and Benzotrifuroxane[J].Journal of Solid State Chemistry,113(02)：229-238.
[15] 刘英哲,毕福强,来蔚鹏,等.5,5’-联四唑-1,1’-二氧二羟铵在不同生长条件下的晶体形貌预测[J].含能材料,2018,26(03):210-217.
Liu Y Z,Bi F Q,Lai W P,et al.Crystal morphology prediction of Dihydroxylammonium 5,5’-Bistetrazole-1,1’-diolate under different growth conditions[J].Energetic Materials,2018,26(03):210-217.
[16] Murray J S,Lane P,Politzer P.Relationships between impact sensitivities and molecular surface electrostatic potentials of nitroaromatic and nitroheterocyclic molecules[J].Molecular Physics,1995,85(01)：1-8.
[17] Murray J S,Lane P,Politzer P.Effects of strongly electron-attracting components on molecular surface electrostatic potentials:application to predicting impact sensitivities of energetic molecules[J].Molecular Physics,1998,93(02)：187-194.
[18] 许晓娟,肖继军,黄辉,等.ε-CL-20基PBX结构和性能的分子动力学模拟——HEDM理论配方设计初探[J].中国科学:B辑,2007,37(06):556-563.
Xu X J,Xiao J J,Huang H,et al.Molecular dynamics simulation of the structure and properties of ε-CL-20-based PBX:A preliminary study on the theoretical formulation design of HEDM[J].Science in China:Series B,2007,37(06):556-563.
[19] 肖继军,朱卫华,朱伟,等.高能材料分子动力学[M].北京:科学出版社,2013.
Xiao J J,Zhu W H,Zhu W,et al.Molecular dynamics of high energy materials[M].Beijing:Science Press,2013.
[20] HILL R.The elastic behaviour of a crystalline aggregate[J].Proceedings of the Physical Society.Section A,1952,65(05):349-354.
[21] 杨犁,余庚泽,余晨,等.基于TKX-50的PBXs含能材料力学性能计算模拟[J].实验技术与管理,2021,38(03):61-68.
Yang L,Yu G Z,Yu C,et al.Computational simulation of mechanical properties of PBXs energetic materials based on TKX-50[J].Experiment Technology and Management,2021,38(03):61-68.
[22] 孙婷,肖继军,赵锋,等.CL-20/DNB共晶基PBXs相容性、界面作用和力学性能的MD模拟[J].含能材料,2015,23(04):309-314.
Sun T,Xiao J J,Zhao F,et al.MD simulation of compatibility,interfacial interaction and mechanical properties of CL-20/DNB eutectic-based PBXs[J].Energetic Materials,2015,23(04):309-314.