萘熱解生焦過程的反應分子動力學模擬

楊海洋， 范啟明， 任 強， 周 涵， 王麗新

(中國石化 石油化工科學研究院，北京 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模擬中，最常見的是二聚反應?；诜磻肿觿恿W模擬結果，構建萘受熱生成二聚物的反應路線如圖5所示。并選用搜索過渡態方法，計算反應過程中的能壘，標于圖內。

由圖5可觀察到，萘在氫自由基作用下脫去1個氫原子，生成·C10H7?！10H7有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′-聯萘再脫氫縮合生成苝。苝的反應歷程與萘類似，重復進行脫氫縮合反應，芳核也隨之增大。