基于浓度敏感性分析和遗传算法的甲烷燃烧机理简化与优化

董清丽，蒋 勇，邱 榕

（中国科学技术大学火灾科学国家重点实验室，合肥，230026）

0 引言

甲烷是天然气的主要成分，是最简单的碳氢燃料，也是碳氢燃料燃烧过程中重要的中间组分。目前，国际通用的描述甲烷燃烧的化学动力学机理是GRIMECH 3.0，该机理包含53个组分，325步基元反应，将其应用于流动／燃烧耦合数值模拟时，会产生极大的刚性问题和计算效率问题［1］。因此，需要在保证燃烧特性精度的要求下，对甲烷燃烧反应的详细机理进行简化。

常用的化学机理简化方法主要有敏感性分析方法（Sensitivity Analysis，SA）［2］、浓度敏感性分析（Principal Component Analysis of Matrix S，PCAS）［3］、准稳态假设（Quasi Steady State Assumption，QSSA）［4］、固有低维流形（Intrinsic Low－Dimensional Manifolds，ILDM）［5］、直接关系图法（Directed Relation Graph，DRG）［6］等。其中浓度敏感性分析PCAS通过对反应系统矩阵的特性值和特征向量进行分析，从不同时刻反应系统的组分线性敏感系数的计算中获得重要的动力学信息，找出化学反应系统中的主要成分和限制速率反应，分析各个基元反应之间的关联性以及这些相互关联的反应对反应系统的重要程度，从而对反应机理进行简化。这种方法能够直观地分析组分、基元反应、反应条件等因素对系统反应参数变化的敏感程度，得到广泛的认可和应用。

简化机理提高了计算效率，却带来了计算结果精度的损失。为减少精度损失，一般以详细机理的计算结果为标准值，以与标准值偏差最小为目标，对简化机理的反应参数进行修正。美国教授Holland在1975年提出遗传算法［7］，这种方法模拟达尔文的优胜劣汰的物种进化机制，是一种通过优势个体的繁衍进化来搜索最优解的方法。

遗传算法从一个群体多点开始搜索，是一种全局优化搜索算法，能以较大的概率搜索到最优解。遗传算法的搜索过程依赖于对目标函数的重复求值计算，不存函数求导和函数连续性的限定，具有很大的灵活性，可广泛应用于适合具有复杂目标函数的特殊问题［8－10］。

近几年，Elliott，L［11］使用遗传算法对煤油燃烧的详细机理进行了优化，Federico［12］使用遗传算法对乙醇燃烧的详细机理进行优化，但是还没有对甲烷／空气的详细机理进行优化的研究。本文采用PCAS方法和元素流通法，对甲烷／空气燃烧机理进行简化，然后利用遗传算法，对简化机理基元反应的反应参数进行优化，并将简化机理和优化后的简化机理的计算结果与详细机理的计算结果进行比较。

1 计算模型和简化方法

1.1 计算模型

层流预混火焰模型（PREMIX）通过守恒方程和状态方程，把已燃气和未燃气的状态联系起来，模拟一维层流火焰在燃料／氧化剂预混燃料中传播的燃烧现象，可用来研究火焰结构和火焰速度等火焰特征，以及这些火焰特征和燃料参数之间的关系。

良搅拌反应器（Perfectly Stirred Reactor，PSR）代表一个重要的极限情况模型，对于高强度湍流，假设湍流混合的足够快，以至于燃烧区内的温度场和浓度场是均匀的，结合质量控制方程和气体状态方程，可用于计算反应温度、组分摩尔浓度和组分净生成速率。

1.2 元素流通法EF

元素流通法（Element Flux Analysis，EF）［13］是一种便于理解燃料氧化路径的有效工具，元素流通法测定某个特定元素从一个组分到另外一个组分的瞬时流通通量，其通量可由下式计算获得：

其中，Aijk为反应i中元素A从 系数为j的物质流入系数为k的物质的瞬时通量，nA.j和nA.j分表表示系数为j和系数为k的物质中原子A的数目，nA.i表示反应i中原子A的总数目，qi表示反应进度变量。以甲烷燃烧反应中碳元素流通为例，从组分甲烷出发，确定其主要产物，并对其产物进行元素流通分析，确定二次产物，依次类推直到最终产物CO2，以此类推，便可确定甲烷燃烧中占主导地位的氧化路径。

1.3 浓度敏感性分析PCAS

浓度敏感性分析是研究反应参数摄动对多个目标组分浓度的影响，其影响可通过下列目标函数式来衡量：

