带有局部熄火现象的部分预混火焰大涡模拟研究

曾海翔， 王 平， SHROTRIYA Prshnt， 姜霖松 MURUGESAN Meentchidevi

(江苏大学 a. 能源与动力工程学院; b. 能源研究院，江苏 镇江 212013)

部分预混燃烧的特点是混合气在空间和时间上具有可变的当量比，且反应发生在若干子区域，因而同时兼具扩散火焰及预混火焰的特性[1].通过部分预混燃烧技术既可以降低燃烧过程中污染物的排放，又避免了预混火焰容易发生的回火问题[2-3].

早期工作中，对部分预混火焰主要通过实验测量的方法加以研究.文献[4]在高低不同雷诺数下的甲烷部分预混对冲火焰实验中，揭示了火焰熄灭主要受大尺度的涡流控制.文献[5]研究了正常和微重力条件对部分预混层流火焰的举升特性及火焰结构的影响.文献[6]在同心锥形燃烧器中产生丙烷射流部分预混湍流火焰，研究了伴流对火焰稳定性和火焰结构的影响.文献[7]就进口速度和当量比振荡对部分预混旋流火焰的扰动进行了实验探究，结果表明通过改变喷油位置、喷油器阻抗和平均喷嘴速度，可以有效地控制燃烧室进口处的当量比扰动和进口速度之间的相位差，进而抑制或放大自激压力振荡的强度.文献[8]采用粒子图像测速法和平面激光诱导荧光(PLIF)技术，研究了稀薄部分预混湍流旋流火焰中流动与反应之间的非稳态相互作用，结果表明非定常流场的主导结构受剪切层中的进动涡核控制，并对火焰的稳定起着重要作用.

得益于计算机性能的提高和数值模拟技术的发展，对湍流火焰的数值模拟逐渐成为一种常用的手段.由于燃烧过程中的火焰锋面区域发生强烈化学反应[9]，同时伴随着流场的输运和混合，这对火焰面的建模提出了极大的挑战.文献[10]利用进度变量、混合分数和标量耗散率的联合概率密度函数(PDF)统计将非预混湍流燃烧模拟中层流火焰的概念推广到部分预混区域，并通过湍流扩散举升火焰和预混本生火焰的数值模拟证明了这一方法的广泛适用性.文献[11]将大涡模拟与条件矩封闭(CMC)方法进行耦合，通过计算湍流部分预混二甲醚射流火焰，证明了CMC方法在复杂化学条件下仍可以精确地预测湍流-化学之间的相互作用.文献[12]在小火焰控制方程中通过引入火焰指数这一概念来识别不同区域的火焰状态，建立了适用于湍流喷雾燃烧大涡模拟的部分预混小火焰模型.

文献[13]首次提出“人工增厚火焰”(ATF)模型，多年以来，世界各地的专家和学者广泛对其进行了论证和改进.文献[14]通过引入效率函数解决了增厚后的火焰褶皱问题.文献[15]考虑到需要增厚的火焰位置只发生在火焰的锋面处，对于远离火焰区的位置无需增厚处理，于是引进了动态增厚因子，这极大地拓展了增厚火焰模型的适用性.文献[16]使用动态增厚火焰(DTF)模型结合小火焰生成流形(FGM)化学制表方法对旋流预混燃烧器进行模拟验证，结果与实验值吻合良好.文献[17]为研究燃气轮机内部的燃烧不稳定性与燃料混合状态的关联，分别采用完全预混和部分预混两种假设，并利用DTF模型对PRECCINSTA燃烧器当量比为0.7和0.83两种工况进行验证；在部分预混方案中，针对部分预混火焰的特性对化学反应指前因子进行修正，进而调整火焰的传播速度和厚度.结果显示在当量比为0.7(接近贫燃极限)的工况下，采用部分预混假设方案捕捉到了自激不稳定现象和频谱，而完全预混假设方案与实验结果差距较大.

