基于网格化Damköhler 数的燃烧室贫熄边界预测方法探究

王中豪 ，张军华 ，邓爱明，胡 斌 ，赵庆军,

(1.中国科学院工程热物理研究所，北京 100190；2.中国科学院大学航空宇航学院，北京 100190；3.中国科学院轻型动力重点实验室，北京 100190)

因贫油熄火导致的飞行安全问题一直被研究者所重视，在航发燃烧室设计阶段，需要贫熄性能评估和燃烧室结构改进的多次循环，而发展一套能够高效、高精度预测航发燃烧室贫熄边界的方法，对于缩短燃烧室设计周期，降低设计成本具有重要意义[1].

从早期Longwell 等人提出均匀搅拌反应器模型[2]以及Ballal 等人提出回流区点火模型开始[3]，基于能量平衡和时间平衡的半经验模型就在大量贫熄实验的基础上不断被发展完善[4-9].随着计算机技术的发展，大涡模拟方法被广泛应用于分析熄火前的火焰行为，部分研究者将大涡模拟方法直接应用于燃烧室的贫熄边界预测.通过大涡模拟方法对燃烧室贫熄边界进行预测，一般是通过再现贫油熄火的瞬态过程实现的，而基于不同的燃烧类型及运行工况恰当地选择燃烧模型，对于能否直接通过数值模拟准确再现贫熄过程有决定性作用.

Kim 等[10]应用涡破碎模型针对钝体稳燃预混火焰进行大涡模拟，在当量比0.52 时模拟出熄火，但是在油气比为化学恰当比时模拟的钝体后回流区长度却与实验有一定偏差，Kim 等将此偏差归因于该工况下发生的强烈的非平衡化学反应.Gokulakrishnan等[11]采用涡耗散概念模型和层流燃烧模型对钝体稳燃预混火焰进行了大涡模拟，在当量比为0.45 时，采用涡耗散概念模型能够模拟出全局熄火，与实验观察到的结果基本一致，而采用层流燃烧模型的模拟中燃烧依然维持稳定，与实验结果不符.Gokulakrishnan等人认为这是由于涡耗散概念模型假定燃烧发生在微尺度涡团中，而层流燃烧模型在亚格子尺度下假定标量不发生脉动，前者更好地模拟了湍流和化学反应的相互作用.Black 等[12]采用概率密度输运模型对低排放旋流雾化燃烧进行了大涡模拟.在实验获得的贫熄当量比下，数值模拟结果中燃烧依然稳定，未能获得预期的贫熄结果.Black 等认为是流场中局部富油而产生的高温区影响了预测结果的准确性，强调了精确模拟燃油轨迹及蒸发对于雾化燃烧贫熄预测的重要性.此外还有Garmory 等[13]和Zhang 等[14]采用条件矩方法对Sandia D&F 火焰进行的大涡模拟.模拟能够再现近贫熄工况频繁发生的局部熄火和再点火过程，Zhang 等[14]通过数值模拟获得了不同来流流速下的贫熄边界，与实验结果误差在25%～60%.Menon 等[15]采用线性涡模型对旋流预混火焰进行了大涡模拟.模拟结果展示出了火焰拉伸对火焰传播的影响，但是由于一些尺度的涡没有准确求解，未能捕捉到准确的贫熄边界.

在当前计算机发展水平下，尽管能够对燃烧室中的流场进行大涡模拟，其计算成本依然很高，并不十分适用于工程计算.在Black 等人的工作中[12]，尽管只采用了一步总包反应，每一个工况仍然需要在16 900 MHz 的计算机上计算20 d.考虑到大涡模拟仍然需要较高昂的计算成本，将贫油熄火半经验模型与雷诺平均数方法结合，直接获得不同工况下的贫油熄火边界，不失为一种高效的贫熄预测方法，也有部分研究者在此方向上做出了一些尝试.

