萘熱解生焦過程的反應(yīng)分子動(dòng)力學(xué)模擬

楊海洋， 范啟明， 任 強(qiáng)， 周 涵， 王麗新

(中國(guó)石化 石油化工科學(xué)研究院，北京 100083)

針狀焦是大分子稠環(huán)芳烴(Polycyclic aromatic hydrocarbons, PAHs)的類石墨晶體[1]，其中，稠環(huán)芳烴以層狀堆積，石墨狀微晶單元取向度高。針狀焦具有機(jī)械強(qiáng)度高、熱膨脹系數(shù)小、電阻率低、耐腐蝕等優(yōu)點(diǎn)，主要用于制造高級(jí)石墨和高功率電極。油系針狀焦由石油加工而成，屬于高附加值石油產(chǎn)品，其市場(chǎng)價(jià)格相比普通石油焦高很多。

1950年美國(guó)大湖炭素公司首先發(fā)明了石油系針狀焦[2]。1965年Brooks和Taylor提出的炭質(zhì)中間相理論[3]推動(dòng)了整個(gè)炭素科學(xué)的發(fā)展，也推動(dòng)了針狀焦生產(chǎn)技術(shù)的發(fā)展。我國(guó)對(duì)針狀焦的研究始于20世紀(jì)70年代末，中國(guó)石化石油化工科學(xué)研究院先后完成油系針狀焦的實(shí)驗(yàn)室生成實(shí)驗(yàn)和中試，確定原料和焦化工藝，建立相應(yīng)的測(cè)試方法[4-5]。目前，我國(guó)針狀焦生產(chǎn)技術(shù)與國(guó)際先進(jìn)水平還有差距[6]。

為了提升針狀焦生產(chǎn)工藝，需從針狀焦生成機(jī)理著手研究。研究發(fā)現(xiàn)：在針狀焦生成前期，相對(duì)分子質(zhì)量較小的芳香烴經(jīng)過脫氫縮合逐漸生成大平面稠環(huán)芳烴[7-8]；而針狀焦的生成過程，主要就是液相炭化的過程[9]。Wang[10]分析針狀焦原料的組成發(fā)現(xiàn)，針狀焦主要成分是直鏈烷烴、萘、蒽、菲、芘，及其衍生物。為了探究針狀焦的生成機(jī)理，筆者選用分子模擬方法研究萘模型化合物液相炭化生成二聚物的反應(yīng)分子動(dòng)力學(xué)，并進(jìn)行量子化學(xué)分析，推斷出了萘分子生焦的主要反應(yīng)歷程。

1 計(jì)算模擬方法

1.1 模型的創(chuàng)建

選用Materials Studio 2017 R2軟件(簡(jiǎn)稱MS)中的Amorphous Cell模塊創(chuàng)建一個(gè)包含40個(gè)萘分子的三維周期性初始模型，密度初步設(shè)定為0.001 g/cm3，然后選用Forcite模塊進(jìn)行結(jié)構(gòu)弛豫。先用正則系綜(NVT)在溫度2000 K下模擬200 ps，促使萘分子充分地散開；在1500～2500 K的溫度區(qū)間，用NVT系綜進(jìn)行退火(Annealing)；用等溫等壓系綜(NPT)在775 K和1 GPa下模擬 200 ps；在500～1000 K的溫度區(qū)間，用NPT系綜進(jìn)行退火；用NPT系綜在775 K和0.2 MPa下模擬200 ps；用NVT系綜在775 K和0.2 MPa下模擬200 ps，最終得到密度為0.658 g/cm3的純萘模型。力場(chǎng)選用Compass II，精度為Medium。

1.2 反應(yīng)分子動(dòng)力學(xué)模擬

從MS中導(dǎo)出純萘模型，導(dǎo)入MAPS軟件(Materials and processes simulations 4.1，SCIENOMICS公司產(chǎn)品)中，選用Lammps[11]模塊，運(yùn)用反應(yīng)分子動(dòng)力學(xué)方法ReaxFF[12]模擬萘的熱解生焦反應(yīng)。力場(chǎng)參數(shù)選自含有C、H、O 3種元素的反應(yīng)力場(chǎng)文件CHO.ff[13-14]。針狀焦生產(chǎn)溫度一般在700～775 K[15]。根據(jù)溫度加速反應(yīng)動(dòng)力學(xué)理論[16]，將模擬溫度提升至3000 K[17]，選用正則系綜(NVT)，步長(zhǎng)0.1 fs，運(yùn)行時(shí)長(zhǎng)3000 ps，每1 ps取一幀畫面。使用MAPS中ReaxFF分析插件對(duì)反應(yīng)分子動(dòng)力學(xué)模擬結(jié)果進(jìn)行分析。

