基于特征值分析的骨架机理获取方法

文 斐 钟北京

(清华大学航天航空学院,北京100084)

基于特征值分析的骨架机理获取方法

文 斐 钟北京*

(清华大学航天航空学院,北京100084)

基于特征值分析法,建立了一套复杂化学反应动力学模型简化方法,并采用该方法对甲烷空气燃烧的详细化学反应动力学模型进行了简化.从GRI1.2得到了21个组分,83个反应的骨架机理.该机理与详细的GRI1.2机理和DRM19机理在不同化学计量比和不同压力下对比了点火延迟时间,结果表明简化机理能有效地再现详细反应动力学模型的反应机理,并具有更高的计算精度.从GRI3.0简化得到两种骨架机理分别为26个组分、120个反应和30个组分、140个反应.这两个机理都能很好地对火焰传播速度以及主要组分和NO浓度分布进行反应动力学模拟.

化学机理简化;骨架机理;特征值分析;甲烷;点火延迟;火焰传播速度

1 引言

随着计算机技术能力的发展,数值方法应用越来越广泛.但燃烧数值模拟仍然存在着很多难题,主要包括两方面,一是计算量太大,由于燃烧反应的复杂性,燃烧的详细反应模型往往包含着上百、上千种反应,流动方程组与化学反应动力学模型耦合求解形成巨大的计算量,尤其是工程中需要对复杂结构中湍流非预混反应过程进行数值模拟时,庞大的计算量超出了现有计算能力.二是计算稳定性,由于化学反应过程的特殊性,不同组分的浓度和停留时间都有着数量级上较大的差别,尤其是一些中间组分,详细化学反应动力学中包含着特征时间量级为100

-10-10的不同的基元反应,造成微分方程的强刚性,1增大求解计算难度和计算的稳定收敛.因此必须对反应动力学机理进行一定程度上的简化,以达到减少计算量提高计算效率的目的.

数值模拟过程中,研究者关心是一些重要的表征燃烧过程的参数,如点火延迟时间、火焰传播速度,一些主要组分的浓度变化及其对反应过程的影响,简化的目标是采用简化机理进行计算时,在这些重要参数与详细机理误差较小的前提下,能大幅减少计算量和计算时间.根据国内外主要采用的几种反应机理简化方法的研究现状,2简化方法主要有三类,第一类是骨架机理简化,删除冗余组分和反应来减少机理的大小,如灵敏度分析(SA)法3,4-6和主成分分析(PCA)法.7,8灵敏度分析研究反应机理中反应参数变化对整体计算结果会有影响,它需要迭代求解微分方程得到每个组分对应于其它组分及反应的灵敏度系数,当机理较大时计算量非常大,需要存储的数据量多.主成份分析作为灵敏度分析的辅助手段,是通过寻找反应过程的路径,将主要的反应路径在简化机理保留,忽略次要路径.这种方法的缺点是一些对反应过程影响比较大的组分却可能不在主要路径中出现而被忽略.主成份分析法与灵敏度分析方法结合使灵敏度分析具有方向性,提高了分析效率,但是当详细机理规模较大时,仍然会出现灵敏度系数计算量过大且难以得到收敛的灵敏度系数解.Lu和Law9提出了直接关系图(DRG)法,DRG考察组分之间的直接关系,将与重要组分关系较弱的组分视为冗余组分,通过减少这些组分,保留必需组分及相关反应来简化详细机理.蒋勇和邱榕10采用关系图法对甲烷进行了简化,在原来关系图方法的基础上进行了修正,并简化得到甲烷的187个反应,27个组分的燃烧反应机理.此外还有采用奇异摄动法,它是通过删减对所有组分影响很小的反应来达到简化的效果;11-13第二类是通过数学方法寻找简化机理的途径,如准稳态假设(QSSA)法,14低维流形技术(ILDM),15奇异摄动简化(CSP)法16,17等;低维流形技术认为反应由不同流形轨道组成的,反应系统由慢流形轨道控制,从平衡或稳态解开始,通过迭代找到慢流形轨道.由于该方法需要制表求解快的亚空间域值,对存储空间有较高的要求.CSP方法基于准稳态和准平衡假设,采用数学方法有效地确定出稳态组分,稳态组分的求解通过迭代非线性方程组,这种数学变换在提高了简化精度的同时也容易带来整个反应过程的计算不稳定性.因此为了将这两类办法的优势相结合,发展了第三类办法,综合前面两种简化方法的优势来最大程度减少计算量和提高计算稳定性,如DRG与QSSA方法相结合,18先生成骨架机理再得到总包机理,最大程度保证简化机理准确性的同时,减少机理的计算量,提高计算可靠性.对详细机理在保留重要组分和反应路径的基础上进行有效简化,是进行下一步总包机理简化的重要基础,而得到最小且合理的骨架机理则对总包机理的计算量和计算稳定性有着直接的影响.本文基于特征值分析提出了一种新的骨架简化法,目的是得到能准确反映详细机理特征的骨架机理作为总包反应简化的基础.