数值模拟与半经验公式相结合的方法主要可分3 类.第一类是将燃烧室内的流场划分为不同的子区域，将流场简化为由不同反应器联接成的网络，并对每个反应器求解详细的化学反应过程，从而获得熄火边界.比较典型的有 Sturgess，Mongia 等人的工作.在Sturgess 等[16]的工作中，在基于整个燃烧室内部流动建立的反应器网络中定义一个关键反应器，当此反应器温度低于1 480 K 时认为发生贫熄.Mongia等[17-18]则是在反应器网络中的每个反应器中应用Lefebvre 的贫熄预测公式[6]，并对每一反应器中的计算结果进行体积加权，判断是否发生全局熄火.第二类方法将局部熄火模型与数值模拟相结合，首先预测局部熄火，再通过局部熄火发展的程度判断是否发生全局熄火.比较典型的是Knaus 等的工作[19].Knaus等通过Damköhler(Da)数场判断局部熄火，当钝体下游的局部熄火区域连成整体时，判断发生全局熄火.此方法的局限是对局部熄火如何引发全局熄火的解释比较粗略，没能明确指出流场中哪些区域是控制全局贫熄的关键区域.第三类方法是在数值模拟的流场中先定义一个关键反应区，将该区域作为控制熄火的关键区域，并在该反应区中应用半经验模型进行贫熄预测.比较典型的工作有胡斌等人基于燃油浓度分布定义燃烧体积并结合优化后的Lefebvre 半经验模型发展的贫熄预测方法[20-22]；Zheng 等[23]基于速度分布定义特征面结合特征面温度变化发展的贫熄预测方法；王中豪等[24-25]发展了基于燃烧室内OH基分布结合Da 数的贫熄预测方法.本研究发展了一种基于网格化Da 数的贫熄预测方法，将基于OH 基浓度定义的反应区均匀划分为尺度一致的子区域，并基于子区域计算局部和整体Da 数.该方法可以降低反应区组分、温度分布不均匀对预测结果的影响；同时基于局部Da 数场所描述的反应区燃烧组织情况，有助于进一步掌握近贫熄状态下燃烧室内流动和化学反应的匹配情况，为燃烧室设计提供更加完善的理论支撑.

1 理论介绍

1.1 关键反应区

本研究中关键反应区是基于OH 基浓度定义的，当OH 基摩尔分数高于0.001 时，认为该区域为关键反应区，该边界值是通过将数值模拟获得的OH 基分布和可视化实验[26]中由高速相机获得的反应区照片定性对比获得(实验中的反应区定义为燃烧发光部分).图1(a)为实验中获得的火焰区照片，图1(b)为数值模拟中基于OH 基定义的关键反应区，将OH 摩尔分数0.001 作为边界时，二者在形态和大小上呈现出较高的相似度.

图1 实验和数值模拟中的反应区[26]Fig.1 Reaction zone in the experiment and numerical simulation[26]

基于OH 基浓度定义反应区，综合考虑了湍流、化学反应及其相互作用对贫熄的影响.OH 基在燃烧过程中产生并在燃烧结束前消耗，其存在标志着反应的剧烈进行[27].根据链式反应假设，自由基的消灭也是反应终结的原因[28].湍流导致的反应物浓度脉动对燃烧的影响，也可以在OH 浓度分布中得到体现[29].因此将OH 基集中分布的区域作为关键反应区，能够获得反映全场燃烧状态的较为客观的信息，从而对贫熄作出较为准确的预测.

1.2 反应区分区及分区数量无关性验证

