丙炔基双酚A醚硼聚合物热解过程的ReaxFF分子动力学模拟

高宁，王一超，刘育红

（西安交通大学化学工程与技术学院，陕西 西安 710049）

引 言

随着科学技术的发展，耐高温聚合物在交通、能源、建筑和航空航天等领域中的应用越来越广泛[1-3]。传统的耐高温聚合物普遍存在熔融温度高，成型工艺差等缺点。因此，设计新型具有良好工艺性能和优异耐热性能的聚合物成为了研究人员关注的焦点[4-6]。从聚合物结构特性分析，分子链中含有炔基的树脂因其固化前黏度小、固化过程无小分子产生、高温残炭率高等优点备受关注[7-8]。此外，在聚合物分子链中引入无机元素硼，因炭化过程中硼元素的促石墨化作用及高温下B2O3保护层的形成，可明显改善材料的耐高温氧化性能[9-11]。基于上述两方面的考虑，设计并制备了新型硼杂化丙炔基树脂——丙炔基双酚A 醚硼聚合物（PB）[12]，其结构式如图1所示。值得注意的是，PB 固化物的玻璃化转变温度（Tg）为416℃，在氮气氛中800℃时的质量保持率可高达74%。PB 树脂上述优异的耐热性能为其作为高性能热固性树脂的应用提供了可能。

图1 PB 的结构式Fig.1 Schematic structure of PB

PB 树脂在高温下发生热分解，释放出小分子气体物质，结构从有机物逐步转变为无定形态无机物。热分解是由有机态向无机态转化的重要步骤，也是影响热解生成物结构的关键因素。因此，对PB 的热解过程及其生成物结构进行研究十分必要。目前，采用热重-红外光谱联用（TG-FTIR）法研究热解过程较为普遍[13-14]。这一方法可以表征热解生成的气体种类，从而推断出可能发生的反应。但由于热解过程涉及的物理和化学变化较为复杂，TG-FTIR 法只局限于对热解过程初步描述，难以准确探究PB树脂在热解过程中发生的结构变化。

分子模拟技术可以通过计算机在分子水平上呈现分子的构象变化并计算材料的性能，弥补实验手段的不足，被称为第3 种认知自然的方法。其中，反应分子动力学模拟（ReaxFF-MD）因其可模拟聚合物的成键、断键过程而受到广泛的关注。尤其在描述极端条件下的反应现象以及复杂的反应过程时，ReaxFF-MD 法可清晰呈现出物质的结构变化，揭示其反应机理。Jiang 等[15]采用ReaxFF-MD 法计算得出酚醛树脂的主要热解产物H2O 的反应活化能为（144±28）kJ·mol-1。这一数值与实验值相符，证明ReaxFF-MD 在模拟热解过程上是适用的。并且，通过ReaxFF-MD 法，Jiang 等提出H2O 的生成机理为β-H 的脱除，为酚醛树脂的热解理论提供了新的见解。由此可见，ReaxFF-MD 在描述聚合物的热解过程，揭示其热解机理等方面具有实验分析方法不可比拟的优势[16-17]。

本文在实验分析的基础上，采用ReaxFF-MD模拟方法，从微观角度出发，系统研究PB 耐高温树脂的热解过程。为探究聚合物的热解过程及聚合物热稳定性的成因提供新的借鉴。

1 实验材料和方法

1.1 分析仪器

本实验采用德国 NETZSCH 公司的型号为STA449-F3 的热重分析仪及德国 Bruker 公司的Vertex70 型红外光谱仪，升温速率为10℃·min-1，氮气氛流量40 ml·min-1。热重分析（TGA）采用美国TA 公司的TA2950 热重分析仪测定。测试条件：升温范围，室温～800℃；升温速率，10℃·min-1；氮气氛流量为60 ml·min-1。

1.2 模型构建与优化

1.2.1 分子模型的建立[18-19]

（1）建立PB 分子链的结构模型。对于聚合度的选择，要考虑软件能够承受的计算量以及模型是否合理。综合考虑计算量及精确性，还要兼顾合理性与时效性，最终设定单分子链模型的聚合度为3，其模型结构如图2（a）所示。可以看出，此时模型结构处于较高的能量状态，而且远离平衡状态，因此需要对此结构进行松弛和优化。常用的优化方法主要有分子力学和分子动力学。首先，采用智能最小化的方法对此模型进行能量最小化以获得合理的三维结构，如图2（a）所示，分子力场为Dreiding，当能量收敛至0.0010 kcal·mol-1（1 kcal·mol-1=4.18 kJ·mol-1）时即停止优化。然后，将得到的分子模型在300 K 下进行20 ps 的NVT 动力学平衡，完全优化后的模型如图2（b）所示。通过能量优化前后的对比图片可以看出分子模型经过优化后结构会变得更合理，同时模型的能量也降低。

