對異十六烷詳細機理的簡化研究

朱 程，方驍遠，喬信起

（上海交通大學 動力機械及工程教育部重點實驗室，上海 200240）

一直以來，碳氫燃料尤其是重烴燃料的化學反應(yīng)動力學都是汽車發(fā)動機燃燒方向領(lǐng)域的一個研究重點。異十六烷是多種化石或生物衍生燃料中的異構(gòu)烷烴代表性組分，是一種高度支鏈化的烷烴，其十六烷值僅為15，常用作構(gòu)建模型燃料的十六烷值調(diào)節(jié)組分。由于異十六烷的詳細反應(yīng)機理規(guī)模較大，在實際的燃料化學反應(yīng)系統(tǒng)求解中，“剛性”問題相對突出，使其在燃燒過程的數(shù)值模擬中難以直接使用。因此，對其詳細化學反應(yīng)機理進行簡化的意義顯得十分重要[1]。

目前人們對詳細機理的簡化方法通常是在該詳細機理中提取出不重要的物質(zhì)和反應(yīng)，包括敏感性分析法（Sensitivity Analysis，SA）[2]、主成分分析法（Principal Component Analysis，PCA）[3]、生成速率分析法（Rate of Production Analysis，ROP）[4]、直接關(guān)系圖法（Directed Relation Graph，DRG）[3]、反應(yīng)路徑通量分析法[5]等。還有一類方法就是在時間尺度上的簡化，即基于準穩(wěn)態(tài)QSS和部分平衡PE假設(shè)確認時間尺度非常小的組分和反應(yīng)，以代數(shù)方程替換微分方程，降低計算的代價，包括固有低維流形法（Intrinsic Low Dimensional Manifolds，ILDM）[6]、非結(jié)構(gòu)化自適應(yīng)列表法[7]和計算奇異擾動法（Computational Singular Perturbation，CSP）[8]等。

很多研究表明，單一的簡化方法往往不能夠?qū)崿F(xiàn)對化學反應(yīng)機理的有效簡化。PISTCH等[9]提出的基于誤差傳播的直接關(guān)系圖法（DRGEP）能夠快速高效地識別冗余反應(yīng)機組份并刪除，所以本文通過采用DRGEP法對異十六烷的詳細機理進行初步簡化。 ROP法可以僅從反應(yīng)物或生成物反應(yīng)速率的數(shù)量級上識別出機理中的次要反應(yīng)，易于理解及代碼化，所以在對異十六烷機理簡化過程中采用ROP法刪除機理中的次要組分。SA法能夠?qū)Ω鱾€基元反應(yīng)計算目標的重要性提供直觀信息，其對機理的簡化性能可靠，所以將SA法作為最后的簡化方法，對主要的計算目標進行敏感性分析來剔除次要反應(yīng)，保證簡化達到預(yù)期規(guī)模。

簡化機理一方面雖然精簡了機理的規(guī)模并提高了計算效率，另一方面卻在計算結(jié)果精度上造成一定損失。為了克服這一問題，人們一般以計算詳細機理的結(jié)果或者試驗中得到的相關(guān)試驗值作為標準值，以與該標準值的最小偏差為目標，對得到的簡化機理反應(yīng)參數(shù)進行優(yōu)化。近年來，隨著計算機與智能技術(shù)的發(fā)展，智能優(yōu)化算法在各領(lǐng)域受到很大的關(guān)注。1975年美國教授HOLLAND提出遺傳算法；美國電氣工程師EBERHART博士[10]和心理學家KENNEDY博士共同提出粒子群優(yōu)化（PSO）算法；近幾年，PERINI等[5]和董清麗等[11]利用遺傳算法分別對煤油、乙醇和甲烷/空氣燃燒的簡化機理進行了優(yōu)化。而粒子群優(yōu)化作為一種基于啟發(fā)式進化智能計算的優(yōu)化算法，其概念簡明、實現(xiàn)方便、收斂速度快，是一種高效的隨機搜索方法。近年來，雖然該方法已經(jīng)成為新的研究熱點，但在化學反應(yīng)動力學的應(yīng)用相對比較少。本文通過合理設(shè)置PSO優(yōu)化算法的自變量及目標函數(shù)對簡化機理中的反應(yīng)參數(shù)進行優(yōu)化，使簡化機理相應(yīng)的計算精度大大提高。