根据Longwell 等[2]提出的均匀搅拌反应器理论，整个反应区可以被简化为一温度和组分浓度分布处处相等的反应器，该假设将反应区简化为0 维反应器，大大降低了在反应区求解详细的化学反应过程的计算成本，但是实际上反应区内温度与组分分布仍然存在较大的非均匀性，在整个反应区直接应用均匀搅拌反应器假设，精确度难以得到保障.本研究通过将反应区进行分割来最大程度降低反应区内温度和组分分布的不均匀性带来的影响，形成长度尺度相同的子区域，再针对每个子区域进行均匀搅拌反应器假设.分区方法如图2 所示：①统计图2(a)所示燃烧区在x、y、z 方向的坐标范围，依据该范围建立一长方体如图2(b)，其长、宽、高对应燃烧区在x、y、z 方向的跨度，此步骤相当于将不规则的燃烧区投射为规则的长方体；②对该长方体的长、宽、高进行等数目分割(例如，若长方体的每边都分割为5 段，如图2(c)，则长方体共被分割为53＝125 块)；③分割后的长方体投射回不规则燃烧区，每一小长方体对应燃烧区中的一个子区域，如图2(d).

图2 反应区分割过程示意Fig.2 Schematic of reaction zone blocking process

子区域的长度尺度越小，其内部的温度和组分分布就越趋近于均匀，但是考虑到计算效率等因素，不可能对反应区进行无限细分，需要通过分区无关性验证工作，确保分区数目能够达到计算效率和精度的最佳平衡.本研究中选取来流空气流量1.2 kg/s 的某一邻近贫熄工况(燃油流量略高于贫熄工况)作为验证工况，计算不同分区数目下的整体Da 数，从1 开始，逐渐增大分区数目，当整体Da 数不再随分区数目增大而呈趋势性变化时，认为分区数目不再影响整体Da 数的计算.图3 为分区无关性验证结果，可以看出当分区数目从1 增加至216 时，Da 数单调减小，但是减小速度开始变慢，当分区数目大于216 时，Da数不再单调变化，而是呈小幅波动.因此在本研究中认为将反应区分为216 块，能够在保障计算精度的前提下获得最高计算效率.

图3 整体Da 数随分块数目的改变Fig.3 Variation in global Da number with the number of blocks

1.3 时间尺度

在本研究中Da 数被用于贫熄分析及预测.Da数定义为流动时间尺度(τf)与化学时间尺度(τc)的比值.该模型本身没有修正因子，因此具有在宽工况下进行贫熄预测的潜力.

1.3.1 流动时间尺度

流动时间尺度定义为新鲜混气在关键反应区的实际停留时间[24-25].该时间通过对数值模拟中的离散项颗粒进行粒子追踪获得，其过程如图4 所示.考虑到实际的燃油颗粒会蒸发为气态，设置了一个虚拟喷嘴在燃油注射点喷射低密度惰性粒子，流量在10-5kg/s 量级，以保障惰性粒子射流对燃烧不产生影响.通过追踪粒子到达和流出某区域(反应区或其子区域)的对应时刻，确定其在该区域的实际停留时间.数值模拟中惰性粒子轨迹数量为100 条，根据每条轨迹获得的流动时间结果，剔除明显的问题数据之后取平均值，以消除粒子运动轨迹随机性带来的影响.

图4 流动时间尺度获得过程示意Fig.4 Schematic of obtaining the flow time scale

1.3.2 化学时间尺度

化学时间尺度定义为新鲜来流混气能够触发目标区域内反应的最短停留时间[24-25]，通过商业软件包CANTERA 中的均匀搅拌反应器模块计算.化学时间尺度的求解方法具体为：

(1) 将均匀搅拌器初始温度设置为子区域体平均温度，以模拟高温燃烧产物对新鲜来流混气的点燃作用.

(2) 反应器的来流当量比由子区域内的混合分数反推获得，其温度等于燃烧室的入口温度，模拟新鲜混气从冷态到被加热至反应温度的过程.

(3) 反应器的压力设置为子区域内的体平均压力.

(4) 先设置较大的反应物停留时间，确保能够成功点火，再逐步减小停留时间直至点火失败，此时的停留时间即化学时间尺度.