（2）定义一个三维盒子。在周期性边界条件下，基于Dreiding 力场，对聚合物建立了三维周期边界盒子，盒子尺寸为2.514 nm×2.514 nm×2.514 nm，其中放入8 条硼聚合物分子链。将聚合物的初始密度设定为0.88 g·cm-3（实验值），其结构模型如图3（a）所示。

图2 能量优化前和能量优化之后的PB 分子链模型Fig.2 Snapshot of simulation of PB chain before and after equilibration

（3）再次采用智能最小化的方法对体系进行能量优化。然而，智能最小化方法采用的是分子力学方法，严格地讲是在绝对零度下进行的，当体系能量较小时，通常能够较为准确地优化出分子的最低能量从而得到分子的最优几何构型。当分子体系较为庞大时，分子力学优化的结构可能只是体系的某个位置。为了避免这种情况，需要在较高的温度下，对体系进行分子动力学优化，即先将分子“加热”到较高温度，让分子结构在高温下充分舒展和松弛，然后进行“冷却”。因此，对模型进行退火模拟，令体系从300 K 为起点，步长为50 K，在101325 Pa下对每次的温度进行5 ps 的NVT 动力学模拟。直至温度达到700 K，接着再以50 K 为步长，将体系降至300 K，共重复进行10 次。至此，体系的空间结构得到了逐步优化。最后，在300 K 下进行40 ps的NVT 的分子动力学模拟，当体系的能量基本维持恒定时停止优化。此时，得到了如图3（b）所示的完全优化之后的周期性边界盒子。

图3 在300 K 下能量优化前和能量优化之后的PB 周期边界盒子（灰色：碳；红色：氧；白色：氢；黄色：硼）Fig.3 Snapshot of simulation cell equilibrated at 300 K before equilibration and after equilibration（grey:C; red:O; white:H; yellow:B）

1.2.2 ReaxFF-MD 模拟 将优化完全的聚合物模型体系在300 K 下进行40 ps 的NVT 分子动力学平 衡计算后，分别以50、100 和200 K·ps-1的速率升温至4300 K，通过观察聚合物在热解过程中发生的反应、小分子逸出的数量来研究热解过程并探究其机理。积分算法采用Velocity Verlet 法，控温方法为与热浴耦合的方法（Berendsen 法），模拟步长为0.2 fs。

2 实验结果与讨论

2.1 ReaxFF 升温过程模拟

首先探究不同升温速率对PB 固化物热解开始时间的影响。在此，分别模拟了PB 固化物在50、100 和200 K·ps-1升温速率下聚合物的热解过程。不同升温速率下聚合物的个数随模拟时间的变化如图4所示。不同升温速率下热解开始的时间和温度列于表1中。由表1可以看出，当升温速率从50 K·ps-1升高至200 K·ps-1时，热解反应的开始时间从14.0 ps 缩短到4.8 ps；热解开始的温度由700 K 升高到960 K。较高的升温速率缩短了热解反应的起始时间，促进了热解的进行。升温速率越高，热解开始的温度也越高。

表1 不同升温速率硼聚合物热解开始的时间和温度Table 1 Pyrolysis time and temperature as function of heating rate

通过测试不同升温速率下的DTG（图5）曲线可知，随着升温速率的增加，热解开始的温度也升高。这是由于升温较慢时，体系有充足的反应时间，即在低升温速率下经历的热解时间增加，热解开始的温度较低。由此可见，ReaxFF 反应动力学模拟结果与实验结论一致，进一步证明了ReaxFF 反应动力学在模拟硼聚合物热解过程是适用的。

2.2 TG-FTIR 法表征PB 固化物的热解过程

为了探究PB 固化物的高温热解过程，首先采用TG-FTIR 对热分解过程中逸出的尾气进行在线连续监测，结果见图6。在整个热分解过程，3964～3500 cm-1、1854～1500 cm-1处对应于水蒸气的特征峰[14]及2360～2320 cm-1处属于CO2的特征峰一直存在[13]。其中，水可能以自由水、结合水的形式存在，还可能在热解反应的过程中生成。300℃之前裂解气中的水峰可能由产物吸潮而吸收的自由水或 结合水挥发导致的，300℃之后裂解气中的水应该是热裂解反应生成的水分。从400℃开始，热解尾气的组成与之前有较大差别。从400℃的裂解气的红外谱图中可以看到，裂解产物中除了水蒸气和CO2，在3000～2600 cm-1和1500～1325 cm-1处出现了 甲烷的特征峰[20]，在1850～1600 cm-1处出现了的双键的振动峰，这对应于醛类、酮类及羧酸类等大分子物质[14]。此外，在1845～1500 cm-1处出现了对应于烯烃的特征峰，说明此时产物中的苯环已经开始发生降解[21]。600℃之后，烷烃和烯烃的红外特征峰开始逐渐变小最终消失。在2250～2000 cm-1处代表CO 的特征峰[22]，CO 的生成对应体系中的碳原子逐渐芳香化的过程[14,23]。

