正丙苯高溫氧化機(jī)理的分子動(dòng)力學(xué)模擬研究

周子豪， 王思皓， 黃玳川， 劉波， 甯紅波

正丙苯高溫氧化機(jī)理的分子動(dòng)力學(xué)模擬研究

周子豪， 王思皓， 黃玳川， 劉波， 甯紅波

（西南交通大學(xué)材料先進(jìn)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室, 成都 610031）

正丙苯是Jet A、 Jet A-1及國(guó)產(chǎn)RP-3航空煤油中芳香烴的典型替代組分. 本文采用基于反應(yīng)力場(chǎng)的分子動(dòng)力學(xué)模擬研究了正丙苯高溫氧化過(guò)程的主要反應(yīng)網(wǎng)絡(luò)、 主要產(chǎn)物的形成機(jī)理以及在不同溫度、 密度和當(dāng)量比條件下正丙苯氧化主要產(chǎn)物的分布規(guī)律, 并結(jié)合反應(yīng)動(dòng)力學(xué)理論計(jì)算了正丙苯高溫氧化的速率常數(shù). 結(jié)果表明， 正丙苯高溫氧化主要發(fā)生在烷基側(cè)鏈, 包括6種C—C和C—H鍵斷裂單分子分解反應(yīng)以及3種側(cè)鏈氫原子與氧氣或其它小自由基的氫提取反應(yīng), 其中與芐基相連的C—C鍵具有最小斷鍵能, 是最重要的單分子分解反應(yīng), 而不同位點(diǎn)的自由基氫提取反應(yīng)的貢獻(xiàn)相似; 體系的模擬溫度和密度/壓力與正丙苯的氧化速率呈正相關(guān), 但當(dāng)量比對(duì)氧化速率的影響則嚴(yán)重依賴于體系溫度. 計(jì)算所得正丙苯高溫氧化表觀活化能和指前因子與文獻(xiàn)報(bào)道的實(shí)驗(yàn)值相比在可接受的范圍內(nèi).

正丙苯； 反應(yīng)機(jī)理； 高溫氧化； 反應(yīng)力場(chǎng)； 分子動(dòng)力學(xué)模擬

航空煤油作為航空發(fā)動(dòng)機(jī)燃燒的主要能源， 如何提高其燃燒效率、 降低污染物的生成排放越發(fā)受到關(guān)注. 由于航空煤油包含多種烷烴、 烯烴、 環(huán)烷烴和芳香烴等碳?xì)浣M分， 很難直接構(gòu)建其燃燒動(dòng)力學(xué)模型. 通常采用能夠表征航空煤油重要理化特性的替代燃料模型研究航空煤油的燃燒特性［1］. 在航空煤油的替代燃料模型研究中， 正丙苯作為芳香烴組分的典型替代燃料， 它具有足夠長(zhǎng)的烷基側(cè)鏈（正丙基，-C3H7）和苯基（—C6H5）， 可以研究側(cè)鏈脂肪族性質(zhì)和正丙苯的芳香族性質(zhì). 如， Dagaut等［2］構(gòu)建的航空煤油替代燃料模型中正丙苯占26%， 徐佳琪等［3］開(kāi)發(fā)的RP-3三組分替代燃料模型中正丙苯占15.8%， 而Dooley等［4］開(kāi)發(fā)的Jet A燃料和第二代替代燃料模型中正丙苯占22.8%（摩爾分?jǐn)?shù)）. 此外， 與其它單環(huán)芳烴（如甲苯、 乙基苯和二甲苯）相比， 正丙苯氧化可以產(chǎn)生更多的小分子烴類（如甲烷、 乙烯和丙烯）， 與丁苯及其它更大的烷基單取代苯相比， 其側(cè)鏈較短能夠有效減少動(dòng)力學(xué)模型的構(gòu)建規(guī)模. 因此， 為了深入理解正丙苯的燃燒特性， 研究人員對(duì)正丙苯展開(kāi)了系列的實(shí)驗(yàn)、 動(dòng)力學(xué)模擬和理論計(jì)算研究.

在實(shí)驗(yàn)方面， 研究主要關(guān)注正丙苯的層流燃燒速率［5～9］和點(diǎn)火延遲時(shí)間［10～13］. 這兩類實(shí)驗(yàn)數(shù)據(jù)能夠直觀反映正丙苯燃燒過(guò)程的鏈誘發(fā)和鏈傳遞動(dòng)力學(xué)特征， 對(duì)正丙苯燃燒動(dòng)力學(xué)模型的開(kāi)發(fā)和驗(yàn)證至關(guān)重要. 此外， 正丙苯燃燒熱解過(guò)程的主要物種濃度分布為正丙苯燃燒動(dòng)力學(xué)模型中關(guān)鍵反應(yīng)網(wǎng)絡(luò)及熱動(dòng)力學(xué)參數(shù)的合理性驗(yàn)證提供了更為豐富的實(shí)驗(yàn)依據(jù)［14～17］. 在動(dòng)力學(xué)模擬方面， 對(duì)于正丙苯的動(dòng)力學(xué)模型構(gòu)建主要根據(jù)反應(yīng)類規(guī)則， 以甲苯及乙苯的動(dòng)力學(xué)模型為基礎(chǔ)， 建立正丙苯的燃燒動(dòng)力學(xué)模型［17，18］， 在此基礎(chǔ)上添加包括烷基過(guò)氧自由基和氫過(guò)氧烷基自由基的相關(guān)反應(yīng)， 建立正丙苯的低溫氧化動(dòng)力學(xué)模型［12，19］， 此外， 通過(guò)刪除正丙苯燃燒模型中的含氧物種及相關(guān)反應(yīng)， 建立正丙苯的熱解動(dòng)力學(xué)模型［14，15］. 在理論計(jì)算方面， 主要采取量子化學(xué)計(jì)算， 針對(duì)正丙苯燃燒裂解過(guò)程中最關(guān)鍵的初始消耗反應(yīng)［即小自由基（H/CH3/HO2/O2）］與正丙苯的提取反應(yīng)動(dòng)力學(xué)開(kāi)展研究［14，20～22］.

