陳菲兒，邱 越，阮 燦，王文鈺，呂興才

(上海交通大學機械與動力工程學院，上海 200240)

計算流體力學(computational fluid dynamics，CFD)模擬由于開發(fā)成本低、設(shè)計周期短，且能提供實驗研究不能提供的信息，是現(xiàn)代發(fā)動機的有效研究方法.但受限于計算機的運算能力，過于詳細的動力學機理用于數(shù)值計算會導(dǎo)致網(wǎng)格無法精細化[1].反之，過于簡化的模型對反應(yīng)路徑和中間產(chǎn)物的預(yù)測會產(chǎn)生較大的偏差[2].因此，為了更好地理解燃料在缸內(nèi)的流動、噴霧以及燃燒，又兼顧運算成本，準確可靠地進行發(fā)動機的CFD 模擬，對燃料的詳細化學動力學機理進行合理的簡化具有重要意義.

近年來針對機理簡化的理論和方法發(fā)展得非常迅速[1]，研究者已經(jīng)針對汽油[3-5]、柴油[6-7]、生物柴油[8]等模型燃料機理進行了大量的簡化工作.其中，骨架簡化方法由于其簡單高效的特點而被廣泛應(yīng)用，常見的骨架簡化方法包括直接關(guān)系圖(directed relation graph，DRG)法[9]、基于誤差傳播的直接關(guān)系圖(directed relation graph with error propagation，DRGEP)方法[10]、路徑通量分析(path flux analysis，PFA)方法[11]等.此外，Lu 等提出了兩步DRG 法[12]與加上敏感性分析的直接關(guān)系圖(directed relation graph-aided sensitivity analysis，DRGASA)方法[13]，Niemeyer 等[14]比較了DRGEP 方法中不同的路徑搜索算法的影響，發(fā)現(xiàn)使用Dijkstra 搜索算法可以得到更優(yōu)的簡化機理.

發(fā)動機的性能和排放水平與燃料燃燒息息相關(guān)，理解碳氫燃料在低溫和高壓下的化學反應(yīng)動力學對均質(zhì)壓燃(HCCI)和低溫燃燒(LTC)等先進燃燒模式的研究具有指導(dǎo)意義.丙烯作為正庚烷和異辛烷等大分子碳氫燃料燃燒的重要中間產(chǎn)物[15-17]，其分解和氧化機理直接影響大分子碳氫燃料的反應(yīng)活性和反應(yīng)路徑，對火焰?zhèn)鞑ヒ灿幸欢ǖ目刂谱饔?此外，丙烯對多環(huán)芳香烴的形成和長大有著非常重要的影響[18-19].因此，理解全工況下丙烯的燃燒特性有助于深入理解大分子碳氫燃料的反應(yīng)特性以及碳煙的生成和氧化過程，能在一定程度上揭示發(fā)動機不同燃燒階段的火焰結(jié)構(gòu)、中間產(chǎn)物和自由基的發(fā)展歷程.

基于以上背景，本文耦合MATLAB 2017 b 與CHEMKIN-PRO 15131[20]，使用DRG、DRGEP、兩步DRG 以及DRG 聯(lián)合DRGEP 4 種方法，建立了一個詳細動力學機理簡化程序.對丙烯詳細機理進行簡化，構(gòu)建了一個精度可靠、規(guī)模合理的丙烯骨架機理.選取最優(yōu)的簡化結(jié)果進行模擬，通過多個模型驗證了該機理的可靠性與簡化程序的有效性.

1 機理簡化方法

1.1 直接關(guān)系圖法

直接關(guān)系圖法[9]借鑒了圖論的思想，能從耦合的微分方程組中直接有效地量化組分間的相互影響.它將組分看作節(jié)點，組分之間的關(guān)系通過帶權(quán)的有向邊表示，從而形成一個有向圖.組分之間的關(guān)系可以定量表示為

