基于濃度敏感性分析和遺傳算法的甲烷燃燒機(jī)理簡(jiǎn)化與優(yōu)化

董清麗，蔣 勇，邱 榕

（中國(guó)科學(xué)技術(shù)大學(xué)火災(zāi)科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，合肥，230026）

0 引言

甲烷是天然氣的主要成分，是最簡(jiǎn)單的碳?xì)淙剂?，也是碳?xì)淙剂先紵^(guò)程中重要的中間組分。目前，國(guó)際通用的描述甲烷燃燒的化學(xué)動(dòng)力學(xué)機(jī)理是GRIMECH 3.0，該機(jī)理包含53個(gè)組分，325步基元反應(yīng)，將其應(yīng)用于流動(dòng)／燃燒耦合數(shù)值模擬時(shí)，會(huì)產(chǎn)生極大的剛性問(wèn)題和計(jì)算效率問(wèn)題［1］。因此，需要在保證燃燒特性精度的要求下，對(duì)甲烷燃燒反應(yīng)的詳細(xì)機(jī)理進(jìn)行簡(jiǎn)化。

常用的化學(xué)機(jī)理簡(jiǎn)化方法主要有敏感性分析方法（Sensitivity Analysis，SA）［2］、濃度敏感性分析（Principal Component Analysis of Matrix S，PCAS）［3］、準(zhǔn)穩(wěn)態(tài)假設(shè)（Quasi Steady State Assumption，QSSA）［4］、固有低維流形（Intrinsic Low－Dimensional Manifolds，ILDM）［5］、直接關(guān)系圖法（Directed Relation Graph，DRG）［6］等。其中濃度敏感性分析PCAS通過(guò)對(duì)反應(yīng)系統(tǒng)矩陣的特性值和特征向量進(jìn)行分析，從不同時(shí)刻反應(yīng)系統(tǒng)的組分線性敏感系數(shù)的計(jì)算中獲得重要的動(dòng)力學(xué)信息，找出化學(xué)反應(yīng)系統(tǒng)中的主要成分和限制速率反應(yīng)，分析各個(gè)基元反應(yīng)之間的關(guān)聯(lián)性以及這些相互關(guān)聯(lián)的反應(yīng)對(duì)反應(yīng)系統(tǒng)的重要程度，從而對(duì)反應(yīng)機(jī)理進(jìn)行簡(jiǎn)化。這種方法能夠直觀地分析組分、基元反應(yīng)、反應(yīng)條件等因素對(duì)系統(tǒng)反應(yīng)參數(shù)變化的敏感程度，得到廣泛的認(rèn)可和應(yīng)用。

簡(jiǎn)化機(jī)理提高了計(jì)算效率，卻帶來(lái)了計(jì)算結(jié)果精度的損失。為減少精度損失，一般以詳細(xì)機(jī)理的計(jì)算結(jié)果為標(biāo)準(zhǔn)值，以與標(biāo)準(zhǔn)值偏差最小為目標(biāo)，對(duì)簡(jiǎn)化機(jī)理的反應(yīng)參數(shù)進(jìn)行修正。美國(guó)教授Holland在1975年提出遺傳算法［7］，這種方法模擬達(dá)爾文的優(yōu)勝劣汰的物種進(jìn)化機(jī)制，是一種通過(guò)優(yōu)勢(shì)個(gè)體的繁衍進(jìn)化來(lái)搜索最優(yōu)解的方法。

遺傳算法從一個(gè)群體多點(diǎn)開(kāi)始搜索，是一種全局優(yōu)化搜索算法，能以較大的概率搜索到最優(yōu)解。遺傳算法的搜索過(guò)程依賴(lài)于對(duì)目標(biāo)函數(shù)的重復(fù)求值計(jì)算，不存函數(shù)求導(dǎo)和函數(shù)連續(xù)性的限定，具有很大的靈活性，可廣泛應(yīng)用于適合具有復(fù)雜目標(biāo)函數(shù)的特殊問(wèn)題［8－10］。