采用瞬态求解器时，不同停留时间下的燃烧状态如图5 所示.图5(a)为稳定燃烧状态，反应器内温度维持在比较高的温度(该算例下为1 500 K 左右)，持续点燃新鲜混气；图5(b)为熄火状态，反应器内温度由初始温度降低到来流温度(495 K)，表明反应器内发生全局熄火.

图5 点火成功和点火失败时反应器温度随时间的变化Fig.5 Variations in reactor temperature with time under successful and failing ignition conditions

2 实验设置

针对扇形燃烧室模型开展了贫熄性能实验，该实验在中科院轻型动力重点实验室进行.该扇形燃烧室模型含三头部，为环形燃烧室的1/4，每头部装有双径向旋流器，第一级旋流数为1.054 6，第二级旋流数为0.947 1.火焰筒壁面采用气膜冷却，火焰筒两侧壁各有220 个冷却孔，防止侧壁因超温而被破坏，由于从冷却孔流失的空气并未实际参与燃烧，因此在计算熄火油气比时被排除.模型燃烧室设计图及实验件照片如图6 所示，图6(a)为模型燃烧室设计图，图6(b)为模型燃烧室实验件照片.燃料采用国产RP-3航空煤油，其物性参数如表1 所示.

图6 实验件设计图及实验现场照片Fig.6 Schematic and photo of test rig

表1 国产RP-3航空煤油物性Tab.1 Physical properties of Chinese RP-3 aviation kerosene

在不同空气来流流量下获得了模型燃烧室的贫熄边界，不同工况的来流温度均为495 K，来流压力为0.23 MPa，不同来流流量为0.6 kg/s、0.7 kg/s、0.8 kg/s、0.9 kg/s、1.0 kg/s、1.1 kg/s 和1.2 kg/s.在每一运行工况下，先提供充足的燃油流量维持稳定的燃烧，然后以微小调节量梯度减少燃油供给，在每次调节燃油流量后先维持该状态2 min，再进行下一次调节，防止因燃油流量突变引起的扰动影响获得贫熄边界的真实数值.通过火焰筒出口下游处热电偶监控出口温度，当温度突然降低且稳定在燃烧室入口温度(T3)附近，认为发生全局熄火.实验反复进行2 次，以确保获得的贫熄数据的可靠性.

3 数值模拟方法

3.1 湍流及燃烧模型

本研究通过商业软件FLUENT 进行数值模拟，采用可实现的k-ε 方法模拟湍流流场，采用小火焰模型对燃烧进行模化，燃烧机理采用Kundu 等[30]的煤油燃烧16 组分23 步反应机理.该机理简化了煤油燃烧时的热解过程，摒弃了燃烧过程中的一些冗余组分和反应，提升了计算效率，在当量比低于1.4 的燃烧反应模拟中表现尤为突出.该机理也在其他研究中被广泛用于组分场和温度场的获得，其可靠性得到了多次验证[31-32].对于湍流k-ε 方法舍弃了部分雷诺应力的数学限制，使之更符合湍流流动的真实物理现象，在模拟包含旋流和回流的流动时表现良好.小火焰模型将湍流火焰简化为镶嵌在湍流流场中的多个层流对冲火焰的集合，将流场内组分、温度的分布作为混合分数、混合分数脉动以及标量耗散率的函数.由于采用查表插值的方法进行化学源项相关计算，小火焰模型耦合详细机理进行化学计算的效率很高.数值方法验证工作可参阅文献[24]，限于篇幅此处不再赘述.

3.2 计算域及网格

采用循环边界对单头部燃烧室进行模拟，计算结果反映在全环燃烧室中每一头部的燃烧情况.计算域如图7 所示，机匣等与燃烧室内部流场无关的结构被简化.

图7 计算域及来流条件Fig.7 Computational domain and inlet conditions

网格划分情况如图8.考虑到该模型燃烧室结构比较复杂，除在旋流器的两级旋流通道、入口扩压段生成结构化网格，其他部分均生成非结构网格，网格总数约为350 万.

图8 模型燃烧室网格Fig.8 Mesh of model combustor