式中：RAB為組分B 對組分A 生成和消耗的影響；Nr為機理包含的反應(yīng)總數(shù)；νA,i為組分A 在第i 個反應(yīng)的計量系數(shù)；為第i 個反應(yīng)的凈反應(yīng)速率.給定一個閾值ε，若 RABε＞，則認為去除組分B 會對組分A的生成和消耗帶來較大影響，因此若組分A 保留，則組分B 也應(yīng)被保留.通過設(shè)定初始的重要組分和閾值，可以從詳細機理里篩選出與初始組分相關(guān)的重要組分集和相應(yīng)的基元反應(yīng)，剔除不重要的反應(yīng)，最終得到一個骨架機理.

DRG 的實現(xiàn)方法是：對于每個樣本點，由零維定容模型中得到的凈反應(yīng)速率和機理系數(shù)矩陣，依據(jù)公式(1)算出組分間相關(guān)性矩陣；再根據(jù)設(shè)定的篩選閾值ε將矩陣中小于ε的數(shù)置0，構(gòu)建出一個有向圖；然后利用深度優(yōu)先搜索算法(depth first search，DFS)從初始組分集(本文中為燃料、O2、N2、H2O、CO2)出發(fā)進行搜索，得到與初始組分耦合的重要組分集合.將不同樣本點的計算結(jié)果取并集，得到最終的反應(yīng)組分集合，篩選出集合中包含的組分參與的重要反應(yīng)，由此構(gòu)建簡化機理.

1.2 基于誤差傳播的直接關(guān)系圖法

DRG 方法高效簡潔、易于實現(xiàn)，但它僅僅考慮了組分之間的直接關(guān)系，忽略了沿路徑傳播時組分間相關(guān)性的減弱，DRGEP 方法[10]則考慮了路徑長短帶來的影響.設(shè)初始組分A 經(jīng)過一系列中間組分到達B，這些中間組分構(gòu)成路徑p.在路徑p 下，B 對A 的作用rAB,p可以表示為沿程兩兩組分相關(guān)性之積：

考慮組分A 到B 的全部路徑，組分B 對A 的總作用定義為所有路徑中的最大值：

文獻[10，21]中詳細分析對比了不同 rAB計算式帶來的影響.本文在使用DRGEP 方法計算組分間相關(guān)性 rAB時，將組分A 的生成和消耗分開考慮：

DRGEP 的具體實現(xiàn)方法與DRG 相似，主要的不同在于得到相關(guān)性矩陣后，使用了Dijkstra 搜索算法來計算初始組分和其余組分之間的遠近相互作用關(guān)系，之后剔除相關(guān)性小于閾值的組分，得到需要保留的組分集，生成最終的簡化機理.需要注意的是，由于考慮了路徑長短，不同的初始組分得到的組分間相關(guān)性是不同的，在使用DRGEP 方法時，初始組分的選取會對簡化結(jié)果產(chǎn)生較大的影響.本文對丙烯的簡化選取初始組分為燃料(C3H6)、O2、N2、H2O 和CO2.

1.3 詳細動力學機理簡化程序

本文提出的機理簡化程序主要嘗試對詳細機理進行骨架簡化，簡化流程如圖1 所示，簡化過程中需要給出詳細機理，選定用于簡化的樣本工況(溫度、壓力、當量比)以及簡化機理的精度.簡化機理的精度根據(jù)燃料在零維定容模型中著火延遲的最大相對誤差(絕對值)量化，其中，著火時刻定義為反應(yīng)過程中溫度變化最大處.

圖1 機理簡化流程Fig.1 Flowchart of mechanism reduction

首先，基于MATLAB 平臺調(diào)用CHEMKIN-PRO求解器CHEM 進行詳細機理闡釋，將輸入文件中包含的機理系數(shù)矩陣、組分、熱力學數(shù)據(jù)、反應(yīng)動力學數(shù)據(jù)等信息進行編號存儲，以方便后續(xù)調(diào)用；隨后，調(diào)用零維定容模型CK Reactor Generic Closed 求解，得到所選簡化工況下的著火延遲時間以及凈反應(yīng)速率，這些是進行骨架簡化需要的動力學數(shù)據(jù)；接著，調(diào)用不同的骨架機理簡化方法，識別并剔除不重要的組分，并篩選出重要組分參與的重要反應(yīng)，從而構(gòu)建相應(yīng)的簡化機理.