2 计算方法

2.1 基于特征值分析获取骨架机理的理论基础

本文在奇异摄动的理论基础上发展了一种用于获取骨架机理的机理简化分析方法,对于一个组分数为N,反应数为R,包含E个元素的化学反应系统,组分守恒方程的一般形式如下:15

其中,g表示总体净反应速率的N维列向量,Sr是第r个基元反应的化学计量系数,Fr(y)是第r个基元反应的反应速率.

取线性独立的基向量组(a1,a2,…,aN),并取对应的行向量(b1,b2,…,bN)满足关系式:

根据文献19中对特征时间的定义,取雅克比矩阵的特征值倒数作为反应系统中各个模态的特征时间.

对这N个模态进行排序,得到序列

对于燃烧过程,当特征时间远小于燃烧过程中的点火和熄火时间,可认为该模态在反应过程中由于迅速消耗而对整个反应过程不造成影响.因此通过下式来确定快模态个数M:

其中tc取燃料的点火时间,ε为一个小量,为了保证快模态的值远小于反应的时间尺度,一般取为0.01.

慢反应基元指针反映了第n个单位向量在慢模态中的映射,其计算公式如下

为了得到组分重要性的排序而在基元指针的基础上计算了组分重要性指针:

通过对组分重要性排序,筛选出对反应过程影响重要的组分,通过给出误差控制阀值ζ,当In＜ζ时认为该组分只影响快模态而对整个反应过程影响不大.对于反应的筛选,采用直接去掉删除组分所参与的相关反应.将上述方法编写成CSPSK代码,在简化时可以同时给出多个误差控制阀值ζ,生成多个不同精度的简化机理,根据与详细机理的计算误差挑选出最小,计算效率最高,最准确的骨架机理.

2.2 机理验证计算数学模型

为了验证机理的准确性,采用零维均相预混模型和层流预混火焰模型来计算燃烧过程.

一维稳态层流预混火焰模型控制方程如下,连续方程:

能量方程:

组分方程:

状态方程:

式中,x为空间坐标;m为质量流率;ρ为混合气体密度;u为流体混合物速度;A为火焰通过时的横截面面积;T为温度;cp为混合气体定压比热;Yk为物种k的质量分数;Vk为物种k的扩散速度;hk为物种k的比焓;Wk为物种k的分子量;为混合物的平均分子量,p为压力;λ为混合物导热系数;R为气体普适常数;ωk为物种k的生成率;Kg为物种的个数.

3 机理简化与验证分析

3.1 甲烷GRI1.220简化骨架机理

为了验证奇异摄动骨架简化法的合理性,对甲烷燃烧过程在零维均相预混反应器中进行简化.零维均相预混燃烧模型假设组分在任意时刻均匀混合,忽略了扩散对于反应的影响,更有利于从反应动力学角度得到与详细机理相符合的简化结果.取化学计量比(φ)为1,压力为105Pa,详细机理采用GRI1.2,该机理包括32个组分和177个反应.用零维均相预混反应模型中计算得到甲烷燃烧化学计量比(φ)为1,压力105Pa,初始温度为1000 K的初始解,从初始解数据出发,通过差分方法求得雅克比矩阵,甲烷在1000 K的点火延迟时间大约在10-1量级,因此取时间常数tc=0.1确定反应过程的慢模态即控制模态个数.进一步计算出特征向量及特征值,由公式(9)计算得到组分的重要性指针,对反应系统进行分析简化.