DTF模型在保留火焰层流燃烧速度的同时，将火焰面加厚到了网格可解的尺度，通过直接求解化学反应方程即可获得火焰的详细信息，避免了模型封闭的困难.理论上，DTF模型广泛适用于各种燃烧模式，然而，DTF模型对部分预混燃火焰的预测还少有报道.这是由于传统DTF模型在使用过程中需要预设增厚因子、层流火焰传播速度及层流火焰厚度3个参数.在处理扩散火焰时，可以认为着火只发生在化学恰当比的对冲区域，故上述3个参数可以按照化学恰当比的层流火焰进行设置；当处理着火界限内的某个预混火焰时，同样可以根据其具体的化学当量比来设置对应的参数.然而在部分预混火焰中，可燃物的化学当量比随时间和流场位置不断地演化，依靠单一化学当量比的建模方式已不太适合，此时需要对着火界限内的层流火焰燃烧特性进行动态修正.

另外一个问题是对部分预混火焰的火焰位置的捕捉.在之前的研究中，通过监测CO2的生成反应[18]或是基于CO2的质量分数构建进度变量[19-20]所捕捉的火焰面都不甚理想.而由文献[21]提出的基于热放率捕捉火焰面的方法表现良好，故本文所使用的DTF模型同样利用反应放热率方法来识别火焰位置.

由于常温常压下，火焰的层流燃烧速度随当量比波动较为明显，鉴于此，所使用的DTF模型通过判断当地可燃物的混合分数，进而自动调整层流火焰的燃烧速度和厚度.最后，利用动态k方程亚格子模型，结合DTF模型和甲烷多步化学反应机理，对悉尼部分预混燃烧装置下FJ200-5GP-Lr75-103接近吹熄极限的工况进行了模拟验证，下文中此工况根据其速度大小简记为U103，证明了DTF模型在计算部分预混火焰方面具有足够的适应性和准确性.

1 湍流燃烧模型

1.1 湍流模型

对三维纳维-斯托克斯方程进行Favre滤波得到的大涡模拟控制方程参考文献[21]，由滤波带来的亚网格应力项及亚网格标量通量项需要模型进行处理，其中本文计算中采用文献[22]提到的动态k方程模型对亚格子应力项加以封闭：

(1)

1.2 动态增厚火焰燃烧模型

文献[13]考察了层流火焰的燃烧速度和火焰面厚度，发现其均匀与组分的扩散系数和化学反应有某种关联：

(2)

(3)

文献[14]认为增厚火焰模型降低了火焰对流场的敏感性，于是引入褶皱因子E的概念.针对部分预混火焰的特性，首先基于GRII-Mech3.0机理得到了可燃界限内甲烷常温常压下，层流燃烧速度sL和层流火焰面厚度δL随当量比φ的变化情况，如图1所示.随后，通过MATLAB提出二者随当量比变化规律的拟合函数，并植入DTF模型中.在所使用的DTF模型中，火焰面的位置通过多步化学反应的放热量dq来捕捉，使得火焰的增厚系数只在燃烧火焰区内有效，修正后的火焰增厚因子Fdyn和褶皱因子Edyn分别由下式表示：

(4)

式中：ndtf为预设的增厚因子；Δe为网格尺度；a1、a2为调整流场中动态增厚曲线形态的缩放因子，此处a1=10,a2=25.褶皱因子可由下式表示：

(5)

(6)

图1 层流火焰的燃烧特性Fig.1 Combustion characteristics of laminar flame

2 算例设置

2.1 悉尼燃烧器