近幾年，Elliott，L［11］使用遺傳算法對(duì)煤油燃燒的詳細(xì)機(jī)理進(jìn)行了優(yōu)化，F(xiàn)ederico［12］使用遺傳算法對(duì)乙醇燃燒的詳細(xì)機(jī)理進(jìn)行優(yōu)化，但是還沒(méi)有對(duì)甲烷／空氣的詳細(xì)機(jī)理進(jìn)行優(yōu)化的研究。本文采用PCAS方法和元素流通法，對(duì)甲烷／空氣燃燒機(jī)理進(jìn)行簡(jiǎn)化，然后利用遺傳算法，對(duì)簡(jiǎn)化機(jī)理基元反應(yīng)的反應(yīng)參數(shù)進(jìn)行優(yōu)化，并將簡(jiǎn)化機(jī)理和優(yōu)化后的簡(jiǎn)化機(jī)理的計(jì)算結(jié)果與詳細(xì)機(jī)理的計(jì)算結(jié)果進(jìn)行比較。

1 計(jì)算模型和簡(jiǎn)化方法

1.1 計(jì)算模型

層流預(yù)混火焰模型（PREMIX）通過(guò)守恒方程和狀態(tài)方程，把已燃?xì)夂臀慈細(xì)獾臓顟B(tài)聯(lián)系起來(lái)，模擬一維層流火焰在燃料／氧化劑預(yù)混燃料中傳播的燃燒現(xiàn)象，可用來(lái)研究火焰結(jié)構(gòu)和火焰速度等火焰特征，以及這些火焰特征和燃料參數(shù)之間的關(guān)系。

良攪拌反應(yīng)器（Perfectly Stirred Reactor，PSR）代表一個(gè)重要的極限情況模型，對(duì)于高強(qiáng)度湍流，假設(shè)湍流混合的足夠快，以至于燃燒區(qū)內(nèi)的溫度場(chǎng)和濃度場(chǎng)是均勻的，結(jié)合質(zhì)量控制方程和氣體狀態(tài)方程，可用于計(jì)算反應(yīng)溫度、組分摩爾濃度和組分凈生成速率。

1.2 元素流通法EF

元素流通法（Element Flux Analysis，EF）［13］是一種便于理解燃料氧化路徑的有效工具，元素流通法測(cè)定某個(gè)特定元素從一個(gè)組分到另外一個(gè)組分的瞬時(shí)流通通量，其通量可由下式計(jì)算獲得：

其中，Aijk為反應(yīng)i中元素A從 系數(shù)為j的物質(zhì)流入系數(shù)為k的物質(zhì)的瞬時(shí)通量，nA.j和nA.j分表表示系數(shù)為j和系數(shù)為k的物質(zhì)中原子A的數(shù)目，nA.i表示反應(yīng)i中原子A的總數(shù)目，qi表示反應(yīng)進(jìn)度變量。以甲烷燃燒反應(yīng)中碳元素流通為例，從組分甲烷出發(fā)，確定其主要產(chǎn)物，并對(duì)其產(chǎn)物進(jìn)行元素流通分析，確定二次產(chǎn)物，依次類(lèi)推直到最終產(chǎn)物CO2，以此類(lèi)推，便可確定甲烷燃燒中占主導(dǎo)地位的氧化路徑。

1.3 濃度敏感性分析PCAS

濃度敏感性分析是研究反應(yīng)參數(shù)攝動(dòng)對(duì)多個(gè)目標(biāo)組分濃度的影響，其影響可通過(guò)下列目標(biāo)函數(shù)式來(lái)衡量：

其中，k為機(jī)理中基元反應(yīng)的反應(yīng)速率常數(shù)，ci為機(jī)理中第i種組分的摩爾分?jǐn)?shù)（z）表示參數(shù)變化后的變量值，z為特定參數(shù)，z1，z2為反應(yīng)的變化區(qū)間，n為反應(yīng)中組分總數(shù)。

結(jié)合局部敏感性，式（3）可簡(jiǎn)化為：

其中，kj為機(jī)理中第j個(gè)基元反應(yīng)的反應(yīng)速率常數(shù)，正交化的敏感性系數(shù)表示了由反應(yīng)速率常數(shù)kj的微小變化導(dǎo)致的組分濃度ci的變化。