刪除一些組分之后，機理隨之改變，式(1)、(4)中的 RAB值相應(yīng)改變，而考慮多步的簡化法能進一步簡化機理[22].Lu 等[12]和Poon 等[23]的研究表明，兩步DRG 法可以達到較好的效果.因此，除了DRG、DRGEP 算法，本文還使用了兩步DRG 以及DRG 聯(lián)合DRGEP 方法對詳細機理進行簡化.

設(shè)定不同的閾值可以得到一系列不同規(guī)模的簡化機理，再進行準確性的驗證，最終得到一個全局最優(yōu)的骨架機理.在兩步法中，第一步設(shè)定的閾值較小，一般在誤差快速增大處，隨后在第一步DRG 簡化的基礎(chǔ)上再應(yīng)用DRG/DRGEP 方法進行簡化.最后，對構(gòu)建的簡化機理進行驗證，包括各個工況范圍下燃料的著火延遲、LFR 平推流反應(yīng)器中的重要組分濃度變化以及一維層流火焰速度等，進一步評估簡化機理的準確性.

2 丙烯動力學機理的簡化

2.1 簡化工況與樣本點選取

本文中丙烯的詳細化學反應(yīng)動力學模型采用AramcoMech2.0 機理，共包含493 種組分和2 716 個基元反應(yīng).該機理針對激波管、快速壓縮機、射流攪拌反應(yīng)器和流動反應(yīng)器內(nèi)的大量實驗數(shù)據(jù)進行了驗證.相對上個版本，這個版本的機理包含了更多大分子反應(yīng)，對著火延遲和層流火焰速度等的驗證有明顯改善[24].值得一提的是，該機理針對丙烯擴展了實驗的工況，可適用于更高壓和更低溫工況下的預(yù)測[25].

機理的簡化通常是針對特定的工況進行的，因此簡化工況的選取相當重要.本文旨在通過計算獲得一個適用于較寬工況(尤其是較高壓力)的簡化機理，結(jié)合發(fā)動機三維CFD 模擬的常用工況，本文的簡化樣本點如表1 所示，共選取不同溫度、壓力和當量比下的90 個工況進行簡化.

其中，溫度點的選取按每個壓力和當量比下著火延遲曲線的變化，沿程選取約10 個點，覆蓋從低溫到高溫區(qū)間，由于丙烯幾乎沒有NTC 區(qū)間，因此每隔100 K 均勻選取溫度點.壓力1 MPa、當量比1 下的溫度點選取如圖2(a)所示.而對于同樣的壓力、當量比、溫度，零維定容模型反應(yīng)的凈反應(yīng)速率又會隨時間而改變，因此對于每一個樣本工況，選取沿反應(yīng)進程的 60 個時間點的數(shù)據(jù)進行簡化計算.壓力1 MPa、當量比Φ為1、溫度1 000 K 時反應(yīng)沿程選取的樣本點如圖2(b)所示.此外，在溫度和壓力變化較大的時段，樣本點更為密集，這樣能捕捉到更多反應(yīng)進程的信息.

表1 丙烯機理簡化樣本工況Tab.1 Sample conditions of propene mechanism reduction

圖2 p＝1 MPa，Φ＝1下簡化樣本點選取示意Fig.2 Schematic of selecting reduced sampling points at p＝1 MPa，Φ＝1

2.2 不同骨架簡化方法比較

本文使用所構(gòu)建的機理簡化程序，對丙烯詳細機理進行了簡化，為了找出最優(yōu)的丙烯簡化機理，比較了4 種不同簡化方法的簡化效果，包括單步DRG 方法、單步DRGEP 方法、兩步DRG 方法以及DRG 聯(lián)合DRGEP 方法.簡便起見，初始選取p＝1 MPa、Φ＝1、T 為700～1 600 K 的工況進行簡化，設(shè)定簡化機理精度為10%.不同閾值下，使用單步法獲得的丙烯簡化機理組分數(shù)變化和最大著火延遲誤差曲線如圖3 所示.