1 異十六烷氧化簡化機理構(gòu)建

本文以RANZI等[12]構(gòu)建的包含183組分和5 744步反應(yīng)的化石及生物衍生燃料綜合機理為基礎(chǔ)，采用DRGEP識別并刪除與異十六烷反應(yīng)不相關(guān)的冗余反應(yīng)和組分；采用ROP法刪除次要組分及其對應(yīng)的反應(yīng)；針對ROP法對滯燃期計算帶來的誤差，采用PSO算法對敏感反應(yīng)參數(shù)進行優(yōu)化；最后對多種反應(yīng)器的不同目標參數(shù)進行敏感性分析，實現(xiàn)對異十六烷機理的最終簡化。簡化過程中選取的目標參數(shù)包括激波管（Shock Tube，ST）中的滯燃期，射流攪拌反應(yīng)器（Jet Stirred Reactor，JSR）中重要組分的濃度以及層流預(yù)混火焰（PREMIX）的火焰?zhèn)鞑ニ俣取?/p>

1.1 DRGEP法簡化

在傳統(tǒng)的DRG中，采用組分A對組分B生成率的正規(guī)化貢獻rAB來表示組分A和B之間的耦合關(guān)系。2008年，PISTCH等[9]提出基于誤差傳播的直接關(guān)系圖法重新定義了組分之間的耦合關(guān)系式：

式中：

式（1）和（2）中，?υA,iω表示第i個反應(yīng)下關(guān)于組分A的凈化學反應(yīng)計量系數(shù)，而ωi表示第i個反應(yīng)的凈化學反應(yīng)速率。簡化時，以詳細機理計算的滯燃期為目標參數(shù)，采用二分法找出最佳的簡化閾值，最終其被定義為7%，將iC16H34、CO2、CO、H2O、N2、O2等反應(yīng)物和最終生成物設(shè)置為保留組分，得到了包含125組分2 421步的異構(gòu)十六烷簡化機理。圖1為選取φ=1.0，p=2 026.5 kPa下覆蓋750 K到2 000 K溫度范圍的5個工況點詳細機理與簡化機理計算結(jié)果的對比。

可以看出，初步簡化后的機理引起滯燃期相對于目標機理的誤差在自定義閾值內(nèi)，5個工況下的目標參數(shù)與目標機理基本吻合，其中在750 K時的相對誤差略大，表明了簡化機理在低溫工況下滯燃期的敏感性相對較高，更容易引起計算偏差。

圖1 DRGEP簡化機理與目標機理滯燃期比較

1.2 ROP法簡化

ROP法可以用來分析得到不同基元反應(yīng)對反應(yīng)物種生成或消耗的重要性指數(shù)，能快速找出反應(yīng)系統(tǒng)中主要的基元反應(yīng)。組分a的生成（消耗）速率pa（數(shù)值上等于凈反應(yīng)速率ωa），以及不同基元反應(yīng)對組分a生成（消耗）速率的重要性指數(shù)E可以表示為：

式中：αai和qi分別為第i個基元反應(yīng)中的化學當量系數(shù)和反應(yīng)速率；αaiqi為第i個基元反應(yīng)對物種a生成（或消耗）的影響；R為反應(yīng)機理系統(tǒng)中涉及到物種a的反應(yīng)個數(shù)。ROP系數(shù)反映了當前基元反應(yīng)對該物種生成的影響程度，當ROP系數(shù)為正時，表示當前基元反應(yīng)促進了該物種的生成，反之則表示會促進該物種的消耗分解。

本文將組分a（a為任意分析組分）參與的所有消耗反應(yīng)中最小的重要性指數(shù)設(shè)為參考值。設(shè)定閾值為ε進行簡化，即當組分a參與的其它消耗反應(yīng)計算的重要性指數(shù)時（其中都為負數(shù)），對應(yīng)的消耗反應(yīng)被視為“次要反應(yīng)”，需要被刪除，否則被保留在簡化機理中。經(jīng)過多次嘗試后發(fā)現(xiàn)ROP閾值被設(shè)定為1%時，簡化機理關(guān)于滯燃期的計算誤差較小，簡化結(jié)果整體較理想，僅低溫區(qū)的計算誤差較大。ROP閾值1%簡化后的機理規(guī)模為86種組分和759步反應(yīng)。

