0 引言
聚叠氮缩水甘油醚(GAP)是一种侧链上含有叠氮基团的含能聚合物,具有密度大(1.3 g/cm3)、能量高、热稳定性好等优点,是高能推进剂一种常见的含能黏合剂,但由于GAP分子侧链上存在大体积侧基(—CH2N3),阻碍了分子链的运动,降低了主链的柔韧性,导致其力学性能尤其是低温力学性能较差[1],因此,如何提高GAP体系的力学性能是近来的研究热点之一。如王鑫等[2]研究了GAP与多炔基化合物(TPTM)固化反应活性和胶片性能,结果表明,当TPTM质量分数由3%增加到9.7%时,TPTM胶片拉伸强度由0.16 MPa增加到0.82 MPa,断裂伸长率由149%降到17%;胡义文等[3]以三羟甲基丙烷(TMP)为交联剂,制备得到了低交联度GAP聚氨酯弹性体并分析其力学性能,结果表明,经交联后的GAP弹性体常温力学性能显著提高;赵昱等[4]研究了不同固化剂对GAP胶片力学性能的影响,发现多异氰酸酯(N-100)/二官能度的甲苯二异氰酸酯(TDI)双固化体系能够提高GAP黏合剂胶片的强度和延伸率。Agawane等[5]研究了异氰酸酯固化GAP黏合剂的改性与应用,发现N-100的引入对GAP力学性能有显著的增强作用。综上所述,N-100与GAP固化能够显著改善GAP的低温力学性能。
由于推进剂在制备、运输和使用的过程中不可避免的会受到外界刺激作用,严重时会导致事故的发生。同时温度的极端条件和反应的复杂性限制了通过实验对爆炸过程的深入了解。近年来,分子模拟技术被广泛应用于诸如含能材料等复杂反应体系的研究中。如Lu等[6]采用分子动力学模拟的方法研究了改性GAP的玻璃化转变温度以及力学性能。Ma等[7]通过耗散力子动力学和分子模拟的方法研究了GAP/HTPB混合物的相平衡。Wu等[8]则利用分子动力学模拟方法探究GAP的增塑机理。反应力场能够与分子动力学(MD)方法相结合[9-10],可以探究含能材料的热分解过程。为了深入研究GAP/N-100交联体系的性能,本文中使用自编的Perl脚本,通过Materials Studio(MS)软件建立了GAP/N-100交联模型,分析不同交联度下GAP/N-100体系的力学性能等。基于ReaxFF-lg反应力场,使用LAMMPS(large scale atomic/molecular massively parallel simulator)模拟软件,针对高交联度的GAP/N-100体系,研究了其在不同温度下的热分解及初始分解机理、反应速率、活化能以及主要产物演化,为探究GAP推进剂的性能和热解机理提供理论指导。
1 建模与计算
1.1 分子模型构建细节
GAP和N-100的分子结构如图1所示。使用MS软件包中的Amorphous Cell模块建立的GAP/N-100无定型结构如图2所示,其中,包含30条重复单元为10的GAP分子链与10个N-100组成的边长为36.37 Å的10个立方体盒子模型,模型的初始密度为1.3 g/cm3,每个模型的原子数为4 370。
图1 GAP与N-100单体结构
Fig.1 Structure of GAP and N-100 monomers
图2 GAP/N-100无定型结构
Fig.2 Amorphous cell of GAP/N-100
为了使模型结构更加合理,对这10个模型采用Smart Minimize方法进行结构优化后,选择能量最低的模型进行了500 ps的温度为600 K、压力为1个大气压的NPT系综下的MD模拟,为了方便后续交联操作,将异氰酸酯基上的C原子标记为R1,GAP分子链段的羟基上的O标记为R2。
1.2 交联模型建模细节
GAP分子中的羟基与N-100中的异氰酸酯基反应机理如图3所示。模拟过程中主要根据R1与R2的距离判断交联反应能否发生,即当R1与R2的距离小于设定的反应半径时交联反应发生,R1与R2生成新的化学键,否则交联反应不会发生。本文参考文献[11]利用Perl语言在MS中编写了可以实现GAP与N-100自动交联的脚本。
图3 羟基与异氰酸酯基反应机理
Fig.3 Hydroxyl and isocyanate group reaction mechanism
为了能够更快的发生交联反应,节省模拟时间,将模拟反应温度设定为600 K,设定初始的反应半径为3.5 Å,当此范围内没有可以反应的基团时,逐步增大反应半径(每次增加0.5 Å)直到达到最大的反应半径(12.5 Å)。为了使交联结构更加合理每次生成新交联键后都要进行能量最小化和100 ps的NPT系综的MD模拟。
交联完成后首先对交联结构进行300 K升温到600 K,然后再降温到300 K的退火处理,温度间隔为50 K,每个温度进行120 ps的NPT MD模拟,然后又进行了298 K,500 ps的NVT MD的模拟和1.01×105 Pa条件下500 ps的NPT MD模拟,最后500 ps NPT MD模拟的结果用于性能分析。MD模拟过程中采用Nose方法进行控温、Berendsen方法进行控压,范德华和库伦力分别采用Atom Based和Ewald方法,时间步长为1 fs,力场选择COMPASSⅡ力场,周期性边界条件。如图4为最终交联度达到96.7%体系的模型。
图4 交联度为96.7%的体系
Fig.4 System with a crosslinking conversion of 96.7%
1.3 热解机理模拟参数
使用LAMMPS软件研究了在ReaxFF-lg反应力场条件下交联度为96.7%的GAP/N-100交联体系在高温下的热分解情况。具体模拟细节如下:首先对MS得到的交联模型在反应力场条件下进行结构弛豫,为了防止弛豫时出现化学键断裂重排的情况,在弛豫时添加了键约束,通过不断减小bond_coeff的能量项,使得键约束最终逐渐消除。然后在4 000、3 500、3 000、2 500和2 000 K下进行100 ps的NVT系综下的反应分子动力学(RMD)模拟,控温方式选择Nose-Hoover方法,时间步长0.1 fs,键级为0.3,用于判定是否成键。
2 结果与分析
2.1 力学性能
选择不同交联度的模型采用静态恒应变法进行力学性能分析,通过交联MD模拟,最终得到了不同交联度的GAP/N-100交联分子模型。其过程如下:假设系统初始为平衡状态,此时应用一个微小的应变到该系统,重新进行能量优化。将该应变施加到不同方向,重复多次该过程,可以计算刚性矩阵为势能相对于应变的二阶导数,如式(1):
(1)
式(1)中:应力为单位体积势能对应变的一阶导数;“+”代表拉伸;“-”代表压缩。得到刚性矩阵后进一步可以根据式(2)以及式(3)计算Lame常数λ和μ。
(2)
(3)
计算出Lame常数后由式(4)—式(7)得到各项力学性质为
(4)
G=μ
(5)
(6)
(7)
其中,E、G、K、ν分别为杨氏模量、剪切模量、体积模量和泊松比。
通过以上公式,计算了5个交联度下的E、G、K、ν,结果如表1所示。从表1中可以看出,随着交联度的增大,E、G、K大致呈现增大的趋势,在交联度96.7%时达到最大值。
表1 不同交联度时体系的力学性能
Table 1 Mechanical properties of the system at different crosslinking conversion
交联度E/GPaG/GPaK/GPaν80.0%3.429 21.367 02.326 10.254 383.3%3.579 81.389 02.822 70.288 686.7%3.634 81.391 73.121 00.305 993.3%3.647 51.394 93.157 00.307 496.7%4.299 91.647 03.682 20.303 4
2.2 内聚能密度
内聚能密度CED(cohesive energy density,)指的是单位体积内凝聚体为克服分子间作用力汽化时所需的能量,是评价分子间作用力大小的物理量,主要反映基团间的相互作用。一般来说,分子中所含基团的极性越大,分子间的作用力就越大,则相应的内聚能密度就越大;反之亦然[12]。表2列出了5个体系下的CED值,可以看出,随着交联度的增加,CED值逐渐增大,表明经过交联反应后体系的体积收缩,同时生成的新的交联键,增大了分子间相互作用力,使得体系的CED值增大。
表2 5个交联度下的CED值和体积
Table 2 CED at five crosslinking conversions
交联度CED/(×108 J·m-3)V/Å380.0%3.80246 34883.3%3.84745 71186.7%3.86145 63893.3%3.88145 61796.7%3.90045 497
2.3 热解过程中势能及分子数变化情况
本文中进一步采用LAMMPS软件分析了交联度为96.7%的体系在高温下的热解情况。图5给出了不同温度下交联体系在热解过程中的势能变化曲线。从图5中的可以看出,5个温度下(4 000、3 500、3 000、2 500和2 000 K)的热解过程中势能变化情况基本相同。首先系统受热,体系势能迅速增加,分别在0.4、0.5、0.6、1.3和1.9 ps时达到峰值,然后逐渐降低,最终热解反应达到平衡。势能曲线的变化反映热解初始阶段需要先吸收能量,热解开始后不断释放能量,热解完成后,能量趋于平稳。温度越高,势能达到最大值的时间越短、最终平衡时的能量也越高,表明温度越高,体系越容易热解,放出的热量也越高。
图5 GAP/N-100体系热解过程中势能变化曲线
Fig.5 Potential energy change curve of GAP/N-100 system
图6为5个温度下体系热解过程中的分子数量的变化情况,从图6中可以看出,温度越高体系热解反应越充分,生成的分子数目也越多。
图6 GAP/N-100体系热解过程中分子数目变化曲线
Fig.6 Curve of molecular number change during pyrolysis of GAP/N-100 system
2.4 初始分解机理
利用LAMMPS输出的产物信息文件结合python脚本对反应过程中的热解机理进行了分析,表3列出了最开始的几个重要反应以及反应发生时间,可以看出,最开始的反应主要是N3叠氮基团的脱落和碳骨架的断裂,这与
等[13]利用质谱法和傅里叶变换红外光谱得到的GAP热解机理的实验结果一致。
表3 4 000 K温度下GAP/N-100体系最开始的反应
Table 3 The first reactions of GAP/N-100 system at 4 000 K
反应方程式时间 /fsC113H192N94O38→C113H192N91O38 + N399.6,103.6,104.6C113H192N94O38→C113H191N94O38+H 108.5C113H192N94O38→C113H191N94O37+HO109.5C113H191N94O38→C21H36N21O8+C92H155N73O30110.5C113H192N94O38→C18H31N18O7+C95H161N76O31110.6H+C92H155N73O30→C92H156N73O30110.9C113H192N94O38→C113H192N91O38 + N3116.2C113H192N94O38→C102H173N82O34+C11H19N12O4116.5
2.5 反应动力学参数
对GAP/N-100体系势能曲线衰减部分进行指数衰减拟合[14],如式(8)所示。所得的反应速率见表4。
表4 GAP/N-100体系各个温度下反应速率
Table 4 Reaction rate of GAP/N-100 system at each temperature
温度/K反应速率k/ps-14 0000.067 393 5000.062 243 0000.056 722 5000.051 952 0000.044 59
U(t)=U(∞)+ΔU*exp[-k*(t-t0)]
(8)
式(8)中:t0为势能到达峰值时间;U(t0)是t=t0时的势能;U(∞)为势能的平衡值;ΔU为U(∞)与U(t0)的差值。
根据式(9)对得到的反应速率进行线性拟合,如图7所示,按照式(9)得到GAP/N-100体系活化能[15]。
图7 GAP/N-100热解反应速率
Fig.7 GAP/N-100 pyrolysis reaction rates
(9)
式(9)中:A为指前因子;Ea为活化能;R为理想气体常数。经过计算得到GAP/N-100热解反应的活化能Ea为13.411 kJ/mol,指前因子A为0.099 1 /ps-1。
2.6 主要产物分析
GAP基含能推进剂的受热分解是一个极其复杂的过程,Mishra等[16]对GAP的热分解产物进行过分析,发现了20多种产物。本文中使用编写的python脚本对LAMMPS输出文件进行了分析,得到了GAP/N-100体系分解的主要产物,主要为H2、N2、NH3、H2O以及中间产物CH2O。
2.6.1 N2数目变化
图8为热解过程中N2数目变化情况,可见反应的初始阶段,N2的数量迅速增加,表明叠氮基团化学键的断裂生成N2是反应初始阶段的主要反应过程。这一结果与陈智群等[17]在研究GAP热分解过程中所得到的第一步放热分解过程的结果是一致的。
图8 GAP/N-100热解过程中N2数目变化情况
Fig.8 Changes in the number of N2 during pyrolysis of GAP/N-100
在5个温度下,最终N2的数目稳定在大约275左右,表明反应已经接近完全。此外,温度越高,N2的生成速率越快。在2 000 K的温度下,N2的生成速率明显低于其他4个温度条件。在前4个温度下,N2的数量在约20 ps内达到峰值,而在2 000 K的情况下,需要50 ps之后才达到峰值。这说明温度是影响化学反应进程的重要因素。
2.6.2 H2数目变化
图9为热解过程中H2数目变化情况,与N2的情况不同,H2的数量变化呈现出明显的分级现象。在2 000 K的温度下,很少观察到H2的生成。而在4 000、3 500和3 000 K温度下,H2的最终数目在大约400左右波动。涉及H2生成的反应是主要为:NH2+H→N2+H2,NH3+N2H→N2+H2+NH2,NH4→NH2+H2等。
图9 GAP/N-100热解过程中H2数目变化情况
Fig.9 Changes in the number of H2 during pyrolysis of GAP/N-100
2.6.3 NH3数目变化
图10为热解过程中NH3数目变化情况,可见在4 000、3 500、3 000以及2 500 K的温度下,NH3的最终数目在大约85左右波动。而在2 000 K的温度下,NH3的最终数目为44,曲线仍呈上升趋势,表明反应尚未完全转化。涉及NH3生成的反应主要有:NH4→NH3+H,NH2+H→NH3,N2+NH4→NH3+NH2等。
图10 GAP/N-100热解过程中NH3数目变化情况
Fig.10 Changes in the number of NH3 during pyrolysis of GAP/N-100
2.6.4 H2O数目变化
图11为热解过程中H2O数目变化情况,可见随着温度的升高,H2O的生成速率增加。在4 000 K和3 500 K的温度下,H2O的最终数目在大约110左右波动。其余3个温度下生成的H2O数目小于100。涉及H2O生成的反应主要有:NH3+OH→H2O+NH2,H+OH→H2O,主要为游离的H+与OH-与体系内的分子片段发生反应。
图11 GAP/N-100热解过程中H2O数目变化情况
Fig.11 Changes in the number of H2O during pyrolysis of GAP/N-100
2.6.5 CH2O数目变化
图12为热解过程中CH2O数目变化情况,可见CH2O数目先增加后减小,表明反应生成的CH2O在与其他分子片段继续反应后被消耗。整个过程随着温度的增加而加快。 在4 000、3 500和3 000 K的温度下,初始生成的CH2O几乎被完全消耗。而2 500和2 000 K仅消耗了一部分,反应进程相对缓慢。
图12 GAP/N-100热解过程中CH2O数目变化情况
Fig.12 Changes in the number of CH2O during pyrolysis of GAP/N-100
3 结论
本文中使用MS通过自编脚本建立了GAP/N-100交联体系的模型,并预测了不同交联度的力学性能和结合能,采用LAMMPS软件在ReaxFF-lg力场条件下对GAP/N-100交联体系的受热机理进行了模拟,得到以下结论:
1) 随交联度的增大,交联键数目增多,体积收缩加剧,CED值增大,生成分子间作用力增大,E、G、K值越大。
2) 温度越高,体系反应越剧烈;热解最开始的反应主要是N3叠氮基团的脱落和碳骨架的断裂。
3) 热解反应的活化能Ea为13.411 KJ/mol,指前因子A为0.099 1/ps-1,热解的主要产物有N2、H2、H2O以及NH3。
[1] K.HUSSEIN A,ELBEIH A,ZEMAN S.The effect of glycidyl azide polymer on the stability and explosive properties of different interesting nitramines[J].RSC Advances,2018,8(31):17272-17278.
[2] 王鑫,黄振亚,刘丽平.多炔基化合物与GAP的固化反应[J].含能材料,2015,23(7):633-637.WANG Xin,HUANG Zhenya,LIU Liping.Curing of glycidyl azide polymer with multiple acetylene-terminated compound[J].Chinese Journal of Energetic Materials,2015,23(7):633-637.
[3] 胡义文,周伟良,肖乐勤,等.低交联度GAP弹性体的制备及其宽温力学性能[J].南京理工大学学报,2017,41(2):226-231.HU Yiwen,ZHOU Weiliang,XIAO Leqin,et al.Preparation of GAP elastomer under low cross-linking degree and its mechanical property in wide temperature range[J].Journal of Nanjing University of Science and Technology,2017,41(2):226-231.
[4] 赵昱,张晓宏,张伟,等.GAP黏合剂胶片力学性能的影响因素[J].火炸药学报,2016,39(1):79-83.ZHAO Yu,ZHANG Xiaohong,ZHANG Wei,et al.Factors affecting the mechanical properties of GAP binder films[J].Chinese Journal of Explosives &Propellants,2016,39(1):79-83.
[5] AGAWANE N T,SOMAN R R,WAGH R M,et al.Optimization of curing agents for linear difunctional glycidyl azide polymer (GAP),with and without isocyanate,for binder applications[J].Central European Journal of Energetic Materials,2018,15(1):206-222.
[6] LU Y,SHU Y,LIU N,et al.Theoretical simulations on the glass transition temperatures and mechanical properties of modified glycidyl azide polymer[J].Computational Materials Science,2017,139:132-139.
[7] MA S,DU W,LUO Y.Simulation of GAP/HTPB phase behaviors in plasticizers and its application in composite solid propellant[J].e-Polymers,2018,18(6):529-540.
[8] WU Y,CHAI K,CAI L,et al.Preparation,characterization and molecular dynamics of novel pre-plasticized nitrocellulose composite microspheres containing burning rate catalysts[J].Journal of Energetic Materials,2021,39(2):228-245.
[9] VAN D A C T,DASGUPTA S,LORANT F,et al.ReaxFF:A reactive force field for hydrocarbons[J].The Journal of Physical Chemistry A,2001,105(41):9396-9409.
[10] LIU L,LIU Y,ZYBIN S V,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.
[11] 杨潞霞,郭志婧,郭丽媛,等.基于Perl语言的交联环氧树脂分子模型的构建和性能仿真[J].计算机与应用化学,2017,34(2):172-176.YANG Luxia,GUO Zhijing,GUO Liyuan,et al.The construction of the model and properties simulation of crosslinking epoey resin based on Perl[J].Computers and Applied Chemistry,2017,34(2):172-176.
[12] LU Y,SHU Y,LIU N,et al.Molecular dynamics simulations on ε-CL-20-based PBXs with added GAP and its derivative polymers[J].RSC Advances,2018,8(9):4955-4962.
[13] 
J.Thermal decomposition of glycidyl azide polymer by direct insertion probe mass spectrometry[J].Journal of Analytical and Applied Pyrolysis,2002,63(2):327-338.
[14] ZENG T,YANG R,LI J,et al.Thermal decomposition mechanism of nitroglycerin by ReaxFF reactive molecular dynamics simulations[J].Combustion Science and Technology,2021,193(3):470-484.
[15] 苑晓峰,张树海,苟瑞君,等.CL-20/HTPB热分解机理的ReaxFF/lg反应分子动力学模拟[J].装备环境工程,2022,19(10):1-11.YUAN Xiaofeng,ZHANG Shuhai,GOU Ruijun,et al.ReaxFF/lg reaction molecular dynamics simulation of thermal decomposition mechanism of CL-20/HTPB[J].Equipment Environmental Engineering,2022,19(10):1-11.
[16] 施明达.GAP与GAP推进剂的研究进展[J].火炸药学报,1994(1):9-16.SHI Mingda.Research progress on GAP and GAP propellants[J].Chinese Journal of Explosives &Propellants,1994(1):9-16.
[17] 陈智群,刘艳,刘子如,等.GAP热分解动力学和机理研究[J].固体火箭技术,2003(4):52-54.CHEN Zhiqun,LIU Yan,LIU Ziru,et al.An investigation on the kinetics and mechanism of the therma ldecomposition of GAP[J].Journal of Solid Rocket Technology,2003(4):52-54.