图4 在50、100 和200 K·ps-1 升温速率下的硼聚合物（C63H60O12B3）热解过程的碎片分析Fig.4 Average number of fragments as function of heating rate at 50,100,and 200 K·ps-1

图5 不同升温速率时固化后PB 的DTG 曲线Fig.5 Derivative thermogravimetry curves of PB thermosets under different heating rates

图6 PB 固化物的TG-FTIR 谱图Fig.6 TG-FTIR spectra of cured PB

TG-FTIR 测试可以得出热解过程生成的小分子种类，但无法探究热解过程中PB 的结构变化。通过ReaxFF-MD 法，不仅可以验证TG-FTIR 测试结果的可靠性，而且可获得PB 在热解过程中的断键顺序及结构变化。在此，以其中一个结构单元为研究对象，通过观察升温反应过程模拟中聚合物结构的变化有助于研究聚合物的热解机理。首先，以100 K·ps-1的升温速率，将聚合物体系升至4300 K。由图4（b）可知，大约在7.3 ps（730 K）时，聚合 物开始降解。在降解前，可以观察到聚合物链的数量并非一成不变，这是由高温下聚合物链较强的活动性，不断地与相邻的分子链成键、断键导致的。图7描述了PB 固化物在100 K·ps-1热解过程中断键过程。随着温度的升高，大约在12 ps（1300 K）时，开始断裂的是溴丙炔基双酚A 醚的脂肪醚键，这与文献[24]相一致；随后，约在20 ps（2300 K） 时，结构单元中的另外一个脂肪醚键也随之断裂，双酚A结构及双酚A氧自由基结构在模拟体系中大量存在。此外，可以观察到发生芳基醚键断裂的反应较少，这是由于五电子的芳环不稳定；随着温度的进一步升高，约在28 ps（2800 K）时，双酚A结构中的甲基脱除，此时有大量的甲基自由基[28 ps（b）]生成，甲基自由基可与氢结合生成裂解的主要产物CH4；当温度进一步升高时，苯环开始断裂，生成烯烃[31 ps（b）]和酮[36 ps（b）]等有机小分子。

图7 PB 热解过程不同时刻分子片段快照（灰色：碳；红色：氧；白色：氢；黄色：硼）Fig.7 Snapshot of simulated pyolysis of PB at different time（grey :C; red:O; white:H; yellow:B）

可见，ReaxFF-MD 模拟得到的硼聚合物的热解过程与TG-FTIR 的分析结果一致，均是由脂肪醚键的断裂开始热解反应，进而脱除甲基、生成甲烷小分子，随后苯环开始裂解生成烯烃、羧酸类、酮类等有机化合物。因此，ReaxFF-MD 模拟热解过程，不仅验证了TG-FTIR 的结论，而且能够更为直观、从分子水平上“观察”聚合物的结构变化。

2.3 ReaxFF 模拟PB 固化物热解产物的生成机理

通过ReaxFF 反应动力学模拟PB 热解过程，主要产物有CH4、H2O、CO 和H2（图8），前3 种产物的在TG-FTIR 分析中得到了证明。其中，热解生成的 H2是非极性分子，吸收峰很弱，所以在TG-FTIR 谱图中无法体现。此外，TG-FTIR 谱图中显示热解过程中有CO2分子生成，但在模拟过程中无CO2生成。导致这种差异性的原因是在热解过程中模拟体系里没有O2分子，且聚合物结构中没有羧基、酯基的结构[25]。

图8 PB 固化物热解过程生成的小分子CH3•、CH4、CO、H2 和H2O 的分布Fig.8 Average number of fragments of as function of time during pyrolysis of PB thermosets:CH3•,CH4,CO,H2 and H2O

通过追踪产物CH4的生成过程，可得到4 种CH4的主要生成途径（RPW m1—4）（图9）。由此图所列反应式可知，CH4主要是由CH3•和裂解过程生成的碎片结构中的氢原子结合而成，此结论验证了上述推测。其中，氢原子主要来源于羟基[包括双酚A 结构中的羟基（RPW m1）和活泼氢（双酚A结构中苯环的邻位氢（RPW m2）、裂解产物重排生成的苯乙烯的α-H（RPW m3）、裂解产物重排生成的炔烃的α-H（RPW m4）]，而CH3•主要来源于双酚A 结构中甲基的脱离。