圖3(a)為使用單步DRG 方法得到的不同閾值下骨架機理組分數(shù)以及樣本點著火延遲計算的最大誤差變化曲線，在對本文的樣本點計算發(fā)現(xiàn)，隨著閾值的增大，骨架機理組分數(shù)呈現(xiàn)線性減少，而著火延遲的誤差隨閾值的增大呈非線性的變化.當閾值小于0.11 時，誤差一直保持在較小的值，而在閾值大于0.11 后，誤差呈迅速上升趨勢.在不超過設(shè)定誤差的前提下，綜合考慮骨架機理的規(guī)模與誤差，選擇一個全局最優(yōu)的閾值來達到最優(yōu)的簡化效果.最終，基于10%的設(shè)定誤差，選取閾值0.106，簡化機理組分數(shù)203，反應(yīng)數(shù)為1 187，最大誤差為9.95%.

圖3 單步法丙烯骨架機理組分數(shù)和著火延遲最大誤差隨閾值的變化Fig.3 Number of species and maximum error of ignition delay time prediction against different thresholds for propene skeletal mechanism generated by single-stage methods

圖3(b)為使用單步DRGEP 方法得到的丙烯骨架機理組分數(shù)和最大誤差隨閾值的變化.由于考慮了距離長短對兩組分之間相關(guān)性的影響，DRGEP 方法的閾值要小于DRG 方法的閾值.相比于DRG 方法，DRGEP 方法的組分數(shù)在簡化初始階段快速減少，當閾值大于0.035 時，組分減少的進程相對變緩，而著火延遲最大誤差曲線也隨閾值的增大呈非線性變化.值得一提的是，在閾值0.07 前出現(xiàn)了一個波谷，即當組分數(shù)進一步減少時，簡化的骨架機理反而具有更高的精度.Tosatto 等[26]在對JP-8 航空煤油的簡化中也發(fā)現(xiàn)了這種非線性變化現(xiàn)象.對該非線性的一個合理解釋為：RAB是絕對值，僅體現(xiàn)組分之間生成或消耗的相關(guān)性大小，而不提供正負相關(guān)性信息.某些組分的去除可能使著火延遲增大，而在進一步減少組分時，組分的去除又會使著火延遲減小，這樣總的著火延遲誤差反而變小.因此，最大誤差并不是嚴格地隨閾值的增大而增大.最終，單步DRGEP方法下選取閾值為0.068，簡化機理包含220 組分和1 253 個反應(yīng)，最大誤差為9.30%.

圖4(a)、(b)分別為兩步法中第二步使用DRG和DRGEP 方法進行簡化的計算結(jié)果.兩步法中，首先根據(jù)圖3(a)中單步DRG 方法計算的誤差曲線，選取第一步簡化閾值為0.106，得到一個簡化機理，再在新機理基礎(chǔ)上進行第二次簡化計算.由圖可見，使用兩步法進行簡化時，誤差開始隨閾值變化不明顯，而后有一個很明顯的階躍，說明此時有較多的耦合組分，進一步刪除組分會帶來很大的誤差.此外，在本文對丙烯的簡化過程中發(fā)現(xiàn)，第二步使用DRG 方法和DRGEP 方法的簡化效果差異不明顯，這可能是因為用于二次簡化的機理規(guī)模已經(jīng)較小，可簡化的空間不大.考慮到10%的設(shè)定誤差，對兩步DRG 方法最終選取閾值0.1，簡化機理組分數(shù)為189，最大誤差為9.94%.同樣地，對DRG 聯(lián)合DRGEP 方法最終的閾值選取 0.175，簡化機理組分 187，最大誤差為9.94%.表2 列出了分別用四種方法得到的簡化機理的相關(guān)信息.