從反應(yīng)分子動(dòng)力學(xué)模擬的結(jié)果中挑選出合理的反應(yīng)進(jìn)行量子化學(xué)計(jì)算分析。量子化學(xué)選用MS軟件中Dmol3模塊，優(yōu)化模型分子的幾何結(jié)構(gòu)，搜索基元反應(yīng)過渡態(tài)，計(jì)算反應(yīng)能壘和反應(yīng)焓變。量子化學(xué)基于密度泛函理論(Density function theory)。具體方法選用基于廣義梯度近似(GGA)的PBE泛函，在大數(shù)值基組DNP(雙數(shù)值軌道基組+ p軌道極化函數(shù))水平上進(jìn)行全電子計(jì)算。自洽場(chǎng)(SCF)迭代收斂閾值設(shè)為1×10-5Ha(1 Ha=2625.5 kJ/mol)，軌道截?cái)喟霃綖?.37 nm，收斂精度設(shè)為：能量2×10-5Ha，受力0.04 Ha/nm，位移5×10-4nm。過渡態(tài)搜索采用完全線性同步和二次同步變化方法(Complete LSF/QSF)。

2 結(jié)果與討論

2.1 萘熱解生焦反應(yīng)分子動(dòng)力學(xué)模擬

圖1為萘熱解生焦反應(yīng)分子動(dòng)力學(xué)模擬過程中不同時(shí)間的結(jié)構(gòu)分布圖。其中圖(a)、(b)、(c)、(d)分別對(duì)應(yīng)0 ps、1000 ps、2000 ps、3000 ps。按不同分子類型取色，顏色種類越多說明反應(yīng)生成的分子種類越多。由圖1可以觀察到，反應(yīng)模擬過程中，萘分子先出現(xiàn)二聯(lián)、三聯(lián)、開環(huán)，再裂解生成小分子或自由基。分子種類先逐漸增加，后不斷減少。這是由于生成的小分子、自由基相互聚合生成大分子的結(jié)果。

圖1 萘熱解生焦反應(yīng)分子動(dòng)力學(xué)模擬過程中不同時(shí)間的結(jié)構(gòu)分布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為萘熱解生焦反應(yīng)模擬過程不同模擬時(shí)間對(duì)應(yīng)的最大平面結(jié)構(gòu)。平面結(jié)構(gòu)如果很規(guī)整，表明在該反應(yīng)條件下萘生成針狀焦的質(zhì)量較好。由圖2可以觀察到，萘熱解生焦反應(yīng)過程生成分子的環(huán)數(shù)逐漸增多，平面部分逐漸擴(kuò)大。模擬前期有少量萘分子從穩(wěn)定性較高的六元環(huán)轉(zhuǎn)變?yōu)榉€(wěn)定性較低的七元環(huán)和五元環(huán)，再裂解開環(huán)，生成小分子或自由基。隨著模擬時(shí)間的增長(zhǎng)，穩(wěn)定性較低的七元環(huán)和五元環(huán)再逐漸轉(zhuǎn)變?yōu)榱h(huán)，生成的平面分子更穩(wěn)定。

圖2 萘熱解生焦反應(yīng)分子動(dòng)力學(xué)模擬過程中不同時(shí)間對(duì)應(yīng)的最大平面分子結(jié)構(gòu)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

萘的熱反應(yīng)不僅僅是脫氫縮合反應(yīng)，也有與小分子和自由基發(fā)生的加成反應(yīng)，最終生成含有較多芳環(huán)的稠環(huán)大分子[18]。

采用MAPS軟件中的ReaxFF插件對(duì)模擬過程進(jìn)一步分析，得到模擬過程中產(chǎn)生的小分子或自由基數(shù)量隨模擬時(shí)間的變化曲線，如圖3所示。

