蒽热生焦初期的反应分子动力学模拟

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

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

针状焦具有机械强度高、热膨胀系数小、电阻率低、耐腐蚀等优点，主要用于高功率、超高功率石墨电极和负极材料的制造。因而，针状焦是一种高附加值产品，其市场价格相比普通焦高得多。

针状焦是大分子稠环芳烃(PAHs)的类石墨晶体[1]，微晶以层状堆积，单元取向度高。按照原料来源，针状焦可分为石油系针状焦和煤系针状焦。1950年美国大湖炭素公司首先发明了石油系针状焦[2]，开启了针状焦的应用与研究。1965年Brooks等提出的炭质中间相理论[3-4]为炭素科学提供了理论研究基础，也推动了针状焦生产技术的发展。中国对针状焦的研究始于20世纪70年代末，中国石化石油化工科学研究院先后完成了石油系针状焦制备的小试、中试和工业试验，确定原料和焦化工艺，建立了相应的测试方法[5-6]。

目前，中国针状焦生产还存在质量不稳定、焦炭强度低、粉焦多等问题，与国际先进水平仍有差距[7]。为了提升针状焦品质，需要从针状焦生成机理着手研究，考察针状焦生成的反应过程，明确对提升针状焦品质有利的影响因素。针状焦生成过程主要遵循液相炭化机理，在700～780 K的高温下，相对分子质量较小的芳香烃(如：萘、蒽、菲、芘)经过脱氢缩合逐渐生成大平面稠环芳烃[8-14]。

为深入研究针状焦的生成机理，选择蒽作为模型化合物，研究蒽热生焦过程的反应分子动力学。在液相炭化的反应温度下，蒽分子很难直接脱氢进而均裂生成自由基，比较容易发生氢转移反应，进而引发自由基反应[15]。为考察蒽热生焦过程的反应分子动力学和热力学，采用分子模拟方法研究蒽液相炭化生成四聚物的反应历程。

1 计算模拟方法

针状焦生成过程是一个复杂的物理化学过程，反应阶段主要是原料受热后逐渐缩合生成低聚物。因此，需要分阶段进行建模和模拟，本研究主要剖析蒽生成二聚物、二聚物生成三聚物和二聚物生成四聚物等3个阶段的蒽热生焦反应过程。

1.1 创建模型

采用Materials Studio 2017 R2软件(简称MS)创建蒽分子的模型。首先选用MS中的Amorphous Cell模块创建一个包含200个蒽分子的三维周期性初始模型，其密度初步设定为0.001 g/cm3；然后选用Forcite模块对初始模型进行结构弛豫，力场选用Compass II，精度为Medium。先选择正则系综(NVT)在2000 K下模拟分子动力学200 ps，促使蒽分子充分散开，在1500～2500 K温度区间选择NVT系综模拟退火(Annealing)；然后选择等温等压系综(NPT)在770 K和1 GPa下模拟分子动力学200 ps，在500～1000 K温度区间选择NPT系综进行退火模拟；再选择NPT系综在770 K和0.2 MPa下模拟分子动力学200 ps，选择NVT系综在770 K和0.2 MPa下模拟分子动力学200 ps；最终得到密度为0.799 g/cm3的纯蒽模型。采用相同方法创建含有100个蒽与50个蒽二聚物的混合体系模型和含有100个蒽二聚物的模型，结构弛豫后的密度分别为0.846 g/cm3和0.948 g/cm3。3种模型如图1所示。

1.2 反应分子动力学模拟

将模型从MS软件中导出，导入MAPS软件(Materials and Processes Simulations 4.1，SCIENOMICS公司产品)中。选用MAPS软件的Lammps[16]模块，采用反应分子动力学方法ReaxFF[17]模拟蒽热生焦过程中的反应分子动力学，力场参数选自含有C、H元素的反应力场文件CH_aromatics.ff[18-20]。针状焦生产温度一般在700～780 K[21-27]。根据温度加速反应动力学理论[28]将模拟温度提升至2600 K[29]。反应分子动力学模拟选择正则系综(NVT)，步长0.1 fs，运行时长200 ps，每1 ps取1帧画面。使用MAPS中ReaxFF分析插件对反应分子动力学模拟结果进行分析。

