杨海洋， 范启明， 任 强， 周 涵， 王丽新

(中国石化 石油化工科学研究院，北京 100083)

针状焦是大分子稠环芳烃(Polycyclic aromatic hydrocarbons, PAHs)的类石墨晶体[1]，其中，稠环芳烃以层状堆积，石墨状微晶单元取向度高。针状焦具有机械强度高、热膨胀系数小、电阻率低、耐腐蚀等优点，主要用于制造高级石墨和高功率电极。油系针状焦由石油加工而成，属于高附加值石油产品，其市场价格相比普通石油焦高很多。

1950年美国大湖炭素公司首先发明了石油系针状焦[2]。1965年Brooks和Taylor提出的炭质中间相理论[3]推动了整个炭素科学的发展，也推动了针状焦生产技术的发展。我国对针状焦的研究始于20世纪70年代末，中国石化石油化工科学研究院先后完成油系针状焦的实验室生成实验和中试，确定原料和焦化工艺，建立相应的测试方法[4-5]。目前，我国针状焦生产技术与国际先进水平还有差距[6]。

为了提升针状焦生产工艺，需从针状焦生成机理着手研究。研究发现：在针状焦生成前期，相对分子质量较小的芳香烃经过脱氢缩合逐渐生成大平面稠环芳烃[7-8]；而针状焦的生成过程，主要就是液相炭化的过程[9]。Wang[10]分析针状焦原料的组成发现，针状焦主要成分是直链烷烃、萘、蒽、菲、芘，及其衍生物。为了探究针状焦的生成机理，笔者选用分子模拟方法研究萘模型化合物液相炭化生成二聚物的反应分子动力学，并进行量子化学分析，推断出了萘分子生焦的主要反应历程。

1 计算模拟方法

1.1 模型的创建

选用Materials Studio 2017 R2软件(简称MS)中的Amorphous Cell模块创建一个包含40个萘分子的三维周期性初始模型，密度初步设定为0.001 g/cm3，然后选用Forcite模块进行结构弛豫。先用正则系综(NVT)在温度2000 K下模拟200 ps，促使萘分子充分地散开；在1500～2500 K的温度区间，用NVT系综进行退火(Annealing)；用等温等压系综(NPT)在775 K和1 GPa下模拟 200 ps；在500～1000 K的温度区间，用NPT系综进行退火；用NPT系综在775 K和0.2 MPa下模拟200 ps；用NVT系综在775 K和0.2 MPa下模拟200 ps，最终得到密度为0.658 g/cm3的纯萘模型。力场选用Compass II，精度为Medium。

1.2 反应分子动力学模拟

从MS中导出纯萘模型，导入MAPS软件(Materials and processes simulations 4.1，SCIENOMICS公司产品)中，选用Lammps[11]模块，运用反应分子动力学方法ReaxFF[12]模拟萘的热解生焦反应。力场参数选自含有C、H、O 3种元素的反应力场文件CHO.ff[13-14]。针状焦生产温度一般在700～775 K[15]。根据温度加速反应动力学理论[16]，将模拟温度提升至3000 K[17]，选用正则系综(NVT)，步长0.1 fs，运行时长3000 ps，每1 ps取一帧画面。使用MAPS中ReaxFF分析插件对反应分子动力学模拟结果进行分析。

从反应分子动力学模拟的结果中挑选出合理的反应进行量子化学计算分析。量子化学选用MS软件中Dmol3模块，优化模型分子的几何结构，搜索基元反应过渡态，计算反应能垒和反应焓变。量子化学基于密度泛函理论(Density function theory)。具体方法选用基于广义梯度近似(GGA)的PBE泛函，在大数值基组DNP(双数值轨道基组+ p轨道极化函数)水平上进行全电子计算。自洽场(SCF)迭代收敛阈值设为1×10-5Ha(1 Ha=2625.5 kJ/mol)，轨道截断半径为0.37 nm，收敛精度设为：能量2×10-5Ha，受力0.04 Ha/nm，位移5×10-4nm。过渡态搜索采用完全线性同步和二次同步变化方法(Complete LSF/QSF)。

2 结果与讨论

2.1 萘热解生焦反应分子动力学模拟

