趙 浩,張英佳,黃佐華

(西安交通大學(xué) 動力工程多相流國家重點實驗室,陜西 西安 710049)

0 引 言

煤炭超臨界水氣化技術(shù)是一種新型的煤氣化技術(shù),利用超臨界水的強(qiáng)氧化性和高溶解性,將煤直接轉(zhuǎn)化為H2、CO、CO2等產(chǎn)物,同時將煤中氮、硫、重金屬等污染物以固體無機(jī)鹽的形式沉淀,是實現(xiàn)煤高效清潔低碳利用的理想技術(shù)之一[1-3]。近年來針對超臨界水煤氣化過程的試驗研究很多,主要集中在分析溫度[1-4]、壓力[5-7]、流量和催化劑等系統(tǒng)[8-11]運行參數(shù)對氣化效率的影響。研究表明,溫度是超臨界水煤氣化轉(zhuǎn)化的決定因素,在中低溫條件下超臨界水煤轉(zhuǎn)化過程中難以完全氣化,導(dǎo)致反應(yīng)器堵塞[12]。雖然提高反應(yīng)溫度能保證煤完全轉(zhuǎn)化,但會嚴(yán)重腐蝕反應(yīng)器并造成大量不可逆熱損失。實現(xiàn)中低溫條件下煤的完全轉(zhuǎn)化仍存在巨大挑戰(zhàn)。超臨界水煤氣化過程十分復(fù)雜,涉及水解、熱解、脫水、加氫、聚合和異構(gòu)化等一系列反應(yīng)[13],同時由于超臨界水(溫度大于374 ℃、壓力高于22.1 MPa)的特殊環(huán)境,限制了試驗裝置在線測量,導(dǎo)致超臨界水煤氣化過程的動力學(xué)機(jī)制尚未明晰。

多環(huán)芳烴的開環(huán)反應(yīng)被認(rèn)為是超臨界水煤氣化過程關(guān)鍵速控步反應(yīng),目前文獻(xiàn)中普遍認(rèn)為超臨界水是芳香環(huán)開環(huán)的促進(jìn)劑[13],但具體反應(yīng)機(jī)制尚不清楚,對此主要有2種不同解釋。JIN等[14]采用反應(yīng)分子動力學(xué)方法結(jié)合密度泛函理論對煤替代物蒽在超臨界水系統(tǒng)中的氣化過程進(jìn)行模擬,比較了蒽直接熱解、亞臨界水蒸氣氛圍及超臨界水氛圍中氣化過程的差異。研究發(fā)現(xiàn),超臨界水系統(tǒng)中氫產(chǎn)率遠(yuǎn)高于熱解和水蒸氣氣化系統(tǒng),這與試驗得出的規(guī)律一致。軌跡分析顯示,超臨界水通過與萘發(fā)生加成反應(yīng)形成復(fù)合物,使C—C鍵裂解能由776.7 kJ/mol 降至218.5 kJ/mol,加速了芳香族化合物的開環(huán)過程。HAN等[15]采用反應(yīng)分子動力學(xué)方法對煤替代物萘在超臨界水系統(tǒng)中的氣化過程進(jìn)行模擬,研究指出,超臨界水不僅提供了生成H2的氫原子,同時還提供了生成CO中的氧原子。在開環(huán)過程中,超臨界水首先通過提供OH自由基或直接與萘發(fā)生脫氫反應(yīng)生成萘自由基,隨后通過分子碰撞促進(jìn)開環(huán)過程。

由于模擬結(jié)果包含上萬種組分信息和基元反應(yīng),缺乏反應(yīng)動力學(xué)層次的分析技術(shù),導(dǎo)致目前主要通過人工發(fā)現(xiàn)關(guān)鍵反應(yīng),且僅能獲得定性水平的結(jié)論,使超臨界水煤氣化反應(yīng)的關(guān)鍵速控步機(jī)理仍存在爭議。在現(xiàn)有分析程序中,美國桑迪亞國家實驗室開發(fā)的開源LAMMPS程序[16]中REAXC模塊[17]能輸出模擬過程中系統(tǒng)化學(xué)鍵和物種變化,但難以提供任何化學(xué)反應(yīng)層面的信息。D?NTGEN等[18-19]基于Python提出了Chemtrayzer開源程序,該程序通過分析反應(yīng)過程的化學(xué)鍵變化給出了詳細(xì)的化學(xué)反應(yīng)路徑,但計算效率較低,僅適用于小分子體系。中國科學(xué)院過程工程研究所李曉霞等針對煤體系提出了VARxMD程序[20],該程序在對3D化學(xué)結(jié)構(gòu)的唯一物種識別、反應(yīng)位點識別、化學(xué)鍵識別的基礎(chǔ)上,能對模擬結(jié)果進(jìn)行可視化分析,并生成完整的化學(xué)反應(yīng)列表,該軟件非開源,代碼并未公開。華東師范大學(xué)朱通等開發(fā)了ReacNetGenerator開源程序[21],該程序可從反應(yīng)軌跡中自動提取反應(yīng)網(wǎng)絡(luò),并成功預(yù)測了RP-3替代模型燃燒過程的反應(yīng)軌跡。