H2O 分子的3 种主要的生成机理（RPW w1—3）可参考图10。由图10可见，双酚A 结构中的—OH与H•之间的结合以及—OH 之间的结合是生成H2O分子的主要途径（RPW w2 和RPW w3）。由双酚A中的—OH 生成H2O 的反应机理在环氧树脂中也有提及[26]。此外，H2O 的生成还存在分子间的反应，即与硼原子相连的两个—OH 反应脱去水分子。在热解过程中，频繁出现C—B—O 的三元环，RPW w1同样也是此三元环的一种主要的生成途径。

图11列出了模拟热解过程中H2的5 种主要的生成途径（RPW h1—5）。从图11中可以观察到，H•基团与碎片中的活泼氢发生撞击是生成H2的主要机理（RPW h1、RPW h2、RPW h5）。RPW h3 为双酚A 结构—OH 中的氢与甲基中的氢结合生成分子H2。RPW h4 为分子内反应，即甲基中的两个氢原子形成H2，但在整个热解阶段，此类型的反应出现的频率较低。

同时，通过对CO 分子产生过程的追踪，可以获得CO 分子的5 种主要的生成机理（RPW c1—5）（图12）。CO 是由侧链末端羰基或热解产物中含羰基化合物脱羰基产生的。CO 产物的生成带走了分子结构中的氧，有助于聚合物的成炭过程，很好地验证了实验推测的结果。

图9 ReaxFF-MD 模拟过程中CH4 的主要生成机理Fig.9 Typical mechanism of CH4 formation observed in ReaxFF MD

图10 ReaxFF-MD 模拟过程中H2O 的主要生成机理Fig.10 Typical mechanism of H2O formation observed in ReaxFF MD

图11 ReaxFF-MD 模拟过程中H2 的主要生成机理Fig.11 Typical mechanism of H2 formation observed in ReaxFF MD

图12 ReaxFF 模拟过程中CO 的主要生成机理Fig.12 Typical mechanism of CO formation observed in ReaxFF MD

3 结 论

本文采用将实验分析与ReaxFF-MD 相结合方法，研究了PB 固化物的热解过程并分析了其炭化物的结构。通过观察升温过程中PB 固化物结构的变化，可得到其热解过程中的断键顺序以及小分子 生成机理。研究发现硼聚合物热解过程中首先发生脂肪醚键的断裂，然后双酚A 中—CH3脱除生成CH4，当温度进一步升高时，苯环裂解生成烯烃和酮等有机小分子。在整个热分解过程中，生成的小分子主要有H2、H2O、CO 和CH4。裂解过程中的CH4主要由甲基和裂解过程生成的碎片结构中的活泼氢结合生成，H2主要由H•与碎片中的活泼氢撞击生成的，双酚A 结构中的—OH 和氢的结合是生成H2O 分子的主要途径，而CO 的生成对应裂解的碎片结构中羰基的脱除。模拟结果与实验结论吻合良好，同时从分子角度描述了PB 固化物热解过程。

[1]Liu Jingang(刘金刚),Shen Dengxiong(沈登雄),Yang Shiyong(杨士勇).Recent research and development of matrix resins for high temperature polymer composites [J].Aerospace Mater.Technol.(宇航材料工艺),2013,44 (4):8-13

[2]Hirano K,Asami M.Phenolic resins—100 years of progress and their future [J].React.Funct.Polym.,2013,73 (2):256-269

[3]Cheng R,Zhou Q,Ni L,Chen Y,Wang J.Synthesis and thermal property of boron-silicon-acetylene hybrid polymer [J].J.Appl.Polym.Sci.,2011,119(1):47-52

[4]Chung T C M,Jeong Y,Chen Q,Kleinhammes A,Wu Y.Synthesis of microporous boron-substituted carbon (B/C) materials using polymeric precursors for hydrogen physisorption [J].J.Am.Chem.Soc.,2008,130(21):6668-6669

[5]Devapal D,Packirisamy S,Sreejith K J,Ravindran P V,George B K.Synthesis,characterization and ceramic conversion studies of borosiloxane oligomers from phenyltrialkoxysilanes [J].J.Inorg.Org.Polym.Mater.,2010,20(4):666-674

[6]Zhou Q,Mao Z,Ni L,Chen G,Zhan M S.Chemistry and application of dicarbadodecaborane and its derivatives [J].J.Appl.Polym.Sci.,2007,104(4):2498-2503