1.3 PSO算法優(yōu)化反應(yīng)參數(shù)

針對ROP法簡化后機理在低溫區(qū)存在較大誤差的問題，采用粒子群優(yōu)化算法對敏感反應(yīng)的速率常數(shù)進行優(yōu)化。本文通過對滯燃期敏感性分析選取的敏感反應(yīng)為燃料高溫裂解反應(yīng)R1、烷基加氧生成過氧烷基的反應(yīng)R15和過氧化氫酮裂解反應(yīng)R19。

PSO是通過初始化產(chǎn)生多個隨機粒子，經(jīng)過不斷迭代找到最佳粒子的過程。圖2是本文采用粒子群算法的流程圖，其中速度和位置信息更新公式為：

式中：Vi( t)為粒子i在第t代中的速度；Xpbesti和Xgbest分別表示粒子自身經(jīng)歷的歷史最好位置和群體所經(jīng)歷的全局最好位置；設(shè)置學習因子c1= c2=2. 0，慣性權(quán)重ω=0. 4；r1和r2為服從均勻分布U[0,1]的隨機數(shù)；X和V維數(shù)為N×n，N為粒子群規(guī)模大小，n為自變量粒子維度。

圖2 粒子群優(yōu)化算法流程圖

本文對粒子種群的定義為X =(X1, X2,… ,XN)，對每個粒子向量Xi=(xi1, xi2,… ,xin)，粒子變量xij=K KK，其中反應(yīng)速率常數(shù)K = A Tnexp ( ? E / R T )，

ijIJKIJ為第i組初始變量對應(yīng)的第j步反應(yīng)的初始反應(yīng)速率常數(shù)，Kij為優(yōu)化前第i組變量對應(yīng)的第j步反應(yīng)的反應(yīng)速率常數(shù)。設(shè)粒子變量的取值范圍為[0.1,10]。

對上述粒子向量中的每個xij都能通過逆推導(dǎo)得出阿累尼烏斯參數(shù)Ai，Bi和Ei，從而每個粒子向量Xi=(xi1, xi2,… ,xin)都能組成n步的修改過參數(shù)的機理，利用該機理能計算出相應(yīng)試驗工況下的計算值。本文中目標函數(shù)定義為包含m個元素的目標向量，表達式為：

本文中隨機生成了N=20組粒子，進化代數(shù)T=120，目標函數(shù)為m=9個正交工況點下的滯燃期計算結(jié)果與ROP法簡化機理計算結(jié)果的偏差。優(yōu)化結(jié)果見表1。

表1 粒子群對敏感反應(yīng)參數(shù)的優(yōu)化結(jié)果

1.4 SA法簡化

上述ROP法簡化后的機理規(guī)模仍然相對較大，為此采用SA分析法進一步簡化。首先將機理中特定組分對應(yīng)的基元反應(yīng)逐個移除，并計算移除帶來的誤差變化：

式中：δB,ind為組分B參與的某個反應(yīng)被移除后機理計算的目標參數(shù)相對原始機理的誤差；δROP為使用ROP機理計算時相對原始機理的誤差。隨后針對組分B按照不同的δB值從大到小依次刪除對應(yīng)的基元反應(yīng)，每移除一個反應(yīng)后評估總誤差，直到最大誤差達到定義的誤差限值。

將燃料iC16H34和燃料在反應(yīng)中主路徑裂解、異構(gòu)、氧化和脫氫反應(yīng)后的重要產(chǎn)物以及氧化劑O2、燃燒產(chǎn)物CO2、CO、H2O等參與的反應(yīng)保留不作分析以節(jié)省計算時間，對其余組分逐個進行敏感性分析。目標工況為覆蓋高中低3個水平下當量比、溫度和壓力的9個滯燃期正交工況點，閾值取該9個點的平均相對誤差為11%進行刪除簡化，得到了187步反應(yīng)的異十六烷機理。