因此,筆者構(gòu)建針對反應(yīng)分子動力學(xué)模擬結(jié)果的化學(xué)反應(yīng)路徑分析方法,用以分析超臨界水氛圍多環(huán)芳烴詳細(xì)反應(yīng)過程,揭示關(guān)鍵自由基的形成和消耗途徑,并確定開環(huán)反應(yīng)的關(guān)鍵反應(yīng)路徑;利用定量化評估方法對比不同工況結(jié)果,實現(xiàn)分析結(jié)論的低溫外推,為超臨界水多環(huán)芳烴的動力學(xué)機(jī)理建模和模型優(yōu)化提供依據(jù)。

1 計算方法

1.1 模擬方法

分子動力學(xué)方法結(jié)合反應(yīng)力場是多相體系化學(xué)反應(yīng)模擬的主要手段,賓夕法尼亞大學(xué)VAN DUIN等[22]提出了反應(yīng)力場方法,該方法采用鍵級描述系統(tǒng)中原子間化學(xué)鍵強(qiáng)弱,鍵級由原子之間的距離和原子對應(yīng)的化學(xué)價決定。鍵級能動態(tài)描述系統(tǒng)鍵能、鍵角能和二面角能的變化,因此可表征系統(tǒng)中化學(xué)鍵的斷裂和生成。此外,反應(yīng)力場用電負(fù)性平衡方法處理極化作用,考慮了范德華和庫倫相互作用,參數(shù)化擬合后的反應(yīng)力場能較好重現(xiàn)密度泛函方法精度下的反應(yīng)勢能面。相比高計算成本的第一性原理模擬方法,反應(yīng)力場方法計算效率可提升3個數(shù)量級以上[22],實現(xiàn)數(shù)十萬原子在納米尺度下和納秒時間內(nèi)模擬,并在原子尺度下獲得系統(tǒng)反應(yīng)過程中的溫度分布、壓力分布、物種濃度和化學(xué)反應(yīng)速率等基本參數(shù),對吸附、擴(kuò)散和相變等多相反應(yīng)過程進(jìn)行準(zhǔn)確描述[23]。分子動力學(xué)結(jié)合反應(yīng)力場的方法是一種原子間反應(yīng)勢驅(qū)動的方法,僅需構(gòu)建反應(yīng)力場信息而無需預(yù)設(shè)反應(yīng)路徑,這是其應(yīng)用于復(fù)雜體系微觀化學(xué)反應(yīng)機(jī)制研究的重要優(yōu)勢。選取ASHRAF等[24]提出的CHO-2016力場勢函數(shù)參數(shù),相比2008年的CHO-2008參數(shù),顯著改進(jìn)了小分子烴類和水的熱力學(xué)和動力學(xué)描述精度,是目前針對烴類分子動力學(xué)模擬最準(zhǔn)確的力場參數(shù)集。各研究案例的詳細(xì)參數(shù)見表1,水的化學(xué)計量系數(shù)通過假定萘與水反應(yīng)完全生成H2和CO2來確定,數(shù)值為20。

表1 各模擬案例的詳細(xì)參數(shù)Table 1 Parameters for each simulation case

對比了400與100個萘分子系統(tǒng)的凈反應(yīng)發(fā)生次數(shù)統(tǒng)計(表2),可知2個系統(tǒng)差異較小?？紤]計算耗時,選取100個萘分子作為模擬參數(shù),并對所有模擬案例進(jìn)行4次模擬取平均,以減少系綜誤差。利用Packmol軟件[25]將萘分子和不同數(shù)目的水分子在空間中隨機(jī)布置完成初始模型構(gòu)建,周期性邊界條件立方體盒大小為10.6 nm×10.6 nm×10.6 nm,保證所有模擬工況中水分壓大于水的臨界壓力22.1 MPa,保證水處于超臨界狀態(tài)。模擬系統(tǒng)中固定萘濃度不變,通過調(diào)整水分子數(shù)實現(xiàn)當(dāng)量比變化。