[7]Zhang Shijie(张世杰),Zhang Wei(张炜),Guo Yalin(郭亚林),Bai Xia(白侠).Research progress in a new ablate resistant material-polyarylacetylene resin [J].Thermosetting Resin(热固性树脂),2007,22(6):42-46

[8]Zhang Dexiong(张德雄),Zhang Yan(张衍).An overview of polyarylacetylene as the resin matrix for high temperature composites [J].Journal ofSolid Rocket Technol.（固体火箭技术）,2001,24(1):53-59

[9]Liu Y H,Qiang J P,Jing X L.Synthesis and properties of a novel hyperbranched borate [J].J.Polym.Sci.Part A:Polym.Chem.,2007,45(15):3473-3476

[10]Zhong D H,Sano H,Uchiyama Y,Kobayashi K.Effect of low-level boron doping on oxidation behavior of polyimide-derived carbon films [J].Carbon,2000,38(8):1199-1206

[11]Wang J G,Guo Q G,Liu L,Song J R.Study on the microstructural evolution of high temperature adhesives for graphite [J].Carbon,2002,40(13):2447-2452

[12]Gao Ning (高宁).Synthesis,characterization and molecular dynamics simulations of a novel dipropargyl ether of bisphenol A based boron-containing polymer[D].Xi’an:Xi’an Jiaotong University,2014

[13]Atkinson J L,Vyazovkin S.Thermal properties and degradation behavior of linear and branched poly(L-lactide)s and poly(L-lactide-co-glycolide)s [J].Macromol.Chem.Phys.,2012,213(9):924-936

[14]Wang Shurong(王树荣),Liu Qian(刘倩),Zheng Yun(郑赟),Wen Lihua(文丽华),Luo Zhongyang(骆仲泱),Cen Kefa(岑可法).Mechanism study of biomass pyrolysis by thermogravimetric analysis coupled with infrared spectroscopy [J].J.Eng.Thermophys(工程热物理学报),2006,27(2):351-353

[15]Jiang D E,van Duin A,Goddard W A,Dai S.Simulating the initial stage of phenolic resin carbonizationviathe ReaxFF reactive force field [J].J.Phys.Chem.A,2009,113(25):6891-6894

[16]Tack J L,Ford D M.Thermodynamic and mechanical properties of epoxy resin DGEBF crosslinked with DETDA by molecular dynamics [J].J.Mol.Graphics Modell.,2008,26(8):1269-1275

[17]Ho S R,Aoyagi M.Molecular-cluster-assembly method for analysis of high-dimensional structures ofp-cresol chains in phenolic polymers [J].Macromol.Theory Simul.,2013,22(9):443-61

[18]Andrew R.Molecular Modelling Principles and Applications [M].Harlow:Pearson Education Limited,2001:66-70

[19]Tang Bing (唐兵).Conformational statistics and computer simulation of restricted polymer chains [D].Chengdu:Sichuan University,2002

[20]Salema A A,Afzal M T,Motasemi F.Is there synergy between carbonaceous material and biomass during conventional pyrolysis? A TG-FTIR approach [J].J.Anal.Appl.Pyrol.,2014,105(1):217-226

[21]Su S P,Wilkie C A.The thermal degradation of nanocomposites that contain an oligomeric ammonium cation on the clay [J].Polym.Degrad.Stabil.,2004,83(2):347-362

[22]Chen Y F,Chen Z Q,Xiao S Y,Liu H B.A novel thermal degradation mechanism of phenol-formaldehyde type resins [J].Thermochimica Acta,2008,476(1/2):39-43

[23]Maschio G,Koufopanos C,Lucchesi A.Pyrolysis,a promising route for biomass utilization [J].Bioresource Technol.,1992,42(3):219-231

[24]Itoh M,Iwata K,Ishikawa J I,Sukawa H,Kimura H,Okita K.Various silicon-containing polymers with Si(H)―C≡C units [J].J.Polym.Sci.Part A:Polym.Chem.,2001,39(15):2658-2669

[25]Chenoweth K,Cheung S,van Duin A,Goddard W A,Kober E M.Simulations on the thermal decomposition of a poly(dimethylsiloxane) polymer using the ReaxFF reactive force field [J].J.Am.Chem.Soc.,2005,127(19):7192-7202

[26]Diao Z J,Zhao Y M,Chen B,Duan C L,Song S.ReaxFF reactive force field for molecular dynamics simulations of epoxy resin thermal decomposition with model compound [J].J.Anal.Appl.Pyrol.,2013,104(1):618-624