为研究入口处燃料的混合情况对火焰的影响，文献[23]精心设计了一套燃烧器，结构如图2(a)所示，其中：Lr为回缩距离.燃烧器的入流管道是由两个壁厚为0.25 mm的同心管组成，内管通入甲烷气体，内管直径为4 mm，外管通入空气，外管直径为7.5 mm.混合气体通过外围导管中的高温燃烧废气点燃，废气导管的直径为18 mm，整个装置放置于25 mm × 25 mm的伴流风洞中.燃烧器具有凹入的内管，该内管可以从外管出口平面向上游回缩.当Lr=0时，内管与外管平齐(内管位于外管出口平面)，此时混合气燃烧模式为完全非预混，即扩散燃烧模式；当Lr≥300 mm(40倍的外管直径)时，此时燃气与空气充分混合，燃烧模式为完全预混燃烧模式.当Lr处于两个极限距离之间时，燃烧模式为部分预混燃烧模式.为验证DTF模型对部分预混火焰的预测能力，计算内管回缩距离选取 75 mm，对应混合气体处于部分预混状态，速度入流参数如表1所示，其中：U为内外管混合段的混合气平均速度；UA为外管入流段空气的平均速度；UF为内管入流段燃料的平均速度.为了点燃射流火焰，在射流火焰的外围制造了一个环形的值班火焰区域，其中注入的是由C2H2、H2、CO2、9N2和空气这5种组分按照甲烷-空气火焰在化学恰当比时的元素配比来配置的混合气.计算过程中值班火焰入口取为实验测量结果，具体信息如表2所示.其中：w(N2)、w(CO2)、w(H2O)、w(O2)和w(CO)分别为对应组分的质量分数；Up为值班火焰环空内混合气的平均速度；T为值班火焰峰值温度.

图2 悉尼部分预混燃烧器及计算网格示意图[23]Fig.2 Schematic diagram of Sydney partially premixed burner and calculation grid[23]

表1 U103工况下入流速度信息Tab.1 Inflow speed information under U103 condition

表2 U103工况下值班火焰尾气中各组分质量分数、流速和温度信息

2.2 计算设置

本文计算基于开源的计算流体力学软件OpenFOAM，其向来以其代码的高效率和优良的可拓展性著称.在原有的rhoReactingFoam 求解器基础上，本课题组植入了DTF模型.计算采用一套分块的包含3.82×106个单元的结构网格，如图2(b)所示.其中，径向、周向、轴向分别部署了76×80×628个网格单元，火焰区域内的网格最大宽度为0.2 mm.

化学反应机理对模拟的精度有着明显的影响，文献[24]考量了5种总包机理与文献[25]提出的骨架机理在贫燃部分预混旋流火焰中的表现，所获得的结论为，总包机理可以大大降低计算成本，但同时也影响结果的质量；而包含73步化学反应的Lu骨架机理表现与实验值吻合的最好，但计算成本也远高于总包机理.本文中化学反应机理选用文献[26]提出的甲烷19步骨架反应机理，该机理同样是由GRI-Mech1.2机理简化而来.计算中燃料和空气入口施加强度为10%的随机扰动，时间离散为二阶Backward格式，对流项采用二阶有界线性格式，拉普拉斯项选用高斯线性正交.

3 结果与讨论

流场在不同轴向位置处的流场信息如图3所示.其中：r为轴向距离；x为径向距离；d为钝体直径；Ux为径向平均速度；Ux_rms为径向均方根速度.从图3中可以看出，可燃混合气刚进入燃烧室时，中心轴向处的速度较高，并向外逐渐衰减，混合气与周围伴流气体接触的剪切层中的速度梯度很大.随着混合气沿着下游方向不断的发展，主流的平均速度不断降低，而周围流体的速度却在不断增大.这是由于剪切层的作用，使得速度较高的可燃气不断地将周围的空气卷入主流中，以降低剪切层中的速度梯度，而周围流体由此受益使其平均速度逐渐增大.观察流场中的均方根脉动，可以看到在整个流场中，均方根脉动的峰值几乎没有发生变化，但脉动的频谱却越来越宽，整个流场越来越趋于一致.

图3 不同轴向位置处的速度场信息Fig.3 Velocity field information at different axial positions

