朱 程，方驍遠(yuǎn)，喬信起

（上海交通大學(xué) 動力機(jī)械及工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 200240）

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

目前人們對詳細(xì)機(jī)理的簡化方法通常是在該詳細(xì)機(jī)理中提取出不重要的物質(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]等。還有一類方法就是在時間尺度上的簡化，即基于準(zhǔn)穩(wěn)態(tài)QSS和部分平衡PE假設(shè)確認(rèn)時間尺度非常小的組分和反應(yīng)，以代數(shù)方程替換微分方程，降低計(jì)算的代價，包括固有低維流形法（Intrinsic Low Dimensional Manifolds，ILDM）[6]、非結(jié)構(gòu)化自適應(yīng)列表法[7]和計(jì)算奇異擾動法（Computational Singular Perturbation，CSP）[8]等。

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

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

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

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

1.1 DRGEP法簡化

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

式中：

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

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

圖1 DRGEP簡化機(jī)理與目標(biāo)機(jī)理滯燃期比較

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)中的化學(xué)當(dāng)量系數(shù)和反應(yīng)速率；αaiqi為第i個基元反應(yīng)對物種a生成（或消耗）的影響；R為反應(yīng)機(jī)理系統(tǒng)中涉及到物種a的反應(yīng)個數(shù)。ROP系數(shù)反映了當(dāng)前基元反應(yīng)對該物種生成的影響程度，當(dāng)ROP系數(shù)為正時，表示當(dāng)前基元反應(yīng)促進(jìn)了該物種的生成，反之則表示會促進(jìn)該物種的消耗分解。

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

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

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

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

式中：Vi( t)為粒子i在第t代中的速度；Xpbesti和Xgbest分別表示粒子自身經(jīng)歷的歷史最好位置和群體所經(jīng)歷的全局最好位置；設(shè)置學(xué)習(xí)因子c1= c2=2. 0，慣性權(quán)重ω=0. 4；r1和r2為服從均勻分布U[0,1]的隨機(jī)數(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ù)的機(jī)理，利用該機(jī)理能計(jì)算出相應(yīng)試驗(yàn)工況下的計(jì)算值。本文中目標(biāo)函數(shù)定義為包含m個元素的目標(biāo)向量，表達(dá)式為：

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

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

1.4 SA法簡化

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

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

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

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

2 簡化結(jié)果驗(yàn)證

2.1 ST滯燃期的驗(yàn)證

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

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

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

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

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

2.2 JSR組分濃度的驗(yàn)證

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

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

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

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

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

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

2.3 PREMIX層流火焰速度的驗(yàn)證

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

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

3 結(jié)論

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

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

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

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