隨著先進(jìn)發(fā)動(dòng)機(jī)工況向極端條件（如超高溫、 超高壓）發(fā)展， 受限于反應(yīng)器本身材料的性質(zhì)、 使用環(huán)境的特殊性以及燃燒反應(yīng)的快速性， 導(dǎo)致現(xiàn)有的燃燒實(shí)驗(yàn)診斷技術(shù)難以直接監(jiān)測(cè)燃料反應(yīng)過(guò)程中多個(gè)特定產(chǎn)物的時(shí)間演變（fs級(jí)）. 雖然前述實(shí)驗(yàn)提供了正丙苯氧化的初始產(chǎn)物和有關(guān)反應(yīng)的必要信息， 然而， 特別是在高溫下， 氧化反應(yīng)進(jìn)展非常快， 并產(chǎn)生大量的中間體和活性自由基， 難以通過(guò)實(shí)驗(yàn)捕捉. 此外， 由于反應(yīng)的時(shí)間尺度和空間尺度都很小， 從檢測(cè)到的有限物質(zhì)中推斷出詳細(xì)的氧化反應(yīng)過(guò)程也不夠全面. 因此， 采用理論方法研究極端條件下的燃料燃燒反應(yīng)特性是一種較好的手段. 原則上， 量子化學(xué)計(jì)算可以提供準(zhǔn)確的反應(yīng)物、 過(guò)渡態(tài)和產(chǎn)物的結(jié)構(gòu)信息， 但主要針對(duì)小分子燃料體系［23］. 對(duì)于復(fù)雜體系特別是大分子燃料體系， 量子化學(xué)計(jì)算所需的時(shí)間成本極其昂貴而受到限制.

由van Duin等［24］基于鍵級(jí)和鍵能關(guān)系開(kāi)發(fā)的ReaxFF反應(yīng)力場(chǎng)， 其力場(chǎng)參數(shù)完全來(lái)自對(duì)力場(chǎng)目標(biāo)體系的量子化學(xué)計(jì)算結(jié)果的擬合， 可以快速準(zhǔn)確地觀察到化學(xué)過(guò)程中反應(yīng)物、 中間體和產(chǎn)物的時(shí)間演化趨勢(shì)， 為研究較大反應(yīng)體系的物理化學(xué)性質(zhì)提供了有效途徑. 并且經(jīng)過(guò)多年的發(fā)展， ReaxFF反應(yīng)力場(chǎng)已被廣泛應(yīng)用于含能材料、 催化、 燃燒和聚合物等多個(gè)研究領(lǐng)域［25］. 為了更加準(zhǔn)確地描述碳?xì)淙剂系娜紵^(guò)程， Chenoweth等［26］基于量子化學(xué)計(jì)算結(jié)果優(yōu)化了ReaxFF反應(yīng)力場(chǎng)的力場(chǎng)參數(shù)， 并對(duì)一系列的碳?xì)淙剂先紵匦赃M(jìn)行模擬， 結(jié)果表明， ReaxFF反應(yīng)力場(chǎng)可以正確預(yù)測(cè)碳?xì)淙剂系姆磻?yīng)活性趨勢(shì)［27］. 此外， 其他一些研究者也使用該力場(chǎng)參數(shù)研究了正十二烷、 甲苯、 RP-3和乙烯的燃燒特性， 其模擬結(jié)果與密度泛函理論方法計(jì)算結(jié)果和可測(cè)量的實(shí)驗(yàn)結(jié)果也相吻合［28～31］. 綜上表明， ReaxFF反應(yīng)力場(chǎng)適合于碳?xì)淙剂系娜紵匦匝芯?

本文采用ReaxFF反應(yīng)力場(chǎng)對(duì)正丙苯高溫氧化的起始反應(yīng)路徑、 主要產(chǎn)物的形成機(jī)理以及不同溫度、 密度和當(dāng)量比下正丙苯高溫氧化的主要產(chǎn)物分布規(guī)律進(jìn)行了分子動(dòng)力學(xué)模擬研究， 并采用反應(yīng)動(dòng)力學(xué)方法計(jì)算了正丙苯的高溫氧化動(dòng)力學(xué)參數(shù).

1 計(jì)算方法與模擬細(xì)節(jié)

1.1　ReaxFF反應(yīng)力場(chǎng)

ReaxFF反應(yīng)力場(chǎng)是由van Duin等［24］基于第一性原理提出的， 以原子間的鍵級(jí)為核心， 通過(guò)鍵級(jí)、 鍵長(zhǎng)和鍵能之間的關(guān)系以及原子間電荷的動(dòng)態(tài)傳遞來(lái)模擬原子之間化學(xué)鍵的生成和斷裂. 在分子動(dòng)力學(xué)的計(jì)算過(guò)程中， 當(dāng)化學(xué)鍵斷裂時(shí)， 與價(jià)鍵相關(guān)的力和能量變?yōu)榱悖?因此該反應(yīng)力場(chǎng)能夠描述化學(xué)鍵的生成和斷裂. ReaxFF反應(yīng)力場(chǎng)的計(jì)算表達(dá)式如下：

上式主要由兩大部分組成： 基于鍵級(jí)的成鍵相互作用能量項(xiàng)（kJ/mol）： 鍵能項(xiàng)（bond）、 配位修正項(xiàng) （over和under）、 鍵角能項(xiàng)（val）、 懲罰能項(xiàng)（pen）、 三體共軛項(xiàng)（coa）、 二面角共軛作用能項(xiàng)（tors）和四面體共軛項(xiàng)（conj）等； 非鍵相互作用能項(xiàng)（kJ/mol）： 范德華相互作用能項(xiàng)（vdW）和庫(kù)侖力相互作用能項(xiàng)（Coulomb）.