为了辨别流场中的涡旋结构与位置，利用无量纲化的Q准则和涡量ω提取流场中的拟序结构，如图4所示.Q准则云图中可以看到在火焰的根部r/d约为0.5的位置处达到峰值，然后向下游缓缓的扩张开来，在x/d=5到x/d=20的区间内涡旋较为密集.涡量图显示了流体微团的旋转规律.从图4中可以看到，沿轴向的涡旋几乎总是与反向涡旋成对出现，在x/d>20的下游位置涡量逐渐减弱.对于不同轴向位置，其涡旋峰值所在的位置正好与图4流场中速度沿x方向的均方根波动峰位相吻合，这进一步表明由剪切层产生的涡旋加剧了流场的掺混和轴向速度的传递.

以二氧化碳的瞬时生成速率云图RCO2及火焰指数ξ来分析火焰面的信息，两种方案的结果如图5所示.其中：虚线为RCO2=3 kg/(m3·s) 的等值线，用来勾勒火焰面的位置.由图5可知，由CO2生成的速率云图中不同时刻的火焰面的底部x/d<10的部分均没有产生中断，这意味着火焰底部附近的燃烧较为稳定.这是由于外部的值班火焰(根部左、右两支角状火焰分支)持续作用的结果.由于过量的甲烷和空气在导管内来不及充分预混便进入到燃烧室，火焰根部的混合分数较高，只有在外层较为稀薄的地方被伴流火焰所点燃，所以根部附近的火焰面较窄，通过分析火焰指数云图可以看到，两种方案在此区域均为预混燃烧模式.然而在火焰的内部区域，混合气体此时由扩散模式主导.在火焰中部1020)火焰继续向下游发展，孔洞出现的越来越多，其范围也逐渐增大，局部熄火现象进一步加剧，这一区域的火焰由扩散模式所主导.

图4 燃烧室内无量纲化的Q准则及涡量ω图Fig.4 Dimensionless Q criterion and vorticity ω in combustion chamber

图5 CO2瞬时生成速率云图和火焰指数云图Fig.5 Cloud maps of CO2 instantaneous generation rate and flame index

图6 不同轴向位置处混合分数均值及均方根波动值对比Fig.6 Comparison of mean mixture fraction and root mean square fluctuation at different axial positions

可燃气的混合状态决定了火焰的温度及燃烧产物，故而计算中必须对混合分数做出精准的预测.图6显示了不同轴向位置处混合物的质量分数对比，其中：Zmean为混合分数均值；ZRMS为混合分数均方根波动值.由图6可以看到，无论是火焰的根部亦或是下游区域，采用动态k方程湍流模型及结构化的382万个网格对可燃气的混合分数的计算结果与实验值吻合良好.在整个流场中，混合气的质量分数呈现出了较好的对称性，在r/d=0的中心线上，混合气的质量分数达到最大，然后向两边逐渐递减.然而混合分数的均方根脉动峰值却没有出现在中心线上，而是在中心线附近并沿中心线对称，呈现出两个峰状.随着高度的增加，混合气的质量分数峰值逐渐降低，火焰面随之也变得越来越宽.这是由于边缘较为稀薄的可燃气不断的被消耗，混合气从浓度高的中心区域不断向两边扩散，整个混合气的浓度逐渐被稀释的结果.

图7为计算结果与实验值的温度对比，其中：Tmean为温度均值；TRMS为温度均方根波动值.整体上看，改进后的DTF模型计算结果较为良好的展现出了火焰中的温度分布及波动趋势.当x/d<10时，平均温度峰值出现在r/d=1附近，而均方根波动的峰值出现在r/d=0.5附近；随着下游熄火位点的增多，温度场的均值不断下降，而其波动也逐渐蔓延到外侧区域.

图7 不同轴向位置处的平均温度及均方根波动值对比Fig.7 Comparison of mean value and root mean square fluctuation of temperature at different axial positions