其中，ΔΨ＝UTΔa成為主成分?；磻?yīng)對(duì)目標(biāo)函數(shù)的影響效果以特征值的方式表現(xiàn)出來(lái)，特征值絕對(duì)值較大的基元反應(yīng)中，特征向量絕對(duì)值也較大的基元反應(yīng)，就是重要反應(yīng)，本文特征值和特征向量臨界值分別取0.2和0.001。

2 遺傳算法

圖1為遺傳算法流程圖，從圖1中可以看出，遺傳算法的基本操作步驟為初始群體的生成、適應(yīng)度的計(jì)算、選擇、交叉和變異。

遺傳算法的基本思想是從一個(gè)可能存在解集的種群開(kāi)始，按照適者生存和優(yōu)勝劣汰原理，進(jìn)化產(chǎn)生最佳個(gè)體。在每一代，以這一代種群中個(gè)體的適應(yīng)度最小為目標(biāo)，選擇出最優(yōu)個(gè)體，結(jié)合自然選擇學(xué)的遺傳算子進(jìn)行組合交叉和變異，產(chǎn)生新一代的種群。像自然進(jìn)化一樣，這個(gè)過(guò)程中產(chǎn)生的后生代種群，比前代種群更加適應(yīng)于環(huán)境，末代種群中的最優(yōu)個(gè)體，可以作為問(wèn)題的近似最優(yōu)解。

2.1 初始化群體

遺傳算法是從初始種群開(kāi)始的，種群由若干可能解個(gè)體組成。本文采用隨機(jī)方法來(lái)生成許多個(gè)個(gè)體，構(gòu)造群體。本文以簡(jiǎn)化機(jī)理的Ai，bi和Ei的原始值作為初始值，在初始值的90%的變化范圍內(nèi)搜尋最優(yōu)的反應(yīng)系數(shù)。每代群體中含有100個(gè)個(gè)體。

2.2 適應(yīng)度

遺傳算法在進(jìn)化搜索以適應(yīng)度函數(shù)為依據(jù)，對(duì)問(wèn)題的解空間不斷進(jìn)行校正，以使其進(jìn)化到包含最優(yōu)解的領(lǐng)域。適應(yīng)度函數(shù)作為遺傳算法“適者生存”準(zhǔn)則的環(huán)境因素，決定著優(yōu)化的目標(biāo)和方向，因此，適應(yīng)度函數(shù)的定義十分重要。本文中適應(yīng)度函數(shù)定義如下：

其中，其中f（Ai，Bi，Ei）為適應(yīng)度函數(shù)，NS為反應(yīng)物組分種類(lèi)的總個(gè)數(shù)，K為不同燃空比下計(jì)算組分濃度的總次數(shù)，N為不同燃空比下計(jì)算層流火焰速度的總次數(shù)和分別為不同燃空比下詳細(xì)機(jī)理和優(yōu)化機(jī)理計(jì)算得到的組分濃度和分別為詳細(xì)機(jī)理和優(yōu)化機(jī)理計(jì)算得到的層流火焰速度，Wjk和為組分濃度權(quán)重和火焰速度權(quán)重，并且滿(mǎn)足

從適應(yīng)度函數(shù)的定義可以看出，通過(guò)優(yōu)化機(jī)理得到的組分濃度和層流火焰速度結(jié)果越接近詳細(xì)機(jī)理計(jì)算得到結(jié)果，f（Ai，Bi，Ei）的值越小，適應(yīng)度越高。

2.3 選擇、交叉和變異

根據(jù)計(jì)算出來(lái)的群體中各個(gè)體的適應(yīng)度，選擇適應(yīng)度最小的個(gè)體保留下來(lái)，作為下一代個(gè)體的父代。

交叉和變異是產(chǎn)生新個(gè)體的重要方式，在遺傳算法中必不可少。交叉是通過(guò)父代和上一代父代之間的線性組合來(lái)產(chǎn)生新的個(gè)體，進(jìn)而產(chǎn)生新一代群體。本文交叉概率為0.8。在遺傳算法后期，群體中的個(gè)體主要集中在某個(gè)較優(yōu)的個(gè)體附近，而變異為新個(gè)體的產(chǎn)生提供了機(jī)會(huì)，使計(jì)算從局部最優(yōu)解中跳出來(lái)。本文變異概率為0.01。