圖4 兩步法丙烯骨架機理組分數(shù)和著火延遲最大誤差隨閾值的變化Fig.4 Number of species and maximum error of ignition delay time prediction against different thresholds for propene skeletal mechanism generated by twostage methods

總的來說，基于DRG 的骨架簡化方法均可以對原機理有很大程度上的簡化，本例中當組分數(shù)減少50%時，最大誤差也不超過0.5%.而當組分數(shù)繼續(xù)減少到一定程度，最大誤差會快速上升.這樣的變化趨勢在單步DRG 方法、兩步DRG 方法以及DRG 聯(lián)合DRGEP 方法中均有所體現(xiàn).綜合考慮，本文選擇第四種方法，即DRG 聯(lián)合DRGEP 算法進行全局工況下的機理簡化.

表2 不同簡化算法得到的簡化機理大小以及最大誤差Tab.2 Size of reduced mechanisms and the corresponding maximum errors obtained by using different reduction methods

2.3 丙烯全局簡化機理構(gòu)建

根據(jù)對不同簡化方法的分析，最終使用DRG 聯(lián)合DRGEP 方法分別對不同壓力、當量比和溫度下的樣本工況進行簡化，將每組得到的簡化組分集取并集，構(gòu)建全局簡化機理.本文發(fā)現(xiàn)簡化機理對壓力和當量比具有較好的可拓展性.經(jīng)過計算發(fā)現(xiàn)，上節(jié)中對p＝1 MPa，Φ＝1 下簡化得到的丙烯骨架機理在全局工況中均有較好的適用性，除p＝1 MPa、Φ＝2、T＝700 K 一個樣本工況外，其余所有簡化工況下的最大誤差都小于10%，各工況下的平均誤差約為2%.這大大節(jié)約了計算機的計算成本.由此可知，機理對溫度的敏感性大于對壓力和當量比的敏感性，在對詳細機理進行簡化時，可以首先對某一壓力、當量比的工況進行簡化，再外推至全局簡化點的驗證，這樣可以有效提高計算效率.最后，本文采用DRG 聯(lián)合DRGEP 算法，將原機理簡化掉約62%的組分數(shù)和58%的反應(yīng)數(shù)，獲得了一個包含187 種組分、1 139 個基元反應(yīng)的丙烯骨架機理，用于最終的驗證.

3 丙烯骨架機理的驗證

3.1 著火特性

為了驗證所得的丙烯骨架機理的有效性，本文計算了零維定容模型反應(yīng)器中不同壓力當量比下的著火延遲隨溫度的變化，并與詳細機理的計算結(jié)果進行比較.圖5(a)～(c)分別展示了在當量比為0.5、1.0、2.0 下骨架機理與詳細機理的著火延遲隨溫度的變化趨勢.由圖可見，丙烯骨架機理在壓力1～4 MPa，溫度680～1 680 K，當量比0.5～2.0 的范圍內(nèi)，簡化機理的計算點幾乎和原始機理計算點完全重合，說明簡化機理在各個工況下對丙烯的著火延遲均有較好的預(yù)測能力.

圖5 不同工況下丙烯骨架機理與詳細機理點火延遲比較Fig.5 Comparison of ignition delay time between propene skeletal mechanism and detailed mechanism in different working conditions

3.2 平推流反應(yīng)器

為了進一步驗證所構(gòu)建的丙烯骨架機理的有效性，本文計算了等壓絕熱假設(shè)下的平推流反應(yīng)器(LFR)模型中，丙烯詳細機理與骨架機理對重要組分濃度隨停留時間的變化(反應(yīng)工況：p＝1 MPa，Φ＝1，初始摩爾分數(shù)0.33% C3H6、1.49% O2、98.18%N2)，并與Burke 等[25]在普林斯頓大學的變壓流動反應(yīng)器(variable pressure flow reactor，VPFR)上獲得的實驗結(jié)果進行比較，如圖6 所示.為了更好地與實驗數(shù)據(jù)進行對比，本文對詳細機理與骨架機理的計算結(jié)果進行了 0.6 s 的時移[27]，以匹配丙烯燃料的消耗.從圖中可以看出，骨架機理能夠較好地預(yù)測丙烯氧化反應(yīng)中組分的變化，且骨架機理與詳細機理的計算結(jié)果比較吻合，僅在C2H4的預(yù)測上出現(xiàn)一些偏差.但計算結(jié)果仍然和實驗結(jié)果存在一定偏離，表明詳細機理還需要進一步修正.