初始分子結(jié)構(gòu)采用密度泛函理論（DFT）方法B3LYP泛函結(jié)合6-311++G（，）基組進(jìn)行優(yōu)化. 在ReaxFF模擬中， 在5 K時(shí)進(jìn)行能量最小化， 時(shí)間步長(zhǎng)為0.1 fs， 原子數(shù)（）、 體積（）和溫度（）恒定， 表示為NVT-MD（分子動(dòng)力學(xué)模擬）. 然后使用NVT-MD模擬， 將這些系統(tǒng)在5 K下馳豫， 使得模擬體系達(dá)到平衡， 平衡50 ps， 時(shí)間步長(zhǎng)為0.1 fs. 時(shí)間步長(zhǎng)的設(shè)置主要參考Chenoweth等［26］的工作. 在所有 NVT-MD模擬過(guò)程中， 采用Berendsen方法［32］對(duì)溫度和壓力進(jìn)行調(diào)節(jié)， 耦合參數(shù)為500 fs， 使溫度在設(shè)定值附近波動(dòng). 將鍵級(jí)設(shè)定為0.3， 當(dāng)鍵級(jí)大于0.3時(shí)， 判定新化學(xué)鍵形成， 產(chǎn)物生成. 當(dāng)量比（）由正丙苯實(shí)際氧化過(guò)程的正丙苯/氧氣的摩爾比與正丙苯完全燃燒所需正丙苯/氧氣的摩爾比之比確定， 體系密度通過(guò)控制立方體體積設(shè)置， 所有輸入條件的匯總見(jiàn)表1， 其范圍從在最低溫度和密度（= 2300 K，=0.05 g/cm3，=2.0）到最高溫度和密度（=3500 K，=0.35 g/cm3，=0.5）， 關(guān)于體系密度設(shè)置主要參考已報(bào)道的類似反應(yīng)體系［27，28，33］， 本文通過(guò)理想氣體狀態(tài)方程， 計(jì)算出的壓力范圍為22～286 MPa. 注意提高模擬溫度研究燃料燃燒和熱解是ReaxFF模擬的常見(jiàn)做法， 旨在通過(guò)提高溫度來(lái)縮短反應(yīng)時(shí)間， 同時(shí)不改變反應(yīng)路徑， 又保留原本溫度下的動(dòng)力學(xué)行為［34］. 本文采用的是Chenoweth等［26］報(bào)道的CHO-2008力場(chǎng)參數(shù). 目前， 該研究小組有更新改進(jìn)的CHO-2016力場(chǎng)參數(shù)［35］， 該力場(chǎng)參數(shù)與 CHO-2008相比， 保持了CHO-2008對(duì)大分子碳?xì)浠衔锏哪M精度， 主要改善了CHO-2008對(duì)C1化學(xué)的預(yù)測(cè)精度. 值得注意的是， 文獻(xiàn)［35］進(jìn)一步驗(yàn)證了使用CHO-2016與CHO-2008研究正丁基苯熱解初始反應(yīng)機(jī)理， 表明CHO-2016與CHO-2008模擬結(jié)果一致. 最近， Lindgren等［33］采用CHO-2008探究環(huán)戊烷和環(huán)戊烯的燃燒初始反應(yīng)機(jī)理， 他們也指出CHO-2016主要改善C1化學(xué)， 對(duì)于研究較大的碳?xì)浠衔铮?CHO-2008也是可以的. 因此， 選擇CHO-2008研究正丙苯的高溫氧化機(jī)理. 此外， 采用Packmol程序［36］構(gòu)建周期性體系， 并且所有的MD模擬都在LAMMPS軟件［37］中完成， 使用ChemTrayzer程序［38］進(jìn)行反應(yīng)路徑的后處理.

Table 1　Temperature, density, and equivalence ratio of n-propylbenzene oxidation system

1.2　DFT計(jì)算

為了驗(yàn)證反應(yīng)力場(chǎng)的合理性和可遷移性， 同時(shí)采用DFT計(jì)算和ReaxFF模擬研究了正丙苯氧化過(guò)程中側(cè)鏈C—C鍵的鍵解離能， 為了分析正丙苯燃燒化學(xué)動(dòng)力學(xué)過(guò)程中的初始?xì)涮崛》磻?yīng)， 采用DFT 計(jì)算了正丙苯不同碳原子位點(diǎn)上的C—H鍵的鍵能. 所有分子和自由基的幾何結(jié)構(gòu)優(yōu)化采用DFT 方法B3LYP泛函結(jié)合6-311++G（，）基組進(jìn)行. DFT計(jì)算均采用Guassian 16量子化學(xué)計(jì)算軟件包［39］完成.

2 結(jié)果與討論

2.1　正丙苯氧化的初始反應(yīng)機(jī)理

為了確定正丙苯氧化的初始反應(yīng)機(jī)理， 盡量減小由于反應(yīng)體系過(guò)小、 反應(yīng)物的數(shù)量過(guò)少帶來(lái)結(jié)果的差異過(guò)小而不能顯著觀察出初始反應(yīng)通道， 將50個(gè)正丙苯分子和600個(gè)氧氣分子放置在溫度為2900 K、 密度為0.35 g/cm3條件下進(jìn)行NVT?MD模擬. 通過(guò)分析原子成鍵情況以及在不同時(shí)刻下的氧化產(chǎn)物， 最終得到了正丙苯高溫氧化的初始反應(yīng)途徑. 正丙苯初始反應(yīng)主要可分為兩大類： 正丙苯側(cè)鏈上發(fā)生裂解反應(yīng)以及與其它自由基發(fā)生的提取反應(yīng). 圖1展示了觀察到的正丙苯發(fā)生裂解的4條反應(yīng)路徑及頻率. 正丙苯鏈上發(fā)生C—C斷裂生成芐基和乙基（CH3CH2）是最主要的裂解反應(yīng)， 反應(yīng)頻率為25次， 占到正丙苯初始消耗的一半. 其次是C—C斷裂生成C6H5CH2CH2和CH3自由基， 反應(yīng)頻率為 5次. 裂解反應(yīng)中除C—C鍵的斷裂外， C—H鍵的斷裂也比較重要， 正丙苯通過(guò)斷裂側(cè)鏈上不同碳位上的C—H鍵生成不同的自由基： C6H5CH2CH2CH2自由基和C6H5CH2CHCH3自由基， 頻率分別為2和1. 如 圖2所示， 正丙苯與不同自由基及氧氣分子發(fā)生反應(yīng)， 模擬中共發(fā)現(xiàn)了9條不同的反應(yīng)路徑， 總頻率為17次. 值得注意的是， 氫提取反應(yīng)頻率就占到了14次. 上述結(jié)果表明正丙苯消耗的初始反應(yīng)主要為裂解反應(yīng)與氫提取反應(yīng). 由于裂解反應(yīng)在正丙苯初始消耗中至關(guān)重要. 因此， 使用DFT方法和ReaxFF方法計(jì)算了正丙苯側(cè)鏈上不同C原子上C—H和C—C鍵的鍵能.