比較圖3中曲線可知，在模擬反應(yīng)過程中，萘(C10H8)分子數(shù)量逐漸減少，在1000 ps左右已經(jīng)全部反應(yīng)，而且萘分子數(shù)量下降先加速然后減緩。原因在于：反應(yīng)前期，自由基(H·、·CH3)的量較少，僅靠萘自身裂解，反應(yīng)速率很慢；隨著反應(yīng)的進(jìn)行，自由基含量增多，促進(jìn)萘的熱解；當(dāng)萘熱解反應(yīng)發(fā)生到一定程度時(shí)，生成了一些活性較強(qiáng)的分子或自由基，與萘競(jìng)爭(zhēng)活性基團(tuán)，導(dǎo)致萘熱解的速率減緩。

此外，乙炔(C2H2)分子數(shù)量變化頻率較高、幅度較大，最多時(shí)可達(dá)到21個(gè)。曲線下降表明乙炔參與萘的熱解反應(yīng)，被消耗掉。分子數(shù)量變化頻率越高，幅度越大，說明乙炔對(duì)萘熱解的促進(jìn)作用越大。氫氣(H2)分子數(shù)量變化頻率較高、幅度也較大，數(shù)量最多時(shí)可達(dá)到13個(gè)。乙烯(C2H4)分子出現(xiàn)得較晚，開始產(chǎn)生于500 ps，數(shù)量最多時(shí)可達(dá)到12個(gè)。甲烷(CH4)分子數(shù)量最多時(shí)可到5個(gè)。另外，模擬過程中生成了自由基·H和·CH3，但因其活性很強(qiáng)，易于反應(yīng)，在模擬過程中累積的數(shù)量很少。這些小分子和自由基會(huì)促進(jìn)萘的熱解反應(yīng)，加速生焦反應(yīng)進(jìn)程。

圖3 小分子或自由基數(shù)目隨模擬時(shí)長(zhǎng)的變化曲線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 萘熱解生焦反應(yīng)量子化學(xué)分析

2.2.1 結(jié)構(gòu)性質(zhì)分析

由萘熱解生焦反應(yīng)分子動(dòng)力學(xué)模擬結(jié)果分析可知，萘分子參與反應(yīng)的部位主要是1號(hào)和2號(hào)位的碳，很少發(fā)生開環(huán)反應(yīng)。根據(jù)量子化學(xué)第一性原理，計(jì)算萘分子結(jié)構(gòu)性質(zhì)和電荷分布如圖4所示。

圖4 萘分子結(jié)構(gòu)性質(zhì)和電荷分布Fig.4 Naphthalene structure features and mulliken charge(a) Carbon number; (b) Bond length， nm; (c) Mulliken charge of atom

由圖4可知：萘分子中1號(hào)位C-H鍵的鍵長(zhǎng)與2號(hào)位C-H鍵的鍵長(zhǎng)幾乎相同(圖4(b))，因而鍵長(zhǎng)不能體現(xiàn)2種C-H鍵強(qiáng)度的差異；1號(hào)位C的凈負(fù)電荷比2號(hào)位C分布的更多(圖4(c))，表明1號(hào)位C的活性更強(qiáng)，與文獻(xiàn)[19]結(jié)論一致；1號(hào)位H與2號(hào)位H的Mulliken正電荷分布相當(dāng)(圖4(c))，表明2個(gè)H原子的活性相近。

2.2.2 熱反應(yīng)分析

萘受熱反應(yīng)的途徑多種多樣[20]，主要有脫氫加乙炔反應(yīng)(HACA機(jī)理[21])、重排反應(yīng)、二聚反應(yīng)、三聚反應(yīng)。在MAPS模擬中，最常見的是二聚反應(yīng)?；诜磻?yīng)分子動(dòng)力學(xué)模擬結(jié)果，構(gòu)建萘受熱生成二聚物的反應(yīng)路線如圖5所示。并選用搜索過渡態(tài)方法，計(jì)算反應(yīng)過程中的能壘，標(biāo)于圖內(nèi)。

由圖5可觀察到，萘在氫自由基作用下脫去1個(gè)氫原子，生成·C10H7?！10H7有2條反應(yīng)路線：其一，·C10H7按照HACA機(jī)理先與乙炔反應(yīng)生成·C12H9；再與氫自由基反應(yīng)，脫去1個(gè)氫氣分子，生成苊烯；其二，2個(gè)·C10H7反應(yīng)生成1,1′-聯(lián)萘C20H14，并轉(zhuǎn)化生成1個(gè)中間體，再脫氫生成苝C20H12。能壘越低，反應(yīng)越容易進(jìn)行。因而萘生成·C10H7后，很容易生成C20H14，再脫氫縮合生成苝，最終生成大平面稠環(huán)分子。