表2 不同分子系統(tǒng)部分凈反應(yīng)發(fā)生次數(shù)比較Table 2 Comparison of net reactions times in different systems

受限于計算機(jī)性能,反應(yīng)分子動力學(xué)模擬的最大時間為納秒數(shù)量級,為保證研究系統(tǒng)充分反應(yīng),首先在2 000～3 000 K進(jìn)行模擬,通過追蹤萘分子濃度,發(fā)現(xiàn)2 500 K下萘分子在8 ns內(nèi)完全消耗,因此,選取2 500 K作為反應(yīng)的最低模擬溫度。此外,在2 600和2 700 K下進(jìn)行分子動力學(xué)模擬,通過對比不同溫度的模擬結(jié)果,對低溫工況進(jìn)行外推。選取0.1 fs作為時間步長避免數(shù)值誤差累積,選用恒定數(shù)目、恒定體積和恒定溫度的NVT系綜進(jìn)行模擬,溫度阻尼為時間步長的200倍即20 fs。模擬過程中,將系統(tǒng)進(jìn)行能量最小化優(yōu)化,消除初始結(jié)構(gòu)中不合理的原子位置,再對所有原子分配滿足溫度為1 000 K的高斯分布初始速度,并在1 000 K下進(jìn)行10 ps初始平衡,200 ps內(nèi)將系統(tǒng)溫度升至目標(biāo)模擬溫度,在目標(biāo)溫度下固定溫度進(jìn)行模擬并輸出系統(tǒng)模擬結(jié)果,結(jié)果分析中的零時刻均指升溫完成時刻。所有模擬均通過開源分子動力學(xué)模擬軟件LAMMPS[16]完成。

1.2 分析方法

鍵級是描述模擬系統(tǒng)反應(yīng)過程化學(xué)信息的重要參數(shù),不僅用于識別物種的唯一結(jié)構(gòu),還可用于描述基元反應(yīng)變化。在Chemtrayzer基礎(chǔ)上,結(jié)合圖論分析庫重現(xiàn),從模擬系統(tǒng)中的鍵級進(jìn)行物質(zhì)識別和基元反應(yīng)識別,并對基元反應(yīng)進(jìn)行統(tǒng)計獲取反應(yīng)路徑,同時利用并行計算技術(shù)提高分析速度,運行流程如圖1所示。程序首先讀取LAMMPS輸出的鍵級數(shù)據(jù),并將鍵級進(jìn)行二值化處理,即鍵級大于閾值0.3時則認(rèn)為2個原子間存在化學(xué)鍵。需要說明的是,本研究僅關(guān)注芳香環(huán)斷鍵過程,不考慮化學(xué)鍵是否為單鍵、雙鍵等鍵類型信息,忽略鍵類型信息可減少同分異構(gòu)體帶來的反應(yīng)路徑噪聲,同時可不丟失芳香環(huán)的結(jié)構(gòu)信息。根據(jù)鍵級通過拓?fù)鋵W(xué)理論進(jìn)行連接性搜索,尋找與某原子存在連接關(guān)系的所有原子,并調(diào)用Open Babel程序庫[26]將連接關(guān)系與各原子對應(yīng)元素類型進(jìn)行識別,轉(zhuǎn)換成SMILES字符串[27-28],同時利用Open Babel程序SSSR模塊分析結(jié)構(gòu)的環(huán)特性,實現(xiàn)某一時刻下系統(tǒng)化學(xué)信息分析。從時間維度分析,如果不同時刻下某一原子的化學(xué)鍵信息發(fā)生變化,則根據(jù)化學(xué)鍵的生成或斷裂情況和對應(yīng)發(fā)生反應(yīng)前后時刻的物種,生成對應(yīng)的基元反應(yīng)化學(xué)式。通過區(qū)分和合并正、逆化學(xué)反應(yīng)式處理反應(yīng)過程中大量逆向反應(yīng),獲得整個反應(yīng)過程各基元反應(yīng)累積發(fā)生的次數(shù)統(tǒng)計,在此基礎(chǔ)上進(jìn)行反應(yīng)動力學(xué)水平分析,獲取從反應(yīng)物到穩(wěn)定產(chǎn)物的詳細(xì)反應(yīng)路徑。此外,原始反應(yīng)路徑中存在大量僅活性位點不同的同分異構(gòu)體,使反應(yīng)路徑稀疏,難以提取有用信息。為發(fā)現(xiàn)關(guān)鍵反應(yīng)通道,分析過程中合并了分子式相同和碳環(huán)結(jié)構(gòu)相同的物種。合并簡化反應(yīng)路徑,如對于反應(yīng)A+B<->C,當(dāng)C物種的生成路徑僅存在C<->D+E時,2個反應(yīng)將合并為A+B<->D+E?？紤]到模擬溫度達(dá)2 700 K,高溫下分子碰撞頻繁,為保證精確捕捉單個化學(xué)鍵斷裂的瞬間,使識別反應(yīng)為僅涉及一個化學(xué)鍵斷裂的基元反應(yīng),分析步長選取了比Chemtrayzer[18-19]建議值20 fs更低的1 fs分析步長,對應(yīng)為模擬時間步長的10倍,通過合并正、逆化學(xué)反應(yīng)式處理小步長帶來的震蕩反應(yīng)。識別鍵級文件中時間戳信息,將鍵級文件分塊,利用多核并行處理不同時刻的鍵級狀態(tài),提升分析效率,后處理耗時相較LAMMPS模擬耗時短1個數(shù)量級;將鍵級文件進(jìn)行實時壓縮,使磁盤占用減少一個數(shù)量級,解決更小分析步長導(dǎo)致的更多分析耗時和磁盤占用等問題。分析程序通過Python語言實現(xiàn)。