Fig.1　Four reaction pathways of n?propylbenzene pyrolysis

Green and white spheres represent carbon and hydrogen atoms， respectively.

Fig.2　Nine reaction pathways of n?propylbenzene reacting with different small radicals and oxygen

Green， red and white spheres represent carbon， oxygen and hydrogen atoms， respectively.

結(jié)果如表2［17，40］所示， 與芐基相連的C—C鍵的鍵能最小. 因此在反應(yīng)過(guò)程中與芐基相連的C—C鍵最易斷裂， 這與上述模擬結(jié)果相符. 碳?xì)浠衔锶紵牡谝徊绞峭ㄟ^(guò)打破單個(gè)C—C鍵來(lái)分解燃料分子［41］， C—C鍵通常比C—H弱［40］. 此外， 可以看出， 不論是B3LYP還是ReaxFF， 與Wang等［17］采用高精度CBS-QB3方法計(jì)算的正丙苯單分子分解主要斷鍵能和Luo等［40］的斷鍵能數(shù)據(jù)庫(kù)相比， 結(jié)果均吻合較好， 說(shuō)明采用ReaxFF研究正丙苯高溫氧化反應(yīng)體系是合適的.

Table 2　Calculated bond dissociation energies of n-propylbenzene using B3LYP and ReaxFF methods and the corresponding literature results

通過(guò)ChemTraYzer程序［38］對(duì)生成的鍵級(jí)信息進(jìn)行后處理， 獲得特定目標(biāo)組分的每一步反應(yīng)， 再將反應(yīng)物和生成物串聯(lián)即得到正丙苯的主要反應(yīng)網(wǎng)絡(luò)， 如Scheme 1所示， 它提供了從最初反應(yīng)到最終產(chǎn)物反應(yīng)機(jī)理的詳細(xì)描述. 結(jié)果顯示， 無(wú)論是正丙苯（PBZ）直接發(fā)生裂解還是與其它自由基反應(yīng)形成苯丙基自由基（PBZJA）， 隨后大部分都會(huì)生成芐基（PHCH2）， 這與實(shí)驗(yàn)結(jié)果一致［42］， 說(shuō)明PHCH2是非常重要的中間體. 之后PHCH2與HO2自由基結(jié)合生成C7H8O2分子（PHCH3O2）， 然后O—O鍵迅速斷裂生成苯甲醇自由基（PHCH2O）和OH自由基， 或直接與O自由基結(jié)合生成PHCH2O. PHCH2O的C—C鍵斷裂生成苯基（PHJA）和甲醛， PHJA氧化生成C6H5O2自由基（PHO2）， 并且進(jìn)一步氧化生成C6H5O自由基（PHO）. 接下來(lái)PHO異構(gòu)化形成五元環(huán)自由基（CPDCO）， 這是苯環(huán)消耗的重要反應(yīng)之一［1］. 之后裂解形成環(huán)戊二烯自由基（CPDJA）和CO. 最后， CPDJA發(fā)生開(kāi)環(huán)反應(yīng)， 并且進(jìn)一步氧化生成CO， CO2和H2O.

Scheme 1?Propylbenzene partial oxidation pathway at 2900 K

2.2　主要氧化產(chǎn)物的形成機(jī)理

盡管已經(jīng)對(duì)正丙苯的氧化進(jìn)行了廣泛研究， 但對(duì)中間機(jī)制的直接描述仍然十分困難. 由于對(duì)主要產(chǎn)物相關(guān)反應(yīng)的詳細(xì)分析對(duì)于理解中間機(jī)理尤其重要， 因此， 對(duì)正丙苯氧化過(guò)程中H2O， CH2O， CO和CO2的反應(yīng)途徑進(jìn)行了分析. 反應(yīng)的體系為3500 K， 0.35 g/cm3，=2.0，=330. 圖3介紹了正丙苯氧化過(guò)程中H2O， CH2O， CO和CO2的反應(yīng)途徑. 在圖3（A）中， H2O主要是通過(guò)體系中的H和OH自由基與其它分子反應(yīng)形成， 最主要的反應(yīng)通道是H2+OH→H2O+H. 此外， 含氧自由基直接裂解生成H2O也是重要的反應(yīng)通道（如C2H3O→H2O+C2HO）. 在模擬過(guò)程中，H2O主要是通過(guò)被其它分子提取H或OH而被消耗， 或H2O的直接分解（H2O→H+OH）.

Fig.3　Formation and consumption reaction pathways of the main oxidation products H2O(A), CH2O(B), CO(C) and CO2(D)

對(duì)于碳?xì)浠衔锏幕瘜W(xué)動(dòng)力學(xué)模擬研究， 甲醛（CH2O）被證實(shí)為碳?xì)浠衔锔邷厝紵闹匾虚g組分， 直接影響到最后產(chǎn)物CO2和CO等的形成［28］. 在圖3（B）中， CH2O主要是由甲氧基（CH3O）和過(guò)氧甲基（CH3O2）的裂解形成（CH3O→CH2O+H， CH3O2→CH2O+OH）， 而CH2O主要是與H自由基反應(yīng)被消耗（CH2O+H→CH3O， CH2O+H→CHO+H2）.