其中，k为机理中基元反应的反应速率常数，ci为机理中第i种组分的摩尔分数（z）表示参数变化后的变量值，z为特定参数，z1，z2为反应的变化区间，n为反应中组分总数。

结合局部敏感性，式（3）可简化为：

其中，kj为机理中第j个基元反应的反应速率常数，正交化的敏感性系数表示了由反应速率常数kj的微小变化导致的组分浓度ci的变化。

其中，ΔΨ＝UTΔa成为主成分。基元反应对目标函数的影响效果以特征值的方式表现出来，特征值绝对值较大的基元反应中，特征向量绝对值也较大的基元反应，就是重要反应，本文特征值和特征向量临界值分别取0.2和0.001。

2 遗传算法

图1为遗传算法流程图，从图1中可以看出，遗传算法的基本操作步骤为初始群体的生成、适应度的计算、选择、交叉和变异。

遗传算法的基本思想是从一个可能存在解集的种群开始，按照适者生存和优胜劣汰原理，进化产生最佳个体。在每一代，以这一代种群中个体的适应度最小为目标，选择出最优个体，结合自然选择学的遗传算子进行组合交叉和变异，产生新一代的种群。像自然进化一样，这个过程中产生的后生代种群，比前代种群更加适应于环境，末代种群中的最优个体，可以作为问题的近似最优解。

2.1 初始化群体

遗传算法是从初始种群开始的，种群由若干可能解个体组成。本文采用随机方法来生成许多个个体，构造群体。本文以简化机理的Ai，bi和Ei的原始值作为初始值，在初始值的90%的变化范围内搜寻最优的反应系数。每代群体中含有100个个体。

2.2 适应度

遗传算法在进化搜索以适应度函数为依据，对问题的解空间不断进行校正，以使其进化到包含最优解的领域。适应度函数作为遗传算法“适者生存”准则的环境因素，决定着优化的目标和方向，因此，适应度函数的定义十分重要。本文中适应度函数定义如下：

其中，其中f（Ai，Bi，Ei）为适应度函数，NS为反应物组分种类的总个数，K为不同燃空比下计算组分浓度的总次数，N为不同燃空比下计算层流火焰速度的总次数和分别为不同燃空比下详细机理和优化机理计算得到的组分浓度和分别为详细机理和优化机理计算得到的层流火焰速度，Wjk和为组分浓度权重和火焰速度权重，并且满足

从适应度函数的定义可以看出，通过优化机理得到的组分浓度和层流火焰速度结果越接近详细机理计算得到结果，f（Ai，Bi，Ei）的值越小，适应度越高。

2.3 选择、交叉和变异

根据计算出来的群体中各个体的适应度，选择适应度最小的个体保留下来，作为下一代个体的父代。

交叉和变异是产生新个体的重要方式，在遗传算法中必不可少。交叉是通过父代和上一代父代之间的线性组合来产生新的个体，进而产生新一代群体。本文交叉概率为0.8。在遗传算法后期，群体中的个体主要集中在某个较优的个体附近，而变异为新个体的产生提供了机会，使计算从局部最优解中跳出来。本文变异概率为0.01。

3 机理简化与优化

3.1 机理简化

甲烷grimech 3.0机理采用闭式反应器模型和一维层流预混火焰模型进行计算。PSR计算条件为压强1atm，温度1500K，燃料组成为CH4／空气。PREMIX计算条件为压强1atm，温度298K，进口质量流率0.04g／cm2.s，初始网格点6个，网格点间最大梯度GRAD＝0.4，最大斜率CURV＝0.6，燃料为CH4／空气。两种模型中，燃空比均采用0.7、0.8、0.9和1.0四种情况。

根据元素流通法，确定甲烷燃烧反应的主要链式氧化路径为：

上述机理中共有12个组分，在计算中发现，C2H3、C2H2及HCCO对火焰速率和点火延迟时间的影响较小，其相关反应可删去。综合四种燃空比下，PSR模型和PREMIX模型PCAS分析结果及层流火焰速度敏感性分析结果，最终获得17个组分、24步反应的简化机理，见表1。

3.2 机理优化

遗传算法调用PSR模型程序和PREMIX模型程序进行计算，PSR计算条件为1atm，入口温度为298K，反应器初始温度为1500K，燃空比分别为0.55、0.6、0.65、0.7、0.8、0.9、1.0，目标组分为简化机理中的17个组分。PREMIX计算条件同机理简化过程的计算条件，目标对象为层流火焰速度。