简化过程中所得到的组分重要性指针如图1所示.图中横坐标为反应过程的时间节点,纵坐标是组分的重要性指针,直接反映的是组分随反应时间的推进对整个反应过程的影响大小,从图中可以比较清晰地看到不同组分在燃烧过程中的活跃程度,在数量级上有着较大差别.图1(a)中组分的重要性指针值的最大值都在10-4以上,但是部分组分在反应的不同阶段变化较大,如CH4和CH3,在点火以及燃烧反应阶段,这两个组分非常活跃,与这两个组分相关的反应起着重要的速率控制作用,而在燃烧结束后,这两个组分由于已经基本被反应完全,相关的反应速率很低,对整个反应的进行起到的控制作用相对变弱.O2在熄火后仍然保持着其活跃性,由于不完全燃烧产物CO等的生成,使燃烧后剩余了部分O2,在高温下的平衡状态中,O2与一些中间产物继续进行着反应,所以O2在燃烧熄火后的重要性指针值仍然较高,而自由基H、O、OH等中间产物在整个反应过程中都体现出了较高的活跃性.图1 (b)列举了C2的相关组分的重要性指针值的分布,从图中可以看出C2的重要性指针值普遍低于图1(a)中的组分,除了组分C2H4、C2H6的重要性指针值大于10-5,C2H5在燃烧段大于10-5,其它组分的影响则相对较弱.取ζ=10-5,去掉了组分C2H、C2H2、C2H3、HCCO、CH2CO、HCCOH、C、CH、CH2、CH2(S)、CH2OH,得到的一个骨架简化机理,将其命名为CSPSK21,包括21个组分,83个反应.

图1 GRI1.2组分重要性指针值Fig.1 Important index of species in GRI1.2

图2 不同初始温度、压力和化学计量比(φ)下不同机理得到的点火延迟时间的对比Fig.2 Comparison of predicted ignition delay time using various mechanisms at different initial temperatures, pressures,and equivalence ratios(φ)

为了比较得到的简化机理的合理性,对GRI1.2详细机理,CSPSK21机理与DRM1921进行对比验证,其中DRM19机理是通过反应流和灵敏度分析方法相结合由GRI1.2得到的一个简化机理.先由GRI1.2经反应流分析得到24个组分、104个反应(组分中包括AR和N2)的机理DRM22,然后又通过计算DRM22的灵敏度系数而简化得到最终的机理DRM19,该机理包括21个组分和84个反应.

取定压零维均相反应模型,求解组分方程和能量方程计算这三个机理在φ=0.5,1,2,压力为0.1、2 MPa下,温度范围为800-1800 K的点火延迟时间,结果对比如图2所示,从图中可以看出,在中高温度下,无论是CSPSK21机理还是DRM19机理其点火延迟时间在各温度下都与详细机理非常接近,但在初始温度比较低,800-900 K时,CSPSK21机理与详细机理吻合很好,但DRM19则产生的误差比较大.同时从不同压力下的点火延迟计算结果来看,在高压下,CSPSK21机理相对来说比DRM19组分机理的结果更为理想,尤其是在φ=0.5的时候.

图3中对比了几种机理计算的点火延迟时间与实验值,22CH4、O2、Ar气分别为4.8%、19.1%、76.1%,温度范围1300-1880 K.从计算结果来看,在高温段,这几个机理都能比较准确地反映实验结果.

图3 不同机理得到的甲烷、氧气、氩气混合气点火延迟时间与实验值的对比Fig.3 Comparison mixture of CH4,O2andAr ignition delay time between the calculations based on different mechanisms and experiment resultsCH4:4.8%,O2:19.1%,Ar:76.1%

除了点火延迟以外,重要组分的浓度计算也是需要关注的重点,取层流预混火焰模型,压力为105Pa,质量流率(m)为0.04 g·(cm-2·s-1),初始网格数41.详细机理与简化机理计算得到的主要反应产物及一些重要的中间反应物的结果对比如图4和图5所示.在相同网格数的情况下,CSPSK21机理计算所需中央处理器(CPU)的时间为6 s,而GRI1.2所需要的时间为11 s.计算效率提高了将近一倍.