图1为萘热解生焦反应分子动力学模拟过程中不同时间的结构分布图。其中图(a)、(b)、(c)、(d)分别对应0 ps、1000 ps、2000 ps、3000 ps。按不同分子类型取色，颜色种类越多说明反应生成的分子种类越多。由图1可以观察到，反应模拟过程中，萘分子先出现二联、三联、开环，再裂解生成小分子或自由基。分子种类先逐渐增加，后不断减少。这是由于生成的小分子、自由基相互聚合生成大分子的结果。

图1 萘热解生焦反应分子动力学模拟过程中不同时间的结构分布Fig.1 Molecular structures at different naphthalene cokingtime based on reactive molecular dynamic simulation(a) 0 ps; (b) 1000 ps; (c) 2000 ps; (d) 3000 ps

图2为萘热解生焦反应模拟过程不同模拟时间对应的最大平面结构。平面结构如果很规整，表明在该反应条件下萘生成针状焦的质量较好。由图2可以观察到，萘热解生焦反应过程生成分子的环数逐渐增多，平面部分逐渐扩大。模拟前期有少量萘分子从稳定性较高的六元环转变为稳定性较低的七元环和五元环，再裂解开环，生成小分子或自由基。随着模拟时间的增长，稳定性较低的七元环和五元环再逐渐转变为六元环，生成的平面分子更稳定。

图2 萘热解生焦反应分子动力学模拟过程中不同时间对应的最大平面分子结构Fig.2 Max plane structure at different naphthalene cokingtime based on reactive molecular dynamic simulation(a) 0 ps; (b) 1000 ps; (c) 2000 ps; (d) 3000 ps Carbon Hyderogen

萘的热反应不仅仅是脱氢缩合反应，也有与小分子和自由基发生的加成反应，最终生成含有较多芳环的稠环大分子[18]。

采用MAPS软件中的ReaxFF插件对模拟过程进一步分析，得到模拟过程中产生的小分子或自由基数量随模拟时间的变化曲线，如图3所示。

比较图3中曲线可知，在模拟反应过程中，萘(C10H8)分子数量逐渐减少，在1000 ps左右已经全部反应，而且萘分子数量下降先加速然后减缓。原因在于：反应前期，自由基(H·、·CH3)的量较少，仅靠萘自身裂解，反应速率很慢；随着反应的进行，自由基含量增多，促进萘的热解；当萘热解反应发生到一定程度时，生成了一些活性较强的分子或自由基，与萘竞争活性基团，导致萘热解的速率减缓。

此外，乙炔(C2H2)分子数量变化频率较高、幅度较大，最多时可达到21个。曲线下降表明乙炔参与萘的热解反应，被消耗掉。分子数量变化频率越高，幅度越大，说明乙炔对萘热解的促进作用越大。氢气(H2)分子数量变化频率较高、幅度也较大，数量最多时可达到13个。乙烯(C2H4)分子出现得较晚，开始产生于500 ps，数量最多时可达到12个。甲烷(CH4)分子数量最多时可到5个。另外，模拟过程中生成了自由基·H和·CH3，但因其活性很强，易于反应，在模拟过程中累积的数量很少。这些小分子和自由基会促进萘的热解反应，加速生焦反应进程。

图3 小分子或自由基数目随模拟时长的变化曲线Fig.3 Profiles of small molecule or radical numbers versus simulation time(a) C10H8， C2H4， ·H and H2； (b) CH4， C2H2 and ·CH3

2.2 萘热解生焦反应量子化学分析

2.2.1 结构性质分析

由萘热解生焦反应分子动力学模拟结果分析可知，萘分子参与反应的部位主要是1号和2号位的碳，很少发生开环反应。根据量子化学第一性原理，计算萘分子结构性质和电荷分布如图4所示。

图4 萘分子结构性质和电荷分布Fig.4 Naphthalene structure features and mulliken charge(a) Carbon number; (b) Bond length， nm; (c) Mulliken charge of atom

由图4可知：萘分子中1号位C-H键的键长与2号位C-H键的键长几乎相同(图4(b))，因而键长不能体现2种C-H键强度的差异；1号位C的净负电荷比2号位C分布的更多(图4(c))，表明1号位C的活性更强，与文献[19]结论一致；1号位H与2号位H的Mulliken正电荷分布相当(图4(c))，表明2个H原子的活性相近。