圖1 后處理程序的分析流程Fig.1 Flowchart of the post-processing program

2 模擬結(jié)果分析

2.1 邊界條件對反應(yīng)過程的影響

2 500 K下當(dāng)量比對反應(yīng)物萘和中間產(chǎn)物CO的影響如圖2所示(N0為初始值,φ為當(dāng)量比,λ為衰減常數(shù),α為擬合系數(shù)),模擬過程中固定萘分子濃度,通過控制水分子濃度調(diào)節(jié)當(dāng)量比?？芍?dāng)量比變化對萘分子消耗速率影響較小,說明萘分子消耗與水分子濃度呈弱相關(guān),高溫條件熱解反應(yīng)主導(dǎo)反應(yīng)進(jìn)程。水分子濃度對CO生成速率影響顯著,由于反應(yīng)系統(tǒng)中氧原子均來自水分子,說明水分子參與了CO生成過程,控制基元反應(yīng)需進(jìn)一步動力學(xué)分析。

圖2 2 500 K不同當(dāng)量比萘和CO分子數(shù)隨時間變化Fig.2 Molecular numbers of naphthalene and carbon monoxide with different equivalence ratios at 2 500 K

溫度是影響反應(yīng)快慢和反應(yīng)路徑變化的重要參數(shù),因此,對比了溫度對反應(yīng)物濃度的影響如圖3所示?？芍磻?yīng)物呈指數(shù)下降,且下降速率隨溫度升高而增加。曲線呈指數(shù)函數(shù)衰減,衰減常數(shù)λ與溫度呈線性相關(guān),衰減常數(shù)越大意味著反應(yīng)物氧化耗時越長。反應(yīng)物消耗95%時,模擬時間分別為4.58、3.02和1.58 ns,說明更低溫度的反應(yīng)分子動力學(xué)模擬耗時呈指數(shù)增長趨勢。反應(yīng)物濃度與指數(shù)函數(shù)較好的吻合行為說明當(dāng)前模擬條件下反應(yīng)過程速控步為一階反應(yīng),該反應(yīng)過程與溫度強(qiáng)相關(guān)。

2.2 反應(yīng)路徑分析