對(duì)于碳?xì)浠衔锏难趸?CO和CO2是主要的含碳穩(wěn)定產(chǎn)物， 該反應(yīng)在整個(gè)氧化機(jī)理中起著重要作用. 由圖3（C）可見(jiàn)， CO的形成主要是通過(guò)體系中的含氧自由基裂解形成（如CHO→CO+H， CHO2→CO+HO， C2H3O→CH3+CO）. CO的消耗主要基于CO與其它自由基復(fù)合形成含氧自由基 （如CO+H→CHO， CO+HO→CHO）. 由圖3（D）可見(jiàn)， 相比于CO， CO2的形成與CHO2更加密切， 主要反應(yīng)方式為CHO2直接發(fā)生裂解（CHO2→CO2+H）以及CHO2與各類自由基及分子發(fā)生氫提取反應(yīng)形成（如C2H+CHO2→C2H2+CO2， C2H2O+CHO2→CO2+C2H3O）， CO2主要是與H自由基反應(yīng)而被消耗（CO2+H→CHO2， CO2+H→CO+HO）.

2.3　溫度對(duì)正丙苯氧化的影響

為了分析溫度對(duì)正丙苯氧化的影響， 在不同條件（2300～3500 K、 間隔300 K，=0.35 g/cm3，=2.0，=330）下進(jìn)行了一系列NVT-MD模擬， 模擬時(shí)長(zhǎng)1 ns. 總分子數(shù)、 正丙苯分子數(shù)、 主要產(chǎn)物以及勢(shì)能隨時(shí)間的變化如圖4所示. 可見(jiàn)， 勢(shì)能衰減率隨溫度的升高下降得更快［圖4（A）］， 表明溫度加速了正丙苯的氧化. 這種現(xiàn)象的原因可以用系統(tǒng)的熱釋放率來(lái)解釋， 取決于溫度［25］. 在模擬初始階段， 尤其是在高溫條件下， 勢(shì)能會(huì)先升高然后迅速降低， 結(jié)合總分子數(shù)和正丙苯分子數(shù)變化情況 ［圖4（B）和（C）］， 可以確定整個(gè)系統(tǒng)開(kāi)始迅速吸收熱量發(fā)生初始熱解反應(yīng). 并且勢(shì)能達(dá)到最大值的時(shí)間隨溫度的降低而增加， 表明熱解反應(yīng)需要更多時(shí)間. 總分子數(shù)的變化也表明， 隨著模擬溫度的升高， 反應(yīng)進(jìn)行得更快， 正丙苯分子的轉(zhuǎn)化率也迅速增加. 從模擬結(jié)果來(lái)看， 主要的初始反應(yīng)是正丙苯的熱解和氫提取反應(yīng). 最初的反應(yīng)產(chǎn)物和副產(chǎn)物隨后與氧氣反應(yīng)， 釋放出大量的能量， 最終導(dǎo)致勢(shì)能急劇下降. 模擬結(jié)果顯示， 最終的主要產(chǎn)物為CO， CO2和H2O. 并且隨溫度的升高， 它們出現(xiàn)得更早， 這表明溫度升高導(dǎo)致產(chǎn)物形成時(shí)間提前. 產(chǎn)物分布表明， H2O和CO的含量隨溫度的升高而增加， 在3500 K時(shí)達(dá)到最大值［圖4（D）～（F）］. 可見(jiàn)， 在當(dāng)量比為0.5和1.0時(shí)， CO2與H2O和CO的分布具有相同趨勢(shì)， 但在當(dāng)量比2.0時(shí)， CO2的分布情況在高溫下卻表現(xiàn)出先增加后降低的趨勢(shì)， 這是因?yàn)殡S著氧化的進(jìn)行， CO2與系統(tǒng)中的自由基反應(yīng)轉(zhuǎn)化為CO和H2O， 這可能是因?yàn)榕c當(dāng)量比為0.5和1.0時(shí)的模擬相比， 氧分子的濃度較低. 當(dāng)溫度在2300～2900 K范圍內(nèi)時(shí)， 主要產(chǎn)物的最終數(shù)量急劇增加， 但隨著溫度的持續(xù)升高， 主要產(chǎn)物的數(shù)量沒(méi)有顯著變化. 這表明當(dāng)溫度達(dá)到較高范圍時(shí)， 溫度升高對(duì)氧化過(guò)程的影響有限， 這一點(diǎn)從勢(shì)能、 總分子數(shù)和正丙苯分子數(shù)的變化情況上同樣得到驗(yàn)證.

Fig.4　Time evolution of the potential energies(A), total number of species(B), n?propylbenzene(C) and main products of CO2(D), CO(E), H2O(F) at different temperatures(ρ=0.35 g/cm3, ?=2.0)

2.4　密度/壓力對(duì)正丙苯氧化的影響

為了探究密度/壓力對(duì)正丙苯氧化的影響， 分析了3個(gè)不同密度的模擬體系在溫度為3500 K時(shí)的結(jié)果（=330，=2.0）. 圖5給出了在3500 K模擬條件下， 勢(shì)能、 總分子數(shù)、 正丙苯分子數(shù)及主要產(chǎn)物隨時(shí)間的變化. 由圖5（A）可見(jiàn)， 隨著模擬體系中密度的增加， 勢(shì)能下降的趨勢(shì)也增加， 并且密度為0.15和0.35 g/cm3的模擬體系能在1 ns達(dá)到熱力學(xué)平衡. 圖5（B）和（C）分別為總分子數(shù)和正丙苯分子數(shù)變化的趨勢(shì)， 密度的增加導(dǎo)致反應(yīng)的提前開(kāi)始和燃料的快速消耗. 最終產(chǎn)物的分布也可以反映不同密度的影響. 由圖5（D）和（E）可見(jiàn)， 產(chǎn)物H2O和CO隨體系密度的增加而增加， 在0.35 g/cm3時(shí)達(dá)到最大值. 由圖5（F）可見(jiàn)， CO2首先在0.35 g/cm3出現(xiàn)， 其最大分子數(shù)出現(xiàn)在0.15 g/cm3， 并且在不同密度下CO2的分子數(shù)始終低于CO， 與Yuan等［43］的實(shí)驗(yàn)結(jié)果一致. 這主要是由于體系中氧分子含量不足， 難以支撐大量CO轉(zhuǎn)化為CO2. 通過(guò)整體觀察得知， 密度/壓力對(duì)于模擬結(jié)果的影響不如溫度大.