图2为优化历程中适应度的变化关系图，从图2中可以看出，优化初始阶段，适应度迅速下降，表示遗传算法能很快的找到近似最优解。随着优化代数的增加，优化过程中适应度最小值逐渐下降，并在一定的代数后基本保持不变，可以认为找到了最优解。优化后的燃烧机理见表2。

表1 24步简化机理Table 1 The reduced mechanism containing 24reactions

表2 优化后的简化机理Table 2 The optimized mechanism containing 24reactions

图2 优化历程中的适应度变化图Fig.2 The fitness in optimization process

图3 基元反应参数的相对变化率Fig.3 Relative difference of the Arrhenius parameters in basic reactions

图3为优化后的基元反应参数相对于简化机理中基元反应参数的变化率，以活化能E为例，相对变化率为（EGA－E）／E，其中EGA和E分别为优化机理和简化机理中的活化能。反应参数的总体变化范围为［－0.5～0.81］，在初始设置［－0.9～0.9］之内，说明计算过程中，设置的最优个体的搜索范围是合理的。A、b、E的平均相对变化率分别为0.05、－0.06、0.1，三个参数的相对变化率都主要集中在［－0.45～0.15］，说明优化后的反应参数还是在简化机理基元反应参数的基础上波动。

4 机理验证

图4 火焰速度随燃空比的变化Fig.4 Laminar flame speed vs.stoichiometric ratio

4.1 层流火焰速度与火焰构

图4为1atm，298K条件下基于GRIMECH 3.0、24步简化机理（24－step）和优化后的简化机理（24－GA）模拟计算的层流火焰速度随燃空比的变化，与未经优化的简化机理相比，优化后的简化机理对火焰速度的预测有了极大的改进，尤其是燃空比较小的贫燃条件下。在燃空比为1时，简化机理和优化后的机理计算得到的层流火焰速度与详细机理基本一致。

图5为1atm、298K条件下的温度和产物浓度随火焰高度的变化曲线图。由图5中可以看出，贫燃条件下，优化后的简化机理在预测产物和温度时，其计算结果更为接近详细机理的计算结果，尤其是火焰高度大于0.2cm，火焰处于稳定状态时。计算精度总体来说明显提高，但是由于适应度函数选取的局限性，无法保证每个点上的精度，某些地方还是有些偏差。

图6为1atm、298K条件下的反应物甲烷和氧气随火焰高度的变化曲线图。由图6中可以看出，贫燃条件下，简化机理和优化后的简化机理在反应物浓度时，都比较接近详细机理的计算结果。相对于详细机理的甲烷浓度和氧气浓度的预测结果，简化机理计算结果的平均计算误差分别为30%和84%，优化后简化机理计算结果的平均误差分别为15%和43%，计算精度均提高一倍。

图5 温度及产物浓度随火焰高度的变化Fig.5 Temperature and mole fraction of main products vs.flame height

图6 反应物浓度随火焰高度的变化Fig.6 Mole fraction of main reactants vs.flame height

4.2 PSR组分

图7为10atm，2000K下，滞留时间为0.0005s时，产物浓度和反应物随燃空比的变化曲线图。从图7来看，简化机理和优化后的简化机理在计算产物及反应物浓度时，都比较接近详细机理的结果，但是优化后的机理精度更高，尤其是在燃空比0.6～0.9之间，基本与详细机理的结果一样，其平均误差均小于0.2%。

图7 产物浓度和反应物浓度随燃空比的变化Fig.7 The mole fraction of main product and reactant vs.stoichiometric ratio

由于简化机理的计算结果中甲烷浓度与详细机理的计算结果相差很大，平均误差为45%，经过一定代数的优化之后，其平均计算误差达到30%，但计算精度仍然有待提高。由于适应度会随着优化代数的增加而降低，所以，这种情况可以通过优化历程的增长来改善。

5 结论

本文运用元素流通法和浓度敏感性分析法，对甲烷的燃烧反应机理GRIMECH 3.0进行简化，获得17个组分、24步基元反应的简化机理。为提高简化机理的计算精度，以甲烷／空气燃烧产生的各种产物的组分浓度和层流火焰速度作为目标，利用遗传算法的全局搜索能力，快速收敛到最优解附近，获得近似最优解，完成机理优化过程。

将优化后的简化机理进行验证，表明优化机理在计算层流火焰速度和产物浓度时，比简化机理更接近详细机理。由于优化过程中没有考虑火焰结构的影响，优化后的简化机理描述火焰结构时，无法保证每个点上的精度，某些地方还是有些偏差。

［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.