2.2.2 热反应分析

萘受热反应的途径多种多样[20]，主要有脱氢加乙炔反应(HACA机理[21])、重排反应、二聚反应、三聚反应。在MAPS模拟中，最常见的是二聚反应。基于反应分子动力学模拟结果，构建萘受热生成二聚物的反应路线如图5所示。并选用搜索过渡态方法，计算反应过程中的能垒，标于图内。

由图5可观察到，萘在氢自由基作用下脱去1个氢原子，生成·C10H7。·C10H7有2条反应路线：其一，·C10H7按照HACA机理先与乙炔反应生成·C12H9；再与氢自由基反应，脱去1个氢气分子，生成苊烯；其二，2个·C10H7反应生成1,1′-联萘C20H14，并转化生成1个中间体，再脱氢生成苝C20H12。能垒越低，反应越容易进行。因而萘生成·C10H7后，很容易生成C20H14，再脱氢缩合生成苝，最终生成大平面稠环分子。

图5 萘生焦反应路线Fig.5 Reactive routes of naphthalene cokingNumbers in the picture are reactive energy barriers

萘分子可能通过热裂解直接脱氢，也可能与自由基反应脱氢。为比较各种反应速率的大小，通过搜索过渡态计算得到反应能垒Ea和反应焓变ΔH，见表1。由表1可知：萘分子直接热解脱氢，反应能垒Ea约为480 kJ/mol，能垒较高；而在自由基H·或·CH3夺氢的情况下，反应能垒Ea不超过50 kJ/mol，脱氢速率大大提高，因此自由基的出现可促进萘热解反应。针状焦生成前期，应当控制温度和压力，抑制气体挥发速率，以保证小分子和自由基能够参与芳核的生长反应。

萘分子脱氢后生成的自由基·C10H7有2种形式，脱氢位置分别为1号位和2号位。这2种自由基两两聚合，生成二聚联萘。分别搜索生成二聚联萘反应的过渡态，计算得到反应能垒Ea和反应焓变ΔH，如表2所示。由表2可知，1号位反应生成1,1′-联萘的反应活化能最低，表明萘分子脱氢后生成1,1′-联萘的反应速率最高。1,1′-联萘再脱氢缩合生成苝；苝再脱氢缩合。芳核逐渐增长，最终生成平面大分子。

表1 萘脱氢反应能量分析Table 1 Reaction energy analysis of naphthalene dehydrogenation

表2 萘缩合反应能量分析Table 2 Reaction energy analysis of naphthalene condensation

2.3 萘生焦反应历程

根据反应分子动力学和量子化学分析，推断萘分子热解生焦的主要反应历程如图6所示。萘分子被自由基夺走氢原子，生成自由基·C10H7；·C10H7两两结合生成1,1′-联萘；1,1′-联萘再脱氢缩合生成苝。苝的反应历程类似于萘，重复进行脱氢缩合反应，芳核也随之增大。

图6 萘生焦反应历程Fig.6 Reactive mechanisms of naphthalene coking

3 结 论

(1)采用反应力场ReaxFF模拟纯萘模型化合物热解生焦的反应分子动力学，萘经过热解、缩合脱氢生成了大平面稠环分子。萘热解过程中生成了一些小分子和自由基，如乙炔、乙烯、氢、甲烷等。这些小分子和自由基促进萘的热解反应，加快反应进程。其中乙炔对萘热解反应的促进作用最大。针状焦生成前期，应当控制温度和压力，抑制气体挥发速率，以保证小分子和自由基能够参与芳核的生长反应。

(2)采用量子化学方法计算萘分子结构性质和电荷分布发现，萘分子1号位和2号位的C-H键强度相当，但1号位C原子上分布的净电荷稍多，活性较强。

(3) 通过量子化学分析萘热解生焦过程的反应分子动力学，表明萘生焦的主要反应历程为：萘分子失氢生成自由基·C10H7；·C10H7两两结合生成1,1′-联萘；1,1′-联萘再脱氢缩合生成苝。苝的反应历程与萘类似，重复进行脱氢缩合反应，芳核也随之增大。