以滯燃期為目標參數(shù)得到的187步簡化機理僅對滯燃期具有較高的計算精度，并不能保證機理關(guān)于燃料其它燃燒特征的計算可靠性。因此，隨后針對上一步被簡化操作刪除的所有反應(yīng)逐個分析其對JSR中重要組分濃度及層流預(yù)混火焰?zhèn)鞑ニ俣鹊拿舾行?。JSR的目標工況同樣為9個正交的工況點，閾值定義為平均相對誤差2.1%。PREMIX由于計算網(wǎng)格較多，計算復(fù)雜比較難出結(jié)果，所以僅針對ωφ=1.1進行分析，并取閾值為8%。最終在原簡化機理的基礎(chǔ)上增補了112步反應(yīng)，得到了83組分299步反應(yīng)的異十六烷簡化機理。

2 簡化結(jié)果驗證

2.1 ST滯燃期的驗證

圖3～5是基于異十六烷原綜合機理和本文構(gòu)建的簡化機理分別在Chemkin-Pro的激波管模型中模擬得到的滯燃期隨初始溫度的變化，其中選取三個計算當量比ωφ=0.5、1.0、1.5，兩個計算壓力ωp =1 013.25 kPa、4 053 kPa，計算溫度選擇在700～2 000 K范圍內(nèi)。

圖3 簡化機理和詳細機理在激波管中滯燃期隨溫度變化預(yù)測值的比較，φ=0.5

圖4 簡化機理和詳細機理在激波管中滯燃期隨溫度變化預(yù)測值的比較，φ=1.0

圖5 簡化機理和詳細機理在激波管中滯燃期隨溫度變化預(yù)測值的比較，φ=1.5

由圖3～5可知，在700～800 K和1 300～2 000 K范圍內(nèi)兩種機理的預(yù)測值稍有偏差，其中在低溫區(qū)比在高溫區(qū)的誤差更大。簡化機理在800～1 300 K范圍內(nèi)與詳細機理的預(yù)測值基本吻合，表明簡化機理在該溫度區(qū)間內(nèi)的計算精度更高。在三個當量比下4 053 kPa曲線和1 013.25 kPa曲線相比，簡化機理與詳細機理的吻合程度都要高，表明該簡化機理在高壓4 053 kPa情況下的計算精度相對于低壓1 013.25 kPa情況下的計算精度更高。但從整體看，簡化機理與原綜合機理關(guān)于滯燃期的預(yù)測結(jié)果基本一致。

2.2 JSR組分濃度的驗證

采用Chemkin-Pro中常壓、等溫的良攪拌器模型對射流攪拌器（JSR）進行模擬，模擬工況根據(jù)DAGAUT等[13]的試驗選取φ=0.5、1.0、2.0三個當量比，計算溫度范圍為770～1 030 K，計算壓力為1 013.25 kPa。圖6～8分別為在φ=0.5、1.0、2.0的情況下，基于299步簡化機理和詳細機理模擬計算的組分濃度隨溫度的變化結(jié)果，其中重點關(guān)注燃料iC16H34、氧化產(chǎn)物CO2/CO和H2O以及其它重要中間產(chǎn)物CH4、CH2O等共6種組分。從圖中可以看出，簡化機理和詳細機理對6種組分的濃度演變過程預(yù)測誤差都很小，在三種工況下兩者的計算結(jié)果都很接近。

圖6 JSR 中主要組分濃度的詳細機理預(yù)測值（實線）、簡化機理預(yù)測值（虛線）與試驗值（實心符號）[16]的比較，p=1 013.25 kPa，φ=0.5

圖7 JSR 中主要組分濃度的詳細機理預(yù)測值（實線）、簡化機理預(yù)測值（虛線）與試驗值（實心符號）[13]的比較，p=1 013.25 kPa，φ=1.0

圖8 JSR 中主要組分濃度的詳細機理預(yù)測值（實線）、簡化機理預(yù)測值（虛線）與試驗值（實心符號）[13]的比較，p=1 013.25 kPa，φ=2.0