Fig.5　Time evolution of the potential energies(A), total number of species(B), n?propylbenzene(C), and main products of H2O(D), CO(E), CO2(F) at different densities(T=3500 K, ?=2.0)

2.5　當(dāng)量比對(duì)正丙苯氧化的影響

為了分析不同當(dāng)量比對(duì)正丙苯氧化的影響， 在3種不同當(dāng)量比（=690，=0.5；=450，=1.0；=330，=2.0）下進(jìn)行了一系列NVT-MD模擬. 并且是在3種不同密度（0.05， 0.15和0.35 g/cm3）、 溫度范圍2300～3500 K（間隔300 K）下進(jìn)行. 在每次模擬中， 確定了燃料初始消耗所需的時(shí)間以及燃料分子消耗一半所需的時(shí)間. 模擬結(jié)果如圖6（A）～（C）所示， 在較低的溫度和密度下， 燃料初始消耗以及燃料分子消耗一半所需的時(shí)間較長(zhǎng). 表明正丙苯在該模擬工況下反應(yīng)活性較小. 上述已經(jīng)討論了溫度和密度對(duì)反應(yīng)的影響， 因此， 將進(jìn)一步關(guān)注當(dāng)量比對(duì)正丙苯燃燒反應(yīng)的影響. 顯然當(dāng)量比對(duì)結(jié)果有一定的影響， 但其影響作用根據(jù)模擬溫度不同會(huì)產(chǎn)生較大變化，在高溫條件下（≥2900 K）， 結(jié)果基本沒(méi)有差別， 而在低溫條件下， 模擬結(jié)果產(chǎn)生較大差異， 與Cheng等［29］模擬的結(jié)果相類似. 圖7為在不同當(dāng)量比下模擬得到的主要產(chǎn)物隨時(shí)間的變化趨勢(shì)（NVT-MD模擬，=3500 K，=0.35 g/cm3）， 由圖7（A）和（B）可見(jiàn)， 隨著當(dāng)量比的降低， 產(chǎn)物CO2的含量分布逐步增加， CO的分子數(shù)呈現(xiàn)一個(gè)先增加后降低的趨勢(shì)（≥1.0）. 與此同時(shí)， CO2增加. 由圖7（C）可見(jiàn)， 當(dāng)=0.5和1.0時(shí)， 產(chǎn)物H2O的含量基本一致， 而當(dāng)=2時(shí)， H2O的含量卻相對(duì)較低， 這主要是氧氣含量不足所致.

Fig.6　Temperature evolution of the initiation time as well as half?time period of n?propylbenzene oxidation at different densities of 0.05(A), 0.15(B), 0.35 g/cm3(C) and equivalence ratios

Fig.7　Time evolution of the main products of CO(A), CO2(B), H2O(C) of n?propylbenzene oxidation at different equivalence ratios(T=3500 K, ρ=0.35 g/cm3)

2.6　正丙苯氧化的動(dòng)力學(xué)

正丙苯完全氧化反應(yīng)（C9H12+12O2=9CO2+6H2O）， 在2300～3500 K（間隔300 K）溫度范圍內(nèi)、 3種不同密度（0.05， 0.15和0.35 g/cm3）、 當(dāng)量比為0.5～2.0（=330， 450和690）下對(duì)正丙苯/氧氣系統(tǒng)進(jìn)行NVT-MD模擬. 每個(gè)溫度下反應(yīng)的總速率常數(shù)（， cm3·mol-1·s-1）可以通過(guò)式（2）來(lái)簡(jiǎn)化， 式（2）為正丙苯氧化的二級(jí)反應(yīng)動(dòng)力學(xué)微分方程， 主要假設(shè)是將反應(yīng)物正丙苯和氧氣均作為一級(jí)反應(yīng)處理， 即

進(jìn)一步采用不同溫度下擬合得到的速率常數(shù)通過(guò)Arrhenius方程式獲得正丙苯氧化的總反應(yīng)的表觀活化能（a）和指前因子（）等Arrhenius參數(shù)， 具體數(shù)值列于表3.

Table 3　Fitted Arrhenius parameters including Ea and A

從表3可見(jiàn)， 擬合得到的Arrhenius參數(shù)與Liang等［42］的實(shí)驗(yàn)結(jié)果相比也在可以接受范圍內(nèi)， 表明基于ReaxFF反應(yīng)力場(chǎng)的分子動(dòng)力學(xué)模擬的合理性.

3 結(jié) 論