圖5 萘生焦反應(yīng)路線Fig.5 Reactive routes of naphthalene cokingNumbers in the picture are reactive energy barriers

萘分子可能通過熱裂解直接脫氫，也可能與自由基反應(yīng)脫氫。為比較各種反應(yīng)速率的大小，通過搜索過渡態(tài)計(jì)算得到反應(yīng)能壘Ea和反應(yīng)焓變?chǔ)，見表1。由表1可知：萘分子直接熱解脫氫，反應(yīng)能壘Ea約為480 kJ/mol，能壘較高；而在自由基H·或·CH3奪氫的情況下，反應(yīng)能壘Ea不超過50 kJ/mol，脫氫速率大大提高，因此自由基的出現(xiàn)可促進(jìn)萘熱解反應(yīng)。針狀焦生成前期，應(yīng)當(dāng)控制溫度和壓力，抑制氣體揮發(fā)速率，以保證小分子和自由基能夠參與芳核的生長(zhǎng)反應(yīng)。

萘分子脫氫后生成的自由基·C10H7有2種形式，脫氫位置分別為1號(hào)位和2號(hào)位。這2種自由基兩兩聚合，生成二聚聯(lián)萘。分別搜索生成二聚聯(lián)萘反應(yīng)的過渡態(tài)，計(jì)算得到反應(yīng)能壘Ea和反應(yīng)焓變?chǔ)，如表2所示。由表2可知，1號(hào)位反應(yīng)生成1,1′-聯(lián)萘的反應(yīng)活化能最低，表明萘分子脫氫后生成1,1′-聯(lián)萘的反應(yīng)速率最高。1,1′-聯(lián)萘再脫氫縮合生成苝；苝再脫氫縮合。芳核逐漸增長(zhǎng)，最終生成平面大分子。

表1 萘脫氫反應(yīng)能量分析Table 1 Reaction energy analysis of naphthalene dehydrogenation

表2 萘縮合反應(yīng)能量分析Table 2 Reaction energy analysis of naphthalene condensation

2.3 萘生焦反應(yīng)歷程

根據(jù)反應(yīng)分子動(dòng)力學(xué)和量子化學(xué)分析，推斷萘分子熱解生焦的主要反應(yīng)歷程如圖6所示。萘分子被自由基奪走氫原子，生成自由基·C10H7；·C10H7兩兩結(jié)合生成1,1′-聯(lián)萘；1,1′-聯(lián)萘再脫氫縮合生成苝。苝的反應(yīng)歷程類似于萘，重復(fù)進(jìn)行脫氫縮合反應(yīng)，芳核也隨之增大。

圖6 萘生焦反應(yīng)歷程Fig.6 Reactive mechanisms of naphthalene coking

3 結(jié) 論

(1)采用反應(yīng)力場(chǎng)ReaxFF模擬純萘模型化合物熱解生焦的反應(yīng)分子動(dòng)力學(xué)，萘經(jīng)過熱解、縮合脫氫生成了大平面稠環(huán)分子。萘熱解過程中生成了一些小分子和自由基，如乙炔、乙烯、氫、甲烷等。這些小分子和自由基促進(jìn)萘的熱解反應(yīng)，加快反應(yīng)進(jìn)程。其中乙炔對(duì)萘熱解反應(yīng)的促進(jìn)作用最大。針狀焦生成前期，應(yīng)當(dāng)控制溫度和壓力，抑制氣體揮發(fā)速率，以保證小分子和自由基能夠參與芳核的生長(zhǎng)反應(yīng)。

(2)采用量子化學(xué)方法計(jì)算萘分子結(jié)構(gòu)性質(zhì)和電荷分布發(fā)現(xiàn)，萘分子1號(hào)位和2號(hào)位的C-H鍵強(qiáng)度相當(dāng)，但1號(hào)位C原子上分布的凈電荷稍多，活性較強(qiáng)。

(3) 通過量子化學(xué)分析萘熱解生焦過程的反應(yīng)分子動(dòng)力學(xué)，表明萘生焦的主要反應(yīng)歷程為：萘分子失氫生成自由基·C10H7；·C10H7兩兩結(jié)合生成1,1′-聯(lián)萘；1,1′-聯(lián)萘再脫氫縮合生成苝。苝的反應(yīng)歷程與萘類似，重復(fù)進(jìn)行脫氫縮合反應(yīng)，芳核也隨之增大。