對比圖6～8可知，在富氧情況下詳細機理與簡化機理的預(yù)測值相比其它兩種工況下的吻合度要高，證明該簡化機理在富氧條件下對詳細機理關(guān)于重要組分濃度的計算精度損失最小。且在T=850～950 K溫度區(qū)間隨著當量比的增加，CO、H2O濃度的計算精度損失增大；CH2O的濃度在T＜900 K時相對詳細機理的濃度計算誤差最大，其次是CH4；而在高溫區(qū)，其濃度隨著貧氧程度增大其計算誤差也增大。整體上，詳細機理和簡化機理在關(guān)于JSR重要組分濃度計算上較為一致，簡化機理的計算相對詳細機理有較高的還原度。

而在與DAGAUT等的試驗數(shù)據(jù)對比上，iC16H34的吻合度最高，其次是CO、CO2、H2O等組分，且隨著當量比增加，iC16H34、H2O、CO計算濃度和試驗值的吻合度變差，而CH4、CO2、CH2O濃度計算值吻合度變高?？傮w而言，本文的簡化機理對射流攪拌器組分濃度的預(yù)測是可靠的。

2.3 PREMIX層流火焰速度的驗證

采用Chemkin-Pro的層流預(yù)混火焰模型，根據(jù)試驗條件在未燃混合氣溫度Tu=443 K、壓力p=101.325 kPa的情況下驗證層流火焰速度隨當量比的變化，圖9是計算結(jié)果和試驗數(shù)據(jù)的對比圖，實心點符號是來自清華大學李博[14]的試驗數(shù)據(jù)，虛線是簡化機理的模擬曲線?？梢钥闯?，在φ=0.9到φ=1.1之間，預(yù)測值略小于試驗結(jié)果，但在φ=0.8到φ=0.9以及φ=1.1到φ=1.4兩段區(qū)間內(nèi)，299步簡化機理計算值與試驗值高度一致，在φ=0.8到φ=1.4之間的平均誤差小于2%，可認為簡化機理的層流火焰速度預(yù)測值較為精確地反映了實際的層流火焰速度。

圖9 火焰?zhèn)鞑ニ俣鹊念A(yù)測值（虛線）與試驗值[14]的比較，p=101.325 kPa，Tu = 443K

3 結(jié)論

（1）本文基于DRGEP、ROP和SA三種方法對183組分5 744步反應(yīng)的綜合機理進行多次簡化，簡化過程中采用PSO算法對簡化機理反應(yīng)參數(shù)進行優(yōu)化使計算結(jié)果偏差達到最小，最終得到83組分299步反應(yīng)的異十六烷簡化機理。關(guān)于構(gòu)建的異十六烷簡化機理動力學文件、對機理計算時使用的熱力學參數(shù)文件、傳質(zhì)參數(shù)文件，以及應(yīng)用到的簡化方法、優(yōu)化算法程序，讀者都可以通過訪問作者GitHub主頁中的iso-cetane-reduction-mechanism倉庫下載參考，該倉庫地址為：https://github.com/dutzhucheng/iso-cetane-reduction-mechanism。

（2）DRGEP方法對剔除冗余反應(yīng)和組分的效果明顯，而對標記的主要和次要組分中，ROP方法邏輯直觀、結(jié)果理想，SA方法在移除次要反應(yīng)過程中占主導(dǎo)地位，PSO優(yōu)化算法能夠?qū)Χ喾磻?yīng)模型計算結(jié)果偏差進行修正優(yōu)化，優(yōu)化效果顯著。其中三種簡化方法按照其簡化效率從小到大的順序進行，有效地減小了單一簡化方法的不確定性和目標依賴性，所以該簡化流程同樣可以有效地對其它燃料尤其是重碳純?nèi)剂辖M分如正十二烷等的詳細機理進行簡化，具有較高的普適性。

（3）在對DRGEP、ROP及SA法的閾值進行選擇時，反應(yīng)機理中存在相當一部分的重要反應(yīng)會被過大的閾值簡化，而閾值選擇過小時的簡化效率不高。采用二分法分別定義各簡化方法的可行閾值，一方面可以減小簡化過程的計算量，另一方面可以使反應(yīng)機理在盡可能簡潔的同時保持相關(guān)的計算精度。

（4）最終得到的異十六烷簡化機理模擬的激波管滯燃期、射流攪拌反應(yīng)器組分濃度和層流預(yù)混火焰速度與原機理或試驗值吻合度較高，表明本文采用的簡化方法及簡化流程是合理可靠的。