从上图中可以看出产物与重要中间组分的浓度分布都与详细机理吻合,CSPSK21机理在低温和中温与详细机理的计算结果相比准确度非常高,但是在高温下相比于机理DRM19与详细机理的误差略有偏大,整体上在一个较宽的温度和压力范围内与详细机理有着很好的计算吻合度.对比CSPSK21和DRM19这两个机理的组分,可以发现虽然两个机理总数相同,但具体组分有所区别.DRM19中包含组分CH2和CH2(S),而没有H2O2和CH3OH.为了验证这些组分在反应点火过程中的不同作用,针对GRI1.2机理分别在初始温度为900、1400、1800 K进行灵敏度分析,并将对温度影响排名前20的反应的灵敏度系数列于图6中.

图4 甲烷、空气层流预混火焰主要组分摩尔分数Fig.4 Main species mole fractions for CH4,air premixed laminar flamep=105Pa,m=0.04 g·cm-2·s-1,φ=1;m:mass flow rate

图5 甲烷、空气层流预混火焰中间组分摩尔分数Fig.5 Intermediate species mole fraction for CH4,air premixed laminar flamep=105Pa,m=0.04 g·cm-2·s-1,φ=1

图6 不同初始温度下使用GIR1.2计算得到的灵敏度系数Fig.6 Normalized sensitivity coefficients calculated using GRI1.2 at different initial temperatures(a)T0=900 K,p=105Pa;(b)T0=1400 K,p=105Pa; (c)T0=1800 K,p=105Pa

对比这几个图可以看出,无论是在900 K温度下还是在1400 K下,羟基参与的反应CH3+O2⇌OH+CH2O,HO2+CH3⇌OH+CH3O,H+O2⇌O+OH都对温度变化有着重要影响,影响点火延迟时间,此外2CH3(+M)⇌C2H6(+M)、CH3+O2⇌O+CH3O由于反应放热或者吸热量大而对温度影响最为敏感,这些反应以及涉及的重要中间组分H、O、OH、CH3、CH2O、CH3O都包括到了CSPSK简化机理中.对比图6(a)和6(b),在较低温度下,还有反应H+O2+N2⇌HO2+N2、2OH(+M)⇌H2O2(+M)、2HO2⇌O2+H2O2、HO2+CH3⇌O2+CH4、CH3+CH3OH⇌CH3O+CH4对温度变化影响很大,而在初始温度1000 K下这些反应的影响很小,这些反应里面包括了组分H、HO2、O2、H2O2、CH3、CH4、CH3OH、CH3O、N2.除了前面的那些组分以外还有组分H2O2、CH3OH,这两种组分正是DRM19机理中所不包括的.而随着初始温度提高,在1800 K下,CH2(S)的影响逐渐增大,反应OH+CH3⇌CH2(S)+H2O、CH2(S)+O2⇌H+OH+CO、CH2(S)+O2⇌CO+H2O中都包含有组分CH2(S),同时也可以看到这20个重要反应中也包含了反应CH3+H2O2⇌HO2+ CH4.因此可以看出,H2O2是在反应过程中无论高温还是低温都是非常重要的组分.从这几个图中反应2HO2⇌O2+H2O2的灵敏度系数也可以看出来,由于组分H2O2、CH3OH对温度影响的重要性,尤其对于低温温度的影响,因此在相当宽的范围内,CSPSK程序所得到的机理计算的温度分布及点火延迟时间都具有与详细机理较高的吻合度.

这种基于特征值分析的骨架机理简化法虽然没有计算温度灵敏度系数,但是通过组分重要性指针准确找到对反应过程有重要影响的组分,避开了对灵敏度系数繁琐的迭代求解难以收敛的问题.同时也降低了反应流和灵敏度方法中对简化者经验的要求,能方便快捷地得到所需要的骨架机理.

3.2 甲烷GRI3.020简化骨架机理

为了验证方法对不同机理进行简化的有效性,采用CSPSK程序对甲烷GRI3.0机理进行简化,取不同的误差控制阈值得到不同大小的骨架机理.当误差小于4×10-4时,得到的简化机理中将不包括含N的相关反应,不能对NOx组分浓度分布进行模拟计算.为了能通过数值模拟考察几种重要的NOx组分的浓度变化,需要选择能保留这些组分,并且能准确模拟这些组分浓度变化的机理.从图6中可以看出,当需要模拟计算NOx组分浓度分布时需要选择26个组分,120个反应,或者更大的30个组分, 140个反应的骨架机理.