為分析反應(yīng)物萘的消耗路徑,闡明萘開環(huán)反應(yīng)路徑主導(dǎo)的基元反應(yīng),對2 600 K、當(dāng)量比為1條件,萘消耗95%時的反應(yīng)路徑進(jìn)行分析,如圖4所示,圖中合并了同分異構(gòu)體,不區(qū)分反應(yīng)位點?？芍?① 萘的消耗路徑包括熱解、脫氫、加成反應(yīng),其中熱解是主要消耗路徑,占比超過60%,熱解反應(yīng)的高占比與模擬溫度直接相關(guān);② 萘外部碳碳鍵比內(nèi)部碳碳鍵更易斷裂,萘外部碳碳鍵斷裂消耗占比為內(nèi)部斷裂占比的5倍以上,但萘脫去一個氫原子后內(nèi)部碳碳鍵斷裂消耗占比與外部相當(dāng);③ 對于開環(huán)反應(yīng),即使在2 600 K下,萘分子脫氫反應(yīng)仍為開環(huán)的關(guān)鍵步驟,占比超60%;④ 水分子難以直接與萘分子發(fā)生加成反應(yīng)形成化學(xué)鍵,相比之下,萘分子脫氫后,水分子中的氧原子加成到萘原來的脫氫位點更容易,此后水分子再消除一個氫原子形成醇類,最后發(fā)生開環(huán)反應(yīng),這與JIN等[14]的結(jié)論不同。盡管本文后處理分析中捕捉到了萘分子與水分子的加成基元反應(yīng),但考慮其逆反應(yīng)過程后,凈反應(yīng)發(fā)生次數(shù)在反應(yīng)路徑中占比不超過0.5%,故為非主要反應(yīng)路徑。

圖4 2 600 K下當(dāng)量比為1時萘消耗95%時的萘開環(huán)反應(yīng)路徑Fig.4 Ring opening reaction pathway diagram at 2 600 K with 95% naphthalene consumption at an equivalence ratio of 1

不同溫度下直接消耗萘分子反應(yīng)的發(fā)生次數(shù)見表3。由于反應(yīng)過程中存在循環(huán)反應(yīng)過程,且不同溫度下模擬進(jìn)行的時間不同,使凈反應(yīng)發(fā)生次數(shù)絕對值對比規(guī)律性不強(qiáng)。

表3 當(dāng)量比為1時不同溫度下萘消耗95%時直接消耗萘分子凈反應(yīng)發(fā)生次數(shù)Table 3 Net reactions times that directly consume naphthalene with 95% naphthalene consumption at an equivalence ratio of 1

將數(shù)據(jù)進(jìn)行歸一化處理,根據(jù)反應(yīng)類型定義反應(yīng)1～3為熱解主導(dǎo)反應(yīng)類,反應(yīng)4～6為自由基主導(dǎo)反應(yīng)類。自由基反應(yīng)類在萘分子總消耗比值與溫度的關(guān)系如圖5所示?？芍S溫度降低,自由基反應(yīng)類占比不斷提升,進(jìn)一步通過線性擬合將結(jié)果進(jìn)行低溫外推發(fā)現(xiàn),溫度低于2 000 K時,自由基反應(yīng)類占比超過50%,溫度低于1 000 K時自由基反應(yīng)類占比將超過75%。盡管實際超臨界水煤反應(yīng)過程中自由基相關(guān)反應(yīng)主導(dǎo)反應(yīng)過程,但該結(jié)論仍需試驗和相關(guān)理論計算進(jìn)一步驗證。

2.3 物種通量分析

圖6 當(dāng)量比為1時不同溫度下生成OH自由基的基元反應(yīng)Fig.6 Main elementary reaction of OH radical formation at an equivalence ratio of 1 with different temperatures

圖7 當(dāng)量比為1、2 500 K下萘消耗95%時生成氫原子的主要基元反應(yīng)Fig.7 Main elementary reaction of generating hydrogen atom when naphthalene consumption is 95% at an equivalence ratio of 1,2 500 K

3 結(jié) 論

1)采用反應(yīng)分子動力學(xué)方法研究了不同溫度和不同當(dāng)量比條件下超臨界水和萘的氣化過程。計算結(jié)果表明,高溫條件下萘分子數(shù)呈指數(shù)下降,消耗速率與當(dāng)量比呈弱相關(guān),呈現(xiàn)明顯的一階反應(yīng)特性。通過分析系統(tǒng)每一時刻的鍵級變化,獲得了定量化的基元反應(yīng)路徑,萘的開環(huán)反應(yīng)路徑包括直接熱解、脫氫后開環(huán)和OH加成開環(huán)。

2)高溫條件下直接熱解占主導(dǎo),1 000 K下脫氫和OH加成反應(yīng)占主導(dǎo)。本研究與文獻(xiàn)結(jié)論并不一致,超臨界水直接提供活性H原子和OH自由基以及超臨界水直接與反應(yīng)物發(fā)生反應(yīng)的占比極小,非主導(dǎo)反應(yīng)。