文献[2]在层流部分预混射流火焰的直接数值模拟中指出，部分预混燃烧火焰中包含两个反应层，内层中的甲烷和氧气不断消耗产生CO、H2、OH等中间产物，然后在外层的氧化层中被反应为CO2和H2O.这里选择CO作为中间产物的代表，并与实验值进行对比，结果如图8所示，其中:w(CO)mean为一氧化碳质量分数均值;w(CO)RMS为一氧化碳质量分数均方根波动值.对于中上游中间产物的计算结果，DTF模型结果较好的展现了中间产物CO的含量及波动特征.同时可以看到在火焰的下游区域，计算结果较实验值出现了一定的偏差，此区域的CO质量分数均方根波动也变得较为剧烈，这可能是由于采用简化机理无法较为全面地描述CO的反应过程.从图8左图可以看到，不同平面高度上CO的峰值均出现在r/d=0.75附近，然后向两边缓慢扩散；到了下游区域，曲线逐渐变得平滑.与图7中温度曲线的峰值位置进行对比，可以发现放热峰值出现在r/d=0.9附近，稍远于CO峰值面，这表明Suresh的假设是比较合理的.由于避免了环境因素的影响，计算值比实验值呈现出了更好的对称性.

将混合气的温度随质量分数分布制成散点图如图9所示，其中：T为温度；Z为混合分数；蓝色虚线为化学恰当比所在线，将火焰分割为左、右两支.左半支的混合分数较低，而对应的温度较高，故其为高温燃烧废气与伴流空气的混合支，右侧的混合分数较高但对应的温度比较低，对应可燃气与高温燃烧废气的混合支.由于入流速度过高(接近吹熄极限)，导致火焰根部x/d=1附近大量的可燃气体来不及引燃便被带到下游.进一步在x/d=5的位置观察火焰燃烧状态，可以发现此时混合气的燃烧状态较为理想.从x/d>10，左、右两个混合分支间开始出现熄火点，越到下游熄火点出现的越为密集.由图3中的火焰指数云图可以看出，在流场的中下游，扩散燃烧模式逐渐占据主导地位.对比实验和计算的散点图可以发现，采用改进后的动态层流火焰燃烧速度及层流火焰厚度的DTF模型整体非常良好的捕捉到了这些变化趋势.散点图中火焰的根部x/d=1处，无论是实验亦或是计算模拟，散点数据都比较少；且大部分散点都集中在化学当量比分割线附近.当混合分数Z超过0.1，300 K附近出现一条水平直线.这是由于混合气刚进入燃烧室时，底部主流区域的可燃气混合分数过高且流速过快，来不及与周围的燃烧废气进行热量交换，这部分散点全部集中在300 K的带状分布附近.由于这部分的散点比例过大，导致整个散点图中的点显得比较稀少.

图8 不同轴向位置处混合物的CO质量分数平均值及均方根波动值对比Fig.8 Comparison of mean mass fraction and root mean square fluctuation of CO at different axial positions

图9 不同轴向位置处混合物的温度随混合分数分布散点图Fig.9 Scatter plot of temperature with mixture fraction distribution at different axial positions

4 结论

为了应用动态增厚火焰模型计算当量比不均匀分布的湍流部分预混火焰，需要考虑到当量比变化对火焰褶皱函数中出现的层流火焰速度、厚度的影响.本文通过计算一系列的层流火焰构造了两个拟合函数，其可以在DTF模型中自适应调整火焰参数，很好地解决了上述问题.结合动态k方程亚格子湍流模型，对接近吹熄极限的悉尼部分预混火焰进行了大涡模拟研究.本文研究所获得的结论如下.

(1) 通过流场中混合分数均值及均方根脉动对比，显示了动态k方程亚格子湍流模型可以很好地预测非均匀入流的射流流动问题.

(2) 观察CO2生成速率云图和OH基渲染的三维瞬时火焰面，不难发现采用改进后的DTF模型对此类带有熄火现象的部分预混火焰具有较好的预测能力，温度-混合分数散点图分布也很好的说明了这一点.

(3) 本文使用的甲烷多步简化机理对某些组分如CO的质量分数预测较实验值偏低，若采用更为详细的化学反应机理或结合化学制表方法，这种现象可能会有一些改善.