采用DFT計(jì)算和ReaxFF模擬相結(jié)合的方法對(duì)正丙苯高溫氧化機(jī)理進(jìn)行研究， 并對(duì)初始反應(yīng)進(jìn)行了分析. 正丙苯的初始反應(yīng)主要通過(guò)兩種途徑： （1） 正丙苯側(cè)鏈上C—C和C—H鍵的斷裂； （2） 正丙苯與氧氣分子或其它自由基發(fā)生氫提取反應(yīng)， 其中， 側(cè)鏈中與芐基相連的C—C鍵的鍵能最小， 因此在反應(yīng)過(guò)程中最易斷裂， 是主要的起始反應(yīng)路徑. 依據(jù)模擬結(jié)果， 給出了正丙苯高溫氧化過(guò)程的主要反應(yīng)網(wǎng)絡(luò)， 并且發(fā)現(xiàn)無(wú)論是通過(guò)哪一種初始反應(yīng)， 隨后大部分都會(huì)生成芐基， 說(shuō)明芐基是非常重要的中間體. 在主要產(chǎn)物CO， CO2和H2O的形成機(jī)理分析中， H2O與體系中H和OH自由基緊密相關(guān)， CH2O主要是由CH3O和CH3O2的裂解形成，直接影響到最后產(chǎn)物CO和CO2等的形成，而主要的含碳穩(wěn)定產(chǎn)物CO和CO2的生成和消耗主要由體系中CHO2自由基決定. 增加體系的密度/壓力和模擬溫度可以很大程度上加快正丙苯的氧化速率、 加速產(chǎn)物的形成. 但當(dāng)溫度達(dá)到較高范圍時(shí)， 溫度升高對(duì)氧化過(guò)程的影響有限. 相對(duì)而言， 密度/壓力對(duì)結(jié)果的影響要小于溫度. 與溫度和密度/壓力的影響不同， ReaxFF模擬結(jié)果表明， 在低溫下當(dāng)量比對(duì)正丙苯消耗的反應(yīng)速率的影響較大， 而高溫條件下對(duì)模擬結(jié)果影響較小. 對(duì)正丙苯高溫氧化進(jìn)行了反應(yīng)動(dòng)力學(xué)分析， 計(jì)算得到的表觀活化能和指前因子與實(shí)驗(yàn)結(jié)果相比也在可接受的范圍內(nèi).

［1］ Liu Y. X.， Wang B. Y.， Weng J. J.， Yu D.， Richter S.， Kick T.， Tian Z. Y.，， 2018，， 53—65

［2］ Dagaut P.， El Bakali A.， Ristori A.，， 2006，（7/8）， 944—956

［3］ Xu J. Q.， Guo J. J.， Liu A. K.， Wang J. L.， Tan N. X.， Li X. Y.，， 2015，（4）， 643—652（徐佳琪， 郭俊江， 劉愛(ài)科， 王健禮， 談寧馨， 李象遠(yuǎn). 物理化學(xué)學(xué)報(bào)， 2015，（4）， 643—652）

［4］ Dooley S.， Won S. H.， Heyne J.， Farouk T. I.， Ju Y.， Dryer F. L.， Brezinsky K.，， 2012，（4）， 1444—1466

［5］ Johnston R. J.， Farrell J. T.，， 2005，（1）， 217—224

［6］ Hui X.， Das A. K.， Kumar K.， Sung C. J.， Dooley S.， Dryer F. L.，， 2012，， 695—702

［7］ Hui X.， Sung C. J.，， 2013，， 191—200

［8］ Ji C.， Dames E.， Wang H.， Egolfopoulos F. N.，， 2012，（3）， 1070—1081

［9］ Mehl M.， Herbinet O.， Dirrenberger P.， Bounaceur R.， Glaude P. A.， Battin?Leclerc F.， Pitz W. J.，， 2015，（1）， 341—348

［10］ Roubaud A.， Minetti R.， Sochet L. R..，， 2000，（3）， 535—541

［11］ Darcy D.， Mehl M.， Simmie J. M.， Würmel J.， Metcalfe W. K.， Westbrook C. K.， Curran H. J.，， 2013，（1）， 411—418

［12］ Darcy D.， Nakamura H.， Tobin C. J.， Mehl M.， Metcalfe W. K.， Pitz W. J.， Curran H. J.，， 2014，（1）， 65—74

［13］ Darcy D.， Nakamura H.， Tobin C. J.， Mehl M.， Metcalfe W. K.， Pitz W. J.， Curran H. J.，， 2014，（6）， 1460—1473

［14］ Gudiyella S.， Brezinsky K.，， 2012，（3）， 940—958

［15］ Gudiyella S.， Brezinsky K.，， 2013，（1）， 1767—1774

［16］ Anderson H.， McEnally C. S.， Pfefferle L. D.，， 2000，（2）， 2577—2583

［17］ Wang Z.， Li Y.， Zhang F.， Zhang L.， Yuan W.， Wang Y.， Qi F.，， 2013，（1）， 1785—1793

［18］ Li Y.， Cai J.， Zhang L.， Yang J.， Wang Z.， Qi F.，， 2011，（1）， 617—624

［19］ Darcy D.， Tobin C. J.， Yasunaga K.， Simmie J. M.， Würmel J.， Metcalfe W. K.， Tidjani N.， Ahmed S. S.， Westbrook C. K.， Curran H. J.，e， 2012，（7）， 2219—2232

［20］ Robinson R. K.， Lindstedt R. P.，， 2013，（12）， 2642—2653

［21］ Altarawneh M.， Dlugogorski B. Z.，， 2015，（4）， 1406—1416

［22］ Zhou C. W.， Simmie J. M.， Somers K. P.， Goldsmith C. F.， Curran H. J.，， 2017，（9）， 1890—1899

［23］ Shang Y.， Ning H.， Shi J.， Luo S. N.，， 2021，（3）， 711—717

［24］ van Duin A. C. T.， Dasgupta S.， Lorant F.， Goddard W. A.，， 2001，（41）， 9396—9409

［25］ Li X.， Zheng M.， Ren C.， Guo L.，， 2021，（15）， 11707—11739

［26］ Chenoweth K.， van Duin A. C. T.， Goddard W. A.，， 2008，（5）， 1040—1053

［27］ Chenoweth K.， van Duin A. C. T.， Dasgupta S.， Goddard W. A.，， 2009，（9）， 1740—1746

［28］ Wang Q. D.， Wang J. B.， Li J. Q.， Tan N. X.， Li X. Y.，， 2011，（2）， 217—226