3 機(jī)理簡(jiǎn)化與優(yōu)化

3.1 機(jī)理簡(jiǎn)化

甲烷grimech 3.0機(jī)理采用閉式反應(yīng)器模型和一維層流預(yù)混火焰模型進(jìn)行計(jì)算。PSR計(jì)算條件為壓強(qiáng)1atm，溫度1500K，燃料組成為CH4／空氣。PREMIX計(jì)算條件為壓強(qiáng)1atm，溫度298K，進(jìn)口質(zhì)量流率0.04g／cm2.s，初始網(wǎng)格點(diǎn)6個(gè)，網(wǎng)格點(diǎn)間最大梯度GRAD＝0.4，最大斜率CURV＝0.6，燃料為CH4／空氣。兩種模型中，燃空比均采用0.7、0.8、0.9和1.0四種情況。

根據(jù)元素流通法，確定甲烷燃燒反應(yīng)的主要鏈?zhǔn)窖趸窂綖椋?/p>

上述機(jī)理中共有12個(gè)組分，在計(jì)算中發(fā)現(xiàn)，C2H3、C2H2及HCCO對(duì)火焰速率和點(diǎn)火延遲時(shí)間的影響較小，其相關(guān)反應(yīng)可刪去。綜合四種燃空比下，PSR模型和PREMIX模型PCAS分析結(jié)果及層流火焰速度敏感性分析結(jié)果，最終獲得17個(gè)組分、24步反應(yīng)的簡(jiǎn)化機(jī)理，見(jiàn)表1。

3.2 機(jī)理優(yōu)化

遺傳算法調(diào)用PSR模型程序和PREMIX模型程序進(jìn)行計(jì)算，PSR計(jì)算條件為1atm，入口溫度為298K，反應(yīng)器初始溫度為1500K，燃空比分別為0.55、0.6、0.65、0.7、0.8、0.9、1.0，目標(biāo)組分為簡(jiǎn)化機(jī)理中的17個(gè)組分。PREMIX計(jì)算條件同機(jī)理簡(jiǎn)化過(guò)程的計(jì)算條件，目標(biāo)對(duì)象為層流火焰速度。

圖2為優(yōu)化歷程中適應(yīng)度的變化關(guān)系圖，從圖2中可以看出，優(yōu)化初始階段，適應(yīng)度迅速下降，表示遺傳算法能很快的找到近似最優(yōu)解。隨著優(yōu)化代數(shù)的增加，優(yōu)化過(guò)程中適應(yīng)度最小值逐漸下降，并在一定的代數(shù)后基本保持不變，可以認(rèn)為找到了最優(yōu)解。優(yōu)化后的燃燒機(jī)理見(jiàn)表2。

表1 24步簡(jiǎn)化機(jī)理Table 1 The reduced mechanism containing 24reactions

表2 優(yōu)化后的簡(jiǎn)化機(jī)理Table 2 The optimized mechanism containing 24reactions

圖2 優(yōu)化歷程中的適應(yīng)度變化圖Fig.2 The fitness in optimization process

圖3 基元反應(yīng)參數(shù)的相對(duì)變化率Fig.3 Relative difference of the Arrhenius parameters in basic reactions

圖3為優(yōu)化后的基元反應(yīng)參數(shù)相對(duì)于簡(jiǎn)化機(jī)理中基元反應(yīng)參數(shù)的變化率，以活化能E為例，相對(duì)變化率為（EGA－E）／E，其中EGA和E分別為優(yōu)化機(jī)理和簡(jiǎn)化機(jī)理中的活化能。反應(yīng)參數(shù)的總體變化范圍為［－0.5～0.81］，在初始設(shè)置［－0.9～0.9］之內(nèi)，說(shuō)明計(jì)算過(guò)程中，設(shè)置的最優(yōu)個(gè)體的搜索范圍是合理的。A、b、E的平均相對(duì)變化率分別為0.05、－0.06、0.1，三個(gè)參數(shù)的相對(duì)變化率都主要集中在［－0.45～0.15］，說(shuō)明優(yōu)化后的反應(yīng)參數(shù)還是在簡(jiǎn)化機(jī)理基元反應(yīng)參數(shù)的基礎(chǔ)上波動(dòng)。