4 结果与讨论

4.1 贫熄实验结果及分析

图9 为不同来流流量下的熄火油气比.由图9可见，随来流空气流量增加，贫熄油气比几乎成线性增加，说明来流流量的增大会使燃烧室贫熄性能恶化.空气流量的增大使反应区的热损失增大，不利于燃油蒸发和燃烧.另外，空气流量增大，燃烧室内气流速度增大，导致火焰被扭曲得更为强烈，反应区更容易破碎，增大的流动速度也需要更高的火焰传播速度与之匹配，使火焰更难稳定.综上所述，贫熄油气比随来流流量的增大而增大.

图9 贫熄油气比随流量的变化Fig.9 Variation in fuel-air ratio with maunder lean blow-off conditions

4.2 局部Da 数计算结果及分析

图10 为不同工况下的Da 数场，图10(a)～(g)均为近贫熄工况，图10(h)为设计工况.图中白色虚线圈起的部分为反应区，燃烧室内主要的燃烧反应发生在该区域，红色区域为Da 数较高的区域，可以视作反应最为集中的燃烧核心.从图中可见，在近贫熄工况，尽管随着来流流量的增大，红色区域所占空间及大小均有波动，但是均集中在回流区上游滞止点附近.在回流区上游滞止点，旋流引起的回流与上游来流对撞发生滞止，此区域附近的空气流速相对较慢，利于火焰稳定.同时，上游来流带来新鲜混气，提供充足的反应物，高温燃烧产物通过回流传递到此处，提供触发反应需要的能量，因此该区域为反应最为集中的区域.从图10(a)～(g)可见，来流流量逐渐增大，燃烧区逐渐扩大，流量大于1 kg/s 后，其近燃烧区尾端表面逐渐由内凹变为外凸.燃烧区形态的改变主要是由贫熄油气比的上升导致，反应区的扩张有利于增大可燃混气的停留时间从而维持燃烧稳定，但是由于图10(a)～(g)均为不同工况下的近贫熄状态，与反应区体积同时改变的还有反应区内的温度、压力、燃料分布，多种因素作用的效果相互耦合，使反应区扩大带来的有利稳燃的因素被抵消.

图10(h)显示的是设计工况下的Da 数场.从图中可以看出，红色区域占据了反应区大部，表明空间上广泛分布的剧烈的化学反应.由于燃油供给充分，反应区甚至扩展到主燃孔下游，大大增加了反应区内可燃混气的停留时间.

图10 不同工况下的Da 数场Fig.10 Da fields under different conditions

4.3 不同油气参数下Da 数场分析

图11 为来流空气流量0.6 kg/s 时，油气比从5.27 g/kg 逐渐增加至19.2 g/kg 的Da 数场图，白色虚线圈起的部分为反应区.从图11(a)和(b)可以看出随着油气比从5.27 g/kg 增加至9.91 g/kg，反应区显著扩张，形状由内凹变为外凸，但是反应区内平均Da 数没有明显增加，也并没有在局部出现高Da 数的剧烈反应，表明燃油增加的初步效果是拓展燃烧反应发生的范围，但是燃烧强度不会大幅度提升；图11(c)表明，当油气比提升至14.56 g/kg 时，反应区进一步扩大的同时开始出现部分表明高Da 数的红色区域，表明进一步增加燃油会导致燃烧区域的增大和燃烧强度的提升；当油气比增大至19.2 g/kg 时，反应区继续小幅扩张，同时高Da 数区域范围显著扩大，表明此时增油带来的主要是燃烧强度的提升.

图11 不同油气比下反应区Da 数分布变化Fig.11 Variation in Da fields in the reaction zone at different fuel-air ratios

4.4 整体Da 数计算结果及分析