［29］ Cheng X. M.， Wang Q. D.， Li J. Q.， Wang J. B.， Li X. Y.，， 2012，（40）， 9811—9818

［30］ Liu X. L.， Li， X. X.， Han S.， Qiao X. J.， Zhong B. J.， Guo L.，， 2016，（6）， 1424—1433（劉曉龍， 李曉霞， 韓嵩， 喬顯杰， 鐘北京， 郭力. 物理化學(xué)學(xué)報(bào)， 2016，（6）， 1424—1433）

［31］ Liu J. X.， Min J.， Xu H. J.， Ren H. S.， Tan N. X.，， 2022，（4），20210834（劉嘉欣， 閔杰， 許華杰， 任海生， 譚寧馨. 高等學(xué)?；瘜W(xué)學(xué)報(bào)， 2022，（4），20210834）

［32］ Berendsen H. J.， Postma J.， Gunsteren W.， DiNola A. D.， Haak J. R.，.， 1984，（8）， 3684—3690

［33］ Lindgren E. B.， Monteiro J. G. S.， dos Santos A. R.， Fleming F. P.， Barbosa A. G. H.，， 2021，， 121205

［34］ So M. R.， Rensen.， Voter A. F.，， 2000，（21）， 9599—9606

［35］ Ashraf C.， van Duin A. C. T.，， 2017，（5）， 1051—1068

［36］ Martínez L.， Andrade R.， Birgin E. G.， Martínez J. M.，， 2010，（13）， 2157—2164

［37］ Plimpton S.，， 1995，（1）， 1—19

［38］ Do?ntgen M.， Przybylski?Freund M. D.， Kro?ger L. C.， Kopp W. A.， Ismail A. E.， Leonhard K.，， 2015，（6）， 2517—2524

［39］ Frisch M. J.， Trucks G. W.， Schlegel H. B.， Scuseria G. E.， Robb M. A.， Cheeseman J. R.， Scalmani G.， Barone V.， Petersson G. A.， Nakatsuji H.， Li X.， Caricato M.， Marenich A. V.， Bloino J.， Janesko B. G.， Gomperts R.， Mennucci B.， Hratchian H. P.， Ortiz J. V.， Izmaylov A. F.， Sonnenberg J. L.， Williams?Young D.， Ding F.， Lipparini F.， Egidi F.， Goings J.， Peng B.， Petrone A.， Henderson T.， Ranasinghe D.， Zakrzewski V. G.， Gao J.， Rega N.， Zheng G.， Liang W.， Hada M.， Ehara M.， Toyota K.， Fukuda R.， Hasegawa J.， Ishida M.， Nakajima T.， Honda Y.， Kitao O.， Nakai H.， Vreven T.， Throssell K.， Montgomery J. A. Jr.， Peralta J. E.， Ogliaro F.， Bearpark M. J.， Heyd J. J.， Brothers E. N.， Kudin K. N.， Staroverov V. N.， Keith T. A.， Kobayashi R.， Normand J.， Raghavachari K.， Rendell A. P.， Burant J. C.， Iyengar S. S.， Tomasi J.， Cossi M.， Millam J. M.， Klene M.， Adamo C.， Cammi R.， Ochterski J. W.， Martin R. L.， Morokuma K.， Farkas O.， Foresman J. B.， Fox D. J.，16.，. 01， Gaussian Inc.， Wallingford CT， 2016

［40］ Luo Y. R.，， CRC Press， London， 2007

［41］ Kirkpatrick A. T.，， John Wiley & Sons， Hoboken， 2020

［42］ Liang J.， Li F.， Cao S.， Li X.， Jia M. X.， Wang Q. D.，， 2022，， 124940

［43］ Yuan W.， Li Y.， Dagaut P.， Wang Y.， Wang Z.， Qi F.，， 2017，， 178—192

Molecular Dynamics Simulation Study on High Temperature Oxidation Mechanism of-Propylbenzene

ZHOUZihao， WANGSihao， HUANGDaichuan， LIUBo， NINGHongbo*

（，，，610031，）

-Propylbenzene is a typical aromatic substitute component of Jet A， Jet A-1 and RP-3 aviation kerosene. In this work， the main oxidation reaction networks and the product distributions of-propylbenzene at different temperatures， densities and equivalence ratios were investigated by ReaxFF based on reactive molecular dynamics simulation. The reaction kinetics theory was also employed to calculate the rate constants of-propylbenzene oxidation. The results show that the consumption of-propylbenzene mainly occurs in the alkyl side chain including six C—C and C—H bond fissions of unimolecular reactions and three H-abstraction reactions by O2and other small radicals. Due to the lowest bond dissociation energy， the C—C bond fission adjacent to benzyl radical is the most important consumption channel but the contributions of all H-abstraction reactions are similar. The simulated temperature and density/pressure are positively correlated with the oxidation rate of-propylbenzene， while the effect of equivalence ratio is heavily dependent on the system temperature. Additionally， the calculated apparent activation energies and pre-exponential factors are acceptable compared to the reported experimental results.

-Propylbenzene； Reaction mechanism； High temperature oxidation； ReaxFF； Molecular dynamics simulation

2023-06-10

甯紅波, 男, 博士, 副教授, 主要從事燃燒反應(yīng)動(dòng)力學(xué)研究. E-mail: hbning@swjtu.edu.cn

四川省自然科學(xué)基金(批準(zhǔn)號(hào): 2023NSFSC1105)、 中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(批準(zhǔn)號(hào): 2682023ZTPY019)和四川省科技廳重點(diǎn)研發(fā)計(jì)劃項(xiàng)目(批準(zhǔn)號(hào): 2022YFG0033)資助.

O641

A

10.7503/cjcu20230276

2023-09-04.

Supported by the Natural Science Foundation of Sichuan Province, China(No.2023NSFSC1105), the Fundamental Research Funds for the Central Universities, China(No.2682023ZTPY019) and the Project of Key R&D Program of Sichuan Province, China(No.2022YFG0033).

（Ed.： Y， K， S）