4 機(jī)理驗(yàn)證

圖4 火焰速度隨燃空比的變化Fig.4 Laminar flame speed vs.stoichiometric ratio

4.1 層流火焰速度與火焰構(gòu)

圖4為1atm，298K條件下基于GRIMECH 3.0、24步簡(jiǎn)化機(jī)理（24－step）和優(yōu)化后的簡(jiǎn)化機(jī)理（24－GA）模擬計(jì)算的層流火焰速度隨燃空比的變化，與未經(jīng)優(yōu)化的簡(jiǎn)化機(jī)理相比，優(yōu)化后的簡(jiǎn)化機(jī)理對(duì)火焰速度的預(yù)測(cè)有了極大的改進(jìn)，尤其是燃空比較小的貧燃條件下。在燃空比為1時(shí)，簡(jiǎn)化機(jī)理和優(yōu)化后的機(jī)理計(jì)算得到的層流火焰速度與詳細(xì)機(jī)理基本一致。

圖5為1atm、298K條件下的溫度和產(chǎn)物濃度隨火焰高度的變化曲線圖。由圖5中可以看出，貧燃條件下，優(yōu)化后的簡(jiǎn)化機(jī)理在預(yù)測(cè)產(chǎn)物和溫度時(shí)，其計(jì)算結(jié)果更為接近詳細(xì)機(jī)理的計(jì)算結(jié)果，尤其是火焰高度大于0.2cm，火焰處于穩(wěn)定狀態(tài)時(shí)。計(jì)算精度總體來(lái)說(shuō)明顯提高，但是由于適應(yīng)度函數(shù)選取的局限性，無(wú)法保證每個(gè)點(diǎn)上的精度，某些地方還是有些偏差。

圖6為1atm、298K條件下的反應(yīng)物甲烷和氧氣隨火焰高度的變化曲線圖。由圖6中可以看出，貧燃條件下，簡(jiǎn)化機(jī)理和優(yōu)化后的簡(jiǎn)化機(jī)理在反應(yīng)物濃度時(shí)，都比較接近詳細(xì)機(jī)理的計(jì)算結(jié)果。相對(duì)于詳細(xì)機(jī)理的甲烷濃度和氧氣濃度的預(yù)測(cè)結(jié)果，簡(jiǎn)化機(jī)理計(jì)算結(jié)果的平均計(jì)算誤差分別為30%和84%，優(yōu)化后簡(jiǎn)化機(jī)理計(jì)算結(jié)果的平均誤差分別為15%和43%，計(jì)算精度均提高一倍。

圖5 溫度及產(chǎn)物濃度隨火焰高度的變化Fig.5 Temperature and mole fraction of main products vs.flame height

圖6 反應(yīng)物濃度隨火焰高度的變化Fig.6 Mole fraction of main reactants vs.flame height

4.2 PSR組分

圖7為10atm，2000K下，滯留時(shí)間為0.0005s時(shí)，產(chǎn)物濃度和反應(yīng)物隨燃空比的變化曲線圖。從圖7來(lái)看，簡(jiǎn)化機(jī)理和優(yōu)化后的簡(jiǎn)化機(jī)理在計(jì)算產(chǎn)物及反應(yīng)物濃度時(shí)，都比較接近詳細(xì)機(jī)理的結(jié)果，但是優(yōu)化后的機(jī)理精度更高，尤其是在燃空比0.6～0.9之間，基本與詳細(xì)機(jī)理的結(jié)果一樣，其平均誤差均小于0.2%。

圖7 產(chǎn)物濃度和反應(yīng)物濃度隨燃空比的變化Fig.7 The mole fraction of main product and reactant vs.stoichiometric ratio