圖6 丙烯骨架機理(虛線)和詳細機理(實線)對LFR 模型中重要組分濃度變化的預(yù)測及實驗數(shù)據(jù)(點)對比，p＝1 MPa，Φ＝1Fig.6 Comparison of mole fractions of major species in LFR model between predictions using propene skeletal mechanism(dashed line) and detailed mechanism(solid line) and experimental data at p＝1 MPa，Φ＝1

3.3 層流火焰速度

已有的研究[28]表明，機理中大分子的反應(yīng)主要影響著火延遲，而小分子則更多地控制燃料的活化特性，決定燃料的火焰?zhèn)鞑ニ俣?也就是說，以著火延遲為簡化目標獲得的骨架機理并不能保證其對燃料的活化特性也有較好的預(yù)測.為了驗證丙烯骨架機理在模擬火焰的輸運特性時的可靠性，本文分別使用了丙烯詳細/骨架機理進行不同壓力下層流火焰速度的計算.初始溫度T＝298 K，壓力1～4 MPa 下不同稀釋度的丙烯層流火焰速度如圖7 所示，可以看到，簡化機理與詳細機理的模擬結(jié)果基本重合，說明本文發(fā)展的丙烯骨架機理對1-D 層流火焰速度有較好的再現(xiàn)性.

值得一提的是，燃料機理的規(guī)模直接決定了數(shù)值模擬的運算時間，在使用丙烯骨架機理計算層流火焰速度時，發(fā)現(xiàn)計算機的運算速度有明顯提升.Lu 和Law[29]研究了計算機運算成本對機理的組分數(shù)之間的依賴關(guān)系，隱式求解的計算成本是機理包含的組分數(shù)的三次函數(shù).由此可計算得，使用本文所發(fā)展的丙烯骨架機理進行1-D 層流火焰速度預(yù)測，所花的時間縮短至使用詳細機理所花時間約1/18，大大節(jié)約了計算成本，有助于提高后續(xù)CFD 計算的效率.

圖7 不同壓力下丙烯骨架機理與詳細機理預(yù)測層流火焰速度隨當量比變化的比較Fig.7 Comparison of changes in laminar flame speed with equivalence ratio predicted by using propene skeletal mechanism and detailed mechanism at different pressures

4 結(jié)論

(1)基于DRG 方法的骨架簡化方法可以快速、有效地減小詳細機理的規(guī)模.當組分減少50%時，最大著火延遲誤差僅為0.5%.但當組分數(shù)減少到一定程度時，誤差會出現(xiàn)階躍.使用DRGEP 方法時，組分數(shù)在簡化初始階段減少得更快.兩步法相比單步法在誤差范圍內(nèi)可以進一步縮小機理規(guī)模，但可能由于丙烯詳細機理規(guī)模較小，考慮運算成本后使用兩步法的優(yōu)勢不明顯.

(2)選用DRG 聯(lián)合DRGEP 方法，基于p＝1 MPa，Φ＝1.0 工況下簡化的丙烯骨架機理在全局范圍內(nèi)(680～1 680 K，1～4 MPa，Φ為0.5～2.0)均有較好的適用性，能有效減少運算成本.該機理在多個工況下的著火延遲均與原機理較吻合，對較高壓力下層流火焰速度隨當量比變化也有較為準確的預(yù)測效果，并且能較好地模擬丙烯在LFR 模型中重要組分濃度隨停留時間的變化，這也進一步驗證了本文所提出的機理簡化方法的有效性與可行性.