1.3 量子化学计算分析

由于采用升温加速动力学的方法，模拟中可能会发生低温下不会发生的反应而出现失真现象。因此，有必要采用量子化学方法对某些关键反应进行计算分析。量子化学方法的计算分析选用MS软件中Dmol3模块，通过Dmol3模块可以优化模型分子的几何结构，搜索基元反应的过渡态，计算反应能垒和反应焓变。计算基于密度泛函理论(Density Function Theory)的量子化学方法进行，具体方法选用基于广义梯度近似(GGA)的PBE (Perdew-Burke-Ernzerhof)泛函，在大数值基组DNP(双数值轨道基组+p轨道极化函数)水平上选择全电子(All electron)进行计算。自洽场(Self-Consistent Field, 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 纯蒽模型反应分子动力学模拟

2.1.1 纯蒽模型反应分子动力学

图2为纯蒽模型反应分子动力学模拟结果。其中，图2(a)和(b)分别为模拟100 ps和200 ps时模型的结构状态；图2(c)和(d)分别为将图2(a)和(b)中未反应的蒽分子隐藏后，模拟100 ps和200 ps时生成的分子和自由基的结构状态。不同颜色对应不同种类的分子或自由基，颜色种类越多说明生成的分子和自由基种类越多。由图2纯蒽模型反应分子动力学模拟结果可知：纯蒽热生焦过程发生了缩合、氢转移和重排等反应，很少发生开环反应；100 ps时反应生成了·C28H19、·C14H9、·C14H11等自由基；200 ps时反应生成了·C28H20、·C28H19、·C14H9、·C14H11、C28H18等自由基或分子，·C28H19、·C14H11数量增多。

图2 蒽反应分子动力学模拟结果Fig.2 Reactive molecular dynamics simulationof anthracene coking process(a) Structure diagram at 100 ps; (b) Structure diagram at 200 ps;(c) Molecules and free radicals generated in 100 ps;(d) Molecules and free radicals generated in 200 psVarious colors stand for different molecular or radical types

使用MAPS中ReaxFF分析插件对蒽热生焦反应分子动力学模拟结果进行分析。蒽热生焦反应模拟中蒽分子数量和生成物(·C28H20、·C14H9、·C14H11、H2、C28H18等)数量随模拟时间的变化曲线如图3所示。由图3可观察到，蒽分子数量随着模拟时间的增长而逐渐减少，从200个减少到160个左右，总共消耗了大约40个蒽分子，模拟200 ps时转化率约为20%。由图3还可以观察到：·C28H20数量最多时只有3个，约占原料分子数量的1.5%，对后续反应影响很小；·C14H9、·C14H11和·C28H19数量最多时分别有7、10和5个，数量较多，说明这些自由基稳定性较强，对后续反应影响较大；反应生成了1个C28H18分子，为蒽分子脱去2个氢原子后生成的二聚物。

根据模拟所取的帧数统计反应次数，帧数越多，统计到的反应次数可能会越多。模拟中每1 ps取一帧，模拟时长200 ps共取得200帧。其中，发生次数较多的反应如表1所示。由表1可知：蒽热生焦模拟过程中反应次数最多的反应是2个蒽分子的缩合反应(A1)，其生成的·C28H20自由基热稳定性差，通过逆反应(A2)快速生成C14H10；其次是2个蒽分子之间发生氢转移反应(A3)，生成了·C14H9和·C14H11，其稳定性都较强；蒽分子和自由基 ·C14H9结合生成·C28H19的反应(A4)次数也较多，生成·C28H19。

2.1.2 量子化学分析

蒽分子参与反应的部位主要是1号位、2号位和9号位的碳，很少发生在11号位[30-31]。根据量子化学第一性原理计算蒽分子结构性质和净电荷分布，结果如图4所示。由图4可知：蒽分子中1号位、2号位和9号位C-H键的键长相差不大，键长无法体现1号位、2号位和9号位上C-H键强度的差异；而9号位的C比1号位和2号位的C上分布的净电荷稍多，表明9号位的C活性稍强；3个位置上的H电荷分布相当，表明3个H原子的活性相同。Sasaki等[30-31]计算蒽分子中1号位、2号位和9号位C上的自由价分别为0.459、0.408、0.520，认为蒽分子9号位活性最强。综合分析，9号位C活性最强，易发生脱氢反应或者氢转移反应生成自由基，与文献[32]结论一致。

图3 蒽热生焦反应中反应物和生成物数量随模拟时间的变化曲线Fig.3 Numbers of anthracene molecules and ·C28H20, ·C14H9, ·C14H11, ·C28H19, C28H18 versus simulation time(a) Anthracene; (b) ·C28H20, ·C14H9, ·C14H11, ·C28H19 and C28H18

表1 模拟过程中部分反应统计Table 1 Statistics of partial reactions during the simulation

通过量子化学计算可知：1号位、2号位和9号位C-H键的均裂能分别为483.25、484.70、483.60 kJ/mol，均裂能都很高。这说明3个位置上C-H键断裂的难度相当，都需要很高的温度。依据Poutsma规则[33-34]，芳烃很难在780 K以下均裂生成自由基，一般通过氢转移反应引发生成自由基。

图4 蒽分子结构性质和电荷分布Fig.4 Anthracene structure features and Mülliken charge(a) Carbon atom number; (b) Bond length; (c) Mülliken charge of atom

蒽受热以后可能发生多种反应[25]。为避免采用升温加速动力学方法模拟蒽热生焦反应可能导致的模拟失真问题，采用量子化学对反应分子动力学模拟过程中发生的各种反应进行动力学和热力学分析，通过搜索过渡态方法计算反应的能垒和焓变，结果如表2所示。

由表2可知，2个蒽分子之间首先发生氢转移反应(B1)，引发生成自由基。其反应能垒较高，为321.81 kJ/mol，反应焓变296.89 kJ/mol，为吸热反应；因此该反应在温度较高时才能够发生。蒽分子间氢转移反应能垒是蒽热生焦反应过程中的第一道反应能垒，对后续反应影响较大。氢转移反应可引发自由基，因而适当的供氢能够稳定自由基的反应性，使反应中间体保持良好的流动性和溶解性，促使片层分子有序排列[35]。

蒽分子与自由基·C14H9结合生成·C28H19(B2)，反应能垒很低，为10.81 kJ/mol，反应焓变为-165.29 kJ/mol，因而是易于发生的放热反应，且对蒽二聚物的生成有较大影响。生成的·C28H19可能与周围的蒽分子碰撞而失去1个氢原子，生成热力学性质更为稳定的二聚物C28H18(B3)。其反应能垒为207.99 kJ/mol，反应焓变分别为17.29 kJ/mol，反应能垒较高，但在受热时可以发生。

2个蒽分子缩合生成自由基·C28H20的反应(B4)为可逆反应，其正反应的能垒为197.47 kJ/mol，焓变为194.99 kJ/mol。根据微观可逆性原理，正反应和逆反应经过同一活化体，因此该反应逆反应的能垒为2.46 kJ/mol，焓变为-194.99 kJ/mol。该逆反应能垒非常低，且为放热反应，易于发生。纯蒽模型反应分子动力学模拟中正、逆反应次数相近，因而检测到的·C28H20数量非常少。这也说明·C28H20稳定性很差，很容易发生逆反应生成蒽。此外，·C28H20也可能发生分子内氢转移反应(B5)，反应能垒较低，为116.50 kJ/mol；反应焓变为-21.59 kJ/mol，为放热反应，因而也较容易进行。

表2 蒽热生焦反应能量分析Table 2 Energy analysis of anthracene pyrolysis

2.1.3 蒽缩合生成二聚物的反应历程

根据纯蒽模型反应分子动力学模拟和量子化学分析结果，可以推断蒽热生焦过程中生成二聚物的主要反应历程如图5所示，图中数据为反应能垒。蒽先通过分子间氢转移反应引发生成自由基·C14H9；·C14H9与蒽分子结合生成·C28H19；·C28H19被周围的蒽分子夺去1个氢原子生成稳定的二聚物C28H18。反应过程中，自由基引发阶段分子间的氢转移反应能垒最高，是蒽热生焦过程初期的控速步骤；相比自由基引发反应，自由基链传递反应能垒较低，较容易进行。

2.2 混合模型反应分子动力学模拟

图6为混合模型反应分子动力学模拟结果。其中，图6(a)为模拟100 ps时生成的分子和自由基的结构状态；图6(b)为模拟200 ps时生成的分子和自由基的结构状态。图中不同颜色对应不同种类的分子或自由基；颜色种类越多说明生成的分子和自由基种类越多。由图6可知：混合模型模拟过程中发生了缩合、氢转移、重排和开环等反应；100 ps时生成了·C28H17、·C28H19、·C28H20、·C14H9、·C14H11、·C42H27、C42H26等自由基或分子；200 ps时生成了·C28H17、·C28H19、·C14H9、·hC14H11、·C42H28等自由基，·C28H17、·C14H11数量有所增多。

图5 蒽缩合生成二聚物的反应历程Fig.5 Reactive process of anthracene to dimerThe data in the diagram are reactive energy barriers.

图6 混合模型反应分子动力学模拟结果Fig.6 Reactive molecular dynamics simulation ofmix model coking process(a) 100 ps; (b) 200 psVarious colors stand for different molecular or radical types.

模拟中每1 ps取1帧，模拟时长200 ps共取得200帧。使用MAPS中ReaxFF分析插件对反应分子动力学模拟结果进行分析，其中部分反应次数统计如表3所示。蒽和二聚物分子数量随模拟时间的变化曲线如图7所示。

由表3可知：反应次数较多的反应是蒽分子的缩合反应及其逆反应(R1～R6)；缩合生成的·C42H28、·C28H20、·C56H36自由基热稳定性差，很容易发生逆反应；正、逆反应次数相当，累积生成的缩合产物分子数量很少，对后续反应影响很小。对后续反应影响较大的主要是分子间氢转移反应(R7和R9)。氢转移引发生成的自由基·C28H17和 ·C14H9可能与周围蒽分子结合生成·C42H27(R8)和·C28H19(R10)，也可能与周围的二聚物分子反应生成三聚物C42H26(R11)。

由图7可观察到，蒽分子数量随着模拟时间的增长而逐渐减少，从100个减少到85个左右，总共消耗了大约15个分子，模拟200 ps时转化率约为15%。蒽二聚物分子数量从50个减少到40个左右，总共消耗了大约10个分子，模拟200 ps时转化率约为20%。

表3 混合模型模拟过程中部分反应统计Table 3 Statistics of partial reactions during the simulation

图7 蒽和二聚物分子数量随模拟时间的变化曲线Fig.7 Molecular numbers of anthracene anddimer versus simulation time

根据混合模型反应分子动力学结果可以观察到蒽二聚物缩合反应的反应路线，如图8所示。由图8 可知，蒽二聚物先通过氢转移反应引发生成自由基·C28H17，然后自由基·C28H17与周围的二聚物分子反应生成三聚物。

图8 混合模型反应路线Fig.8 Reactive route of mix model in reactive molecular dynamics simulation

2.3 蒽二聚物模型反应分子动力学模拟

图9为蒽二聚物模型反应分子动力学模拟结果。其中，图9(a)为模拟100 ps时生成的分子和自由基的结构状态；图9(b)为模拟200 ps时生成的分子和自由基的结构状态。图中不同颜色对应不同种类的分子或自由基；颜色种类越多说明生成的分子和自由基种类越多。由图9可以观察到：模拟过程中发生了缩合、氢转移、重排和开环等反应；100 ps时生成了·C28H17、·C28H19等自由基；200 ps时生成了·C28H17、·C28H19、·C56H35、·C56H36等自由基。

图9 蒽二聚物模型反应分子动力学模拟结果Fig.9 Reactive molecular dynamics simulation ofanthracene dimer coking process.(a) Molecules and free radicals generated in 100 ps;(b) Molecules and free radicals generated in 200 psVarious colors stand for different molecular or radical types.

使用MAPS中ReaxFF分析插件对反应分子动力学模拟结果进行分析。模拟中每1 ps取1帧，模拟时长200 ps共取得200帧。其中部分反应次数统计如表4所示。蒽二聚物分子数量随模拟时间的变化曲线如图10所示。

由表4可知，反应次数最多的是2个蒽二聚物分子缩合成·C56H36(D1)。但·C56H36属于热稳定性较差的自由基，很快发生逆反应(D2)，因而对后续反应影响很小。对后续反应影响较大的主要是分子间氢转移反应(D3)。氢转移引发生成的自由基·C28H17可能与周围二聚物分子结合生成·C56H35(D4)。由图10可观察到，蒽二聚物分子数量随着模拟时间的增长而逐渐减少，从100个减少到80个，总共消耗了20个分子，模拟200 ps时转化率为20%。

表4 蒽二聚物模型模拟过程中部分反应统计Table 4 Statistics of partial reactions during the simulation

图10 蒽二聚物分子数量随模拟时间的变化曲线Fig.10 Molecular number of anthracene dimerversus simulation time

根据混合模型反应分子动力学结果可以观察到蒽二聚物的反应路线，如图11所示。由图11可知，蒽二聚物先通过氢转移反应引发生成·C28H17自由基。然后·C28H17与周围的二聚物分子结合生成四聚体自由基。

图11 蒽二聚物模型反应路线Fig.11 Reactive route of dimer model in reactive molecular dynamics simulation

2.4 蒽热生焦初期的反应路线

根据以上反应分子动力学模拟结果，可以推测蒽在热生焦初期的反应路线，如图12所示。由图12可知：蒽在热生焦初期，首先通过氢转移反应引发生成自由基·C14H9；自由基·C14H9与周围的蒽分子结合生成自由基·C28H19；·C28H19被周围的蒽分子夺去1个氢原子生成较为稳定的二聚物C28H18；二聚物C28H18通过氢转移反应生成二聚体自由基·C28H17；·C28H17与周围的二聚物分子反应，可能生成三聚物C42H26，也可能生成四聚体自由基·C56H35。

图12 蒽热解生焦初期的反应路线Fig.12 Reactive route of anthracene at the initial stage of coking

3 结 论

通过量子化学计算蒽分子结构性质发现，9号位碳的Mülliken净电荷相比1号位与2号位较多，活性最强，易发生脱氢反应或者氢转移反应生成自由基。

采用反应分子动力学模拟方法模拟蒽热生焦初期的反应过程，推断出蒽热生焦初期的主要反应历程：蒽在热生焦初期，首先通过氢转移反应引发生成自由基·C14H9；·C14H9与周围的蒽分子结合生成自由基·C28H19；·C28H19被周围的蒽分子夺去1个氢原子生成二聚物C28H18；二聚物C28H18通过氢转移反应生成二聚体自由基·C28H17；·C28H17与周围的蒽二聚物分子反应，可能生成三聚物 C42H26，也可能生成四聚体自由基·C56H35。重复进行脱氢缩合反应，稠环芳香核随之增大，最终生成大平面稠环分子。