为了验证这两个机理的准确性,取零维均相预混反应模型,分别计算105Pa压力下,φ=1的情况下甲烷在空气中燃烧过程,与详细机理的计算结果中得到的含氮组分的浓度分布对比如图7所示,从图中可以看出,CSPSK30和CSPSK26机理能很好地预报甲烷燃烧过程中NO组分的浓度分布变化过程,具有较高的准确性,相对于26个组分的反应机理来说,CSPSK30中包含的含N组分更多.为了验证机理的火焰传播速度的计算准确性,取层流预混火焰模型,初始温度298 K,在105Pa压力下,甲烷与空气混合气体在不同的化学计量比下对火焰传播速度进行了计算,并与实验值进行对比(图8),23-25CSPSK26机理计算得到的火焰传播速度在不同的化学计量比下与详细机理及实验值取得了较为吻合的结果.在计算时间上,当φ=1,网格数为25的条件下,CSPSK26所需要的CPU时间为7 s,而GRI3.0所需要的CPU时间为24 s.

图7 甲烷、空气火焰中含氮组分摩尔分数分布Fig.7 Profiles of N-species mole fractions for CH4,air flame

除了26个组分和30个组分的机理外,对比GRI3.0得到22个组分的简化机理与GRI1.2的21个组分简化机理可以发现,这两个机理所包含的组分几乎是一致的,只是多了NO,由于GRI3.0与GRI1.2均为甲烷燃烧的机理,因此从这个角度也体现了简化程序在对于同一种物质但是不同机理时具有同样的适用性.通过设置不同的误差控制阈值可以自动生成大小不同的骨架机理,在其中挑选最优机理则避免了简化不彻底,或者简化程度过大丢失使用者关心的重要组分信息等问题.

图8 甲烷、空气在不同化学计量比下的火焰传播速度Fig.8 Laminar burning velocity of methane,air premixed flames as a function of the equivalence ratioT0=298 K,p=105Pa

4 结论

在特征值分析及奇异摄动理论的基础上建立了一种简单实用的新的骨架机理生成的方法,并编写了用于自动简化的程序包CSPSK,用该程序对甲烷的不同详细机理进行了骨架简化.

简化得到了由GRI1.2机理简化的21个组分, 83个反应的骨架机理CSPSK21,在较宽的范围内对简化机理的合理性进行了验证,并结合灵敏度系数该机理的合理性进行了分析,结果表明,该骨架机理对不同压力不同化学计量比下都有很高的适用性.

为了在更宽范围内验证本文所提出方法的合理性,对甲烷GRI3.0机理进行了简化,在不同的误差控制阀值下,得到了不同大小的机理,经过验证, CSPSK26(26个组分,120个反应)和CSPSK30(30个组分,140个反应)都能对重要的组分和参数进行模拟计算.而使用者可以根据实际需要来选取不同的机理.

在计算负荷上,每增加一个组分会给微分方程组增加一个求解变量,加重计算的负荷.本文中的方法有效减少了组分的个数,使得计算量减少,计算效率大大提高.

基于特征值骨架简化法不需要迭代计算灵敏度系数,避开了灵敏度系数计算量大,计算收敛困难的问题,也避免了筛选组分时对经验的要求.这种方法从整个反应过程的特征来考虑,将对反应进程影响非常小的组分视为冗余组分,使骨架机理中的冗余组分减少到最少,同时满足较高的计算精度,尤其是对点火延迟时间的计算.

(1) Bykov,V.;Maas,U.Proc.Combust.Inst.2007,31(1),465.

(2) Xu,X.G.;Xu,M.H.;Qiao,Y.Coal Conversion 2004,27(4),1. [徐晓光,徐明厚,乔 瑜.煤炭转化,2004,27(4),1.]

(3) Dunker,A.M.Int.J.Chem.Phys.1984,81,2385.

(4) Turanyi,T.;Berces,T.;Vajda,S.Int.J.Chem.Kinet.1989,21, 83.

(5) Turanyi,T.New J.Chem.1990,14,795.

(6) Tomlin,A.S.;Pilling,M.J.;Turanyi,T.;Merkin,J.H.; Brindley,J.Combust.Flame 1992,91(2),107.

(7) Vajda,S.;Valko,P.;Turá nyi,T.Int.J.Chem.Kinet.1985,17, 55.