由于簡(jiǎn)化機(jī)理的計(jì)算結(jié)果中甲烷濃度與詳細(xì)機(jī)理的計(jì)算結(jié)果相差很大，平均誤差為45%，經(jīng)過(guò)一定代數(shù)的優(yōu)化之后，其平均計(jì)算誤差達(dá)到30%，但計(jì)算精度仍然有待提高。由于適應(yīng)度會(huì)隨著優(yōu)化代數(shù)的增加而降低，所以，這種情況可以通過(guò)優(yōu)化歷程的增長(zhǎng)來(lái)改善。

5 結(jié)論

本文運(yùn)用元素流通法和濃度敏感性分析法，對(duì)甲烷的燃燒反應(yīng)機(jī)理GRIMECH 3.0進(jìn)行簡(jiǎn)化，獲得17個(gè)組分、24步基元反應(yīng)的簡(jiǎn)化機(jī)理。為提高簡(jiǎn)化機(jī)理的計(jì)算精度，以甲烷／空氣燃燒產(chǎn)生的各種產(chǎn)物的組分濃度和層流火焰速度作為目標(biāo)，利用遺傳算法的全局搜索能力，快速收斂到最優(yōu)解附近，獲得近似最優(yōu)解，完成機(jī)理優(yōu)化過(guò)程。

將優(yōu)化后的簡(jiǎn)化機(jī)理進(jìn)行驗(yàn)證，表明優(yōu)化機(jī)理在計(jì)算層流火焰速度和產(chǎn)物濃度時(shí)，比簡(jiǎn)化機(jī)理更接近詳細(xì)機(jī)理。由于優(yōu)化過(guò)程中沒(méi)有考慮火焰結(jié)構(gòu)的影響，優(yōu)化后的簡(jiǎn)化機(jī)理描述火焰結(jié)構(gòu)時(shí)，無(wú)法保證每個(gè)點(diǎn)上的精度，某些地方還是有些偏差。

［1］Miller JA，Kee RJ.Chemical kinetics and combustion modeling［J］.Annual Reviews in Physical Chemistry，1990，41（1）：345－387.

［2］Tunlnyi T.Applications of sensitivity analysis to combustion chemistry［J］.Reliability Engineering ＆System Safety，1997，57（1）：41－48.

［3］Brow NJ，et al.Mechanism reduction via principal component analysis［J］.International Journal of Chemical Kinetics，1997，29（6）：393－402.

［4］Montgomery CJ，et al.Reduced chemical kinetic mechanism ffor hydrocarbon fuels［R］，1999，AIAA99－2220.

［5］Maas U，Pope SB.Simplifying chemical kinetics：intrinsic low dimensional manifolds in composition space［J］.Combustion and Flame，1992，88（3）：239－264.

［6］Lu TF，Law CK.A directed relation graph method mechanism reduction［J］.Proceedings of the Combustion Institute，2005，30（1）：1333－1341.

［7］John H.Adaptation in natural and artificial systems：an introductory analysis with applications to biology，control and artificial intelligence［M］.Oxford，England：The University of Michigan Press，1975.

［8］Smigrodzki R.et al.Genetic algorithm for analysis of mutations in parkinson＇s disease［J］.Artificial Intelligence in Medicine，2005，35（3）：227－241.

［9］Sakaguchi T.et al.Genetic algorithm based reactive scheduling in manufacturing system advanced crossover method for tardiness minimization problems in Mechatronics for Safety［M］.Oxford：Security and Dependability in a New Era，Elsevier，2006：207－212.

［10］Hu XB，Paolo ED.An efficient genetic algorithm with uniform crossover for air traffic control［J］.Computers＆ Operations Research，2009，36（1）：245－259.

［11］Peters N，Kee RJ.The computation of stretched laminar methane－air diffusion flames using a reduced fourstep mechanism［J］.Combustion and Flame，1987，68（1）：17－29.

［12］Federico P，Jessica L.Development of reduced and optimized reaction mechanisms based on genetic algorithms and element flux analysis［J］.Combustion and Flame，2012，159（1）：103－119.

［13］Revel J，et al.Derivation of a global chemical kinetic mechanism for methane ignition and combustion［J］.Journal de Chimie Physique，1994，91（4）：365－382.