图12 为不同来流工况下的时间尺度，图12(a)为整体流动时间，图12(b)为整体化学反应时间.整体流动时间尺度通过追踪惰性粒子在整个反应区的停留时间获得，整体化学反应时间为子区域化学反应时间的体积加权平均值.从图12(a)和(b)可以看出，整体流动时间和整体化学反应时间都随来流空气流量的增大而整体呈现上升趋势.整体流动时间主要决定于反应区的体积和反应区内气流的平均流速.来流空气流量的增大及燃油流量的增大，拓展了反应区的范围，增大了反应区体积，而反应区体积的增大会增大混气在燃烧区的停留时间；与此同时，来流流量的增大直接增大了燃烧室内的气流流速，又会减少混气在燃烧区的停留时间.在本研究中，随空气流量的增大，燃烧区体积的增大对流动时间尺度的影响占据主导，但是燃烧区体积的增大会因为受到诸多因素限制而渐渐减慢，气流速度却与来流流量呈一定的线性关系，因此流动时间尺度的增加渐渐减慢.化学时间尺度受到各子区域体平均温度、压力、当量比的影响.体平均温度、压力、当量比的提升，都有利于反应的进行，减小化学反应时间.随着来流流量的增大，燃烧区内整体体平均温度、当量比都呈上升趋势，但是由于压力损失的增大，体平均压力呈下降趋势，最终化学时间受到体平均压力的主导，随来流流量的增大而增大.

图12 近贫熄工况整体流动时间尺度和整体化学反应时间尺度Fig.12 Global flow time scales and chemical reaction time scales under different lean blow-off conditions

图13 为近贫熄工况整体Da 数的计算结果，整体Da 数将整个反应区的燃烧状况用一个Da 数表示，Da 数越大表示燃烧进行得越稳定.图13 所统计的不同来流流量下的Da 数，平均值在7.9 左右，并围绕该值上下轻微波动(平均波动范围在5%左右，最大波动不超过14%).这表明由整体Da 数刻画的贫熄工况对于来流流量的变化不敏感，通过整体Da数判断贫熄具有对不同来流工况的普适性.来流工况的变化会引起燃烧室内温度、组分场变化，而本文中的流动时间尺度和化学反应时间尺度均基于流场获得，因此二者可以在一定程度上反映出由于来流工况变化引起的流场变化，从而使得全局Da 数能够在不同来流流量的近贫熄工况下保持稳定.

按照最理想情况，贫熄工况临界Da 数应等于1，实际计算结果与该值有偏离.初步分析认为，这是因为在本文所涉及的预测方法中对不同状态下的雾化质量未充分区分，未能体现贫熄状态下单油路旋流喷嘴雾化质量恶化对燃烧区的影响，未来需要补充雾化相关的实验，将不同状态下的雾化质量信息整合到预测模型中.另外，需要对分块后的每个反应器入口参数进行更精确的计算.当前方法中每个反应器的入口当量比均近似等于体平均当量比，此种计算来流当量比的方法和实际情况存在一定差距，会造成一定的误差.

图13 近贫熄工况整体Da 数Fig.13 Global Da numbers under different lean blow-off conditions

5 结论

在本研究中，发展了基于网格化Da 数的贫熄预测方法，基于贫熄实验结果通过局部Da 数和整体Da 数进行了熄火机理分析.

(1) 整体流动时间尺度和化学反应时间尺度均随来流流量增大而增大，流动时间的增大主要是因为燃烧区体积随来流流量的增大而增大，而化学反应时间的增大是因为随来流流量增大，反应区内平均压力减小.

(2) 从近贫熄状态逐渐增加燃油供给，反应区向下游扩张，反应区内的反应强度也随之提升.

(3) 基于本文整体Da 数计算结果可知，对于空气流量从0.6 kg/s 增加到1.2 kg/s 的贫熄工况，整体Da 数在7.9 上下波动，平均波动范围在5%左右.该结果表明本文的预测方法可以在一定程度上捕捉到贫油熄火的关键信息，对不同的来流工况具有一定的普适性，可为未来燃烧室贫油熄火边界预测及工程设计提供理论支持.