(8) Gokulakrishnan,P.;Lawrence,A.D.;McLellan,P.J.; Grandmaison,E.W.Comput.Chem.Eng.2006,30,1093.

(9) Lu,T.F.;Law,C.K.Int.J.Chem.Kinet.2005,30,1333.

(10) Jiang,Y.;Qiu,R.Acta Phys.-Chim.Sin.2009,25(5),1019. [蒋 勇,邱 榕.物理化学学报,2009,25(5),1019.]

(11) Massias,A.;Diamantis,D.;Mastorakos,E.;Goussis,D.A. Combust.Flame 1999,117,685.

(12) Massis,A.;Diamantis,D.;Mastorakos,E.;Goussis,D.A. Combust.Theory Model.1999,3,233.

(13) Lu,T.F.;Ju,Y.;Law,C.K.Combust.Flame 2001,126,1445.

(14) Xiao,B.G.;Qian,W.Q.;Yang,S.H.;Le,J.L.Journal of Propulsion Technology 2006,27(2),101.[肖保国,钱炜祺,杨顺华,乐嘉陵.推进技术,2006,27(2),101.]

(15) Mass,U.;Pope,S.B.Combust.Flame 1992,88,239.

(16) Massias,A.;Diamantis,D.;Mastorakos,E.;Goussis,D.A. Combust.Flame 1999,117,685.

(17) Lam,S.H.Combust.Sci.Technol.1993,89,375.

(18) Lam,S.H.;Goussis,D.A.Int.J.Chem.Kinet.1994,26,461.

(19) Liu,J.W.;Xiong,S.W.;Ma,X.S.Journal of Propulsion Technology 2011,32(4),525.[刘建文,熊生伟,马雪松.推进技术,2011,32(4),525.]

(20) Smith,G.P.;Golden,D.M.;Frenklach,M.;Moriarty,N.W.; Eiteneer,B.;Goldenberg,M.;Bowman,C.T.;Hanson,R.K.; Song,S.;Gardiner,W.C.;Lissianski,V.V.;Qin,Z.GRI-Mech Home Page.http://www.me.berkeley.edu/grimech(accessed Nov.25,2009).

(21) Andrei,K.;Michael,F.Reduced Reaction Sets based on GRIMech 1.2.http://www.me.berkeley.edu/drm/(accessed Nov.25, 2009).

(22) Seery,D.J.;Bowman,C.T.Combust.Flame 1970,14,37.

(23) Yasuhiro,O.;Hideaki,K.JSME Int J.Ser.B 2005,48(3),603

(24) Vagelopoulos,C.M.;Egolfopoulos,F.N.;Law,C.K.Proc. Combust.Inst.1994,25,1341.

(25)Hassan,M.I.;Aung,K.T.;Faeth,G.M.Combust.Flame 1998, 115,539.

November 22,2011;Revised:March 12,2012;Published on Web:April 1,2012.

Skeletal Mechanism Generation Based on Eigenvalue Analysis Method

WEN Fei ZHONG Bei-Jing*
(School of Aerospace,Tsinghua University,Beijing 100084,P.R.China)

A new eigenvalue analysis-based method is presented for the construction of skeletal reduced mechanisms from complex chemical reaction mechanisms.A reduced mechanism of 21 species and 83 elementary reactions for methane-air combustion was generated from detailed mechanism GRI1.2.The ignition delay time,obtained for different values of equivalence ratio,initial temperature and pressure on the basis of this reduced mechanism,were compared with those based on the detailed mechanism GRI1.2,and another skeletal mechanism DRM19.The reduced mechanism agreed favorably with the detailed model,and performed more accurately than DRM19.Two reduced mechanisms,the first involving 120 reactions among 26 species,the second,140 reactions among 30 species,were also generated from GRI3.0.They were tested by means of premixed laminar flame calculations.The method very accurately predicted speed of flame propagation and key species concentration and even NO concentration distribution in methane combustion.

Chemical mechanism reduce;Skeletal mechanism;Eigenvalue analysis;Methane; Ignition delay time;Flame propagation speed

10.3866/PKU.WHXB201204012

∗Corresponding author.Email:zhongbj@mail.tsinghua.edu.cn;Tel:+86-10-62772928.

The project was supported by the National Natural Science Foundation of China(51036004).

国家自然科学基金(51036004)资助项目

O641;O642