湍流分层射流火焰的直接数值模拟

任嘉豪，王海鸥，罗 坤，樊建人

(浙江大学能源清洁利用国家重点实验室，杭州 310027)

为了获得较高的燃烧效率，减少污染物排放，先进燃烧装置通常采用贫燃预混燃烧技术[1].然而，在实际应用中，由于混合时间有限，燃料与空气混合不均匀.火焰在不均匀混合物中传播，形成分层火焰.和预混火焰相比，分层火焰中局部放热增强，温度更高，局部火焰结构发生了明显变化.有必要对分层火焰进行更深入的研究，促进高效、清洁工业燃烧装置的设计和优化.

通过实验测量，学者们研究了不同构型的分层火焰，包括槽式燃烧器火焰[2]，旋流火焰[3-4]、V 形火 焰[5-7]、球形火焰[8]和对冲火焰[9].比如Barlow 等[5-7]首次将基于Raman 的多标量激光诊断应用于湍流分层火焰中，发现当火焰传播通过更贫燃的混合物时，由于H2等中间组分和热的优先扩散效应，反应区有更高的温度和中间组分浓度，火焰传播速度更快，火焰厚度更薄.这种火焰称为后支撑火焰.Pasquier 等[8]发现火焰传播通过湍流分层混合物时，贫燃火焰利用局部富燃区提供的过量燃料和高温燃烧气体，获得更大的传播速度；相反，化学当量下的火焰传播通过贫燃混合物时，由于热损失，传播速度下降.这些研究表明混合物分层对火焰传播速度有重要影响.

在数值模拟方面，大涡模拟(LES)仅求解较大的湍流和火焰尺度，计算成本较低.学者们通过不同当量比下自由传播预混火焰建立小火焰表，进行湍流分层火焰的LES 研究[2,10-11].然而，分层火焰中后支撑燃烧所引起的局部燃烧增强等现象通常在LES 中被忽略[12].直接数值模拟(DNS)能够求解所有的湍流和火焰尺度，是分析湍流燃烧机理、验证和发展LES燃烧模型的有效工具[13-20].Poinsot 等[13]最早使用DNS 研究湍流分层火焰，认为分层火焰中燃烧的增强是因为产生了更多火焰面积.Malkeson 等[15]研究了不同湍流强度自由传播的平面分层火焰的位移速度、曲率和切向应变率统计特性.位移速度与曲率呈负相关；在强湍流时，位移速度与切向应变率相关性较高.Richardson 等[16]发现在层流分层对冲火焰中，扩散项对应变的响应很快，而反应和标量耗散项对应变的响应取决于流动和火焰时间尺度；在湍流分层射流火焰中，局部燃烧强度取决于当量比梯度与火焰法向的排列，当产物的火焰速度比反应物的火焰速度快时，燃烧强度增强[17].Wang 等[20]研究了高Karlovitz数分层射流火焰.在主反应区后，出现了反应强度较强的第二反应区.通过分析化学反应路径，发现CO、H2等中间产物在第二反应区积累，提高了该区域的热释放.

由稳燃棒进行稳燃的湍流火焰在工业装置中较为常见，目前还鲜有相关的DNS 研究.本文使用DNS 研究了湍流分层射流火焰，射流周围为空气伴随流.射流火焰通过高温稳燃棒进行稳燃，稳燃棒后产生低速区，促进V 形火焰的生成.该火焰构型类似Barlow 等[5-7]实验中的火焰构型.本文主要内容包括数值模拟细节，火焰结构、火焰分层现象和分层火焰中位移速度统计特性的分析.

1 模拟细节

本文采用DNS 研究了三维湍流分层甲烷/空气射流火焰.DNS 构型示意图如图1(a)所示.射流的宽度H 为1.0 mm.中心射流压强为0.1 MPa，温度为800 K，平均速度ub为120 m/s.周围伴随流为纯空气，温度为800 K，速度为1 m/s.射流中心(-0.2 H＜y＜0.2H)是当量比1.2 的富燃甲烷/空气混合物，两侧(-0.5H＜y＜-0.2H 和0.2H＜y＜0.5H)是当量比0.4的贫燃甲烷/空气混合物，射流的全局当量比为0.7.在这样的条件下，等效的预混火焰速度SL=1.86 m/s，火焰厚度δL=0.35 mm，火焰时间尺度τL=0.19 ms .基于 ub和 H 的射流雷诺数 Rej＝1 400.流向x 使用流入和流出边界条件.横向y 采用无滑移等温壁面边界.展向z 采用周期性边界条件.稳燃棒的温度为1 800 K，直径为0.8H，在x-y 平面内，稳燃棒中心的坐标为(4H，-0.2H).可以看出，稳燃棒中心放置在y＝0 平面的下方，这样能够获得更加显著的混合物分层效果[6].

图1 DNS构型示意，x-y 平面温度T 与速度u 二维分布Fig.1 Schematic of the DNS configuration，and twodimensional distributions of the temperature T and velocity u in an x-y plane

通过产生各向同性的辅助湍流场，得到了约为平均速度8%的脉动速度场，平均速度叠加脉动速度场 作为射流的入口速度.湍流脉动速度u'＝9.24 m/s.湍流积分长度尺度 lt=0.74 mm，其中lt=k3/2ε；k 为湍动能；ε为湍动能耗散率.湍流时间尺度τe=0.08 ms，其中τe=lt/u ′.

计算区域为 Lx× Ly× Lz=30 H ×1 0 H ×8H .各方向均使用均匀网格.网格量为 Nx× Ny× Nz=960 × 400× 240.x 和z 方向网格精度 Δ x (Δ z)为31µm，y 方向网格精度 Δ y为 25µm .考虑到层流火焰厚度为0.35 mm，剪切层处的最小 Kolmogorov 尺度η为20µm，满足ηΔx＞0 .5.所有火焰和湍流结构都足够被DNS 网格求解.图1(b)展示了火焰充分发展时刻x-y 平面温度T 和速度u 的二维分布.稳燃棒下游产生了低速区，火焰附着在稳燃棒两侧，形成V形火焰.

模拟采用基于GRI-Mech 3.0 的简化甲烷/空气燃烧化学反应机理.简化机理包含268 步基元反应和44 个组分，其中28 个组分在DNS 网格上进行输运，16 种组分设为准稳态.简化机理的验证可参考笔者先前的研究[21].

模拟使用的DNS 程序能够求解具有化学反应的可压缩流体的控制方程.DNS 程序采用4 阶显示龙格-库塔格式进行时间积分，以及8 阶显示有限差分格式进行空间离散.该模拟在国家超级计算广州中心天河2 号超算上实现，共花费约20 万机时.

2 结果与讨论

2.1 整体火焰结构

为了分析火焰的分层情况，首先定义混合物局部当量比：

式中：ξ是Bilger 形式[22]计算得到的混合分数，在纯空气中ξ定义为0，在纯燃料中ξ定义为1；ξst为当量混合分数，对于甲烷空气混合物，ξst的值为0.055.此外，还定义了反应过程变量，过程变量的定义基于燃料的质量分数 4CHY .需要注意的是，在分层燃烧中，过程变量同样依赖于混合分数[17,21].因此，过程变量的计算公式为

图2 展示了x-y 平面当量比的分布和最大热释率位置(黑色实线处).图2 中火焰可以分为上下两个分支.在上分支焰中，上游区域(4＜x/H＜8)，反应物侧的当量比大于生成物侧的当量比，属于前支撑(front-supported)燃烧.然而在下游区域(x/H＞8)，反应物侧的当量比小于生成物侧的当量比，属于后支撑(back-supported)燃烧.需要注意的是，下文的研究对象均为上分支火焰.图3 展示了上分支火焰的热释率关于过程变量的条件平均.通过图3 可以发现，热释率最大值对应的过程变量为0.85.因此将c＝0.85视为火焰面位置，将0.01＜c＜0.75 视为预热区，将0.75＜c＜0.95 视为反应区，将0.95＜c＜0.99 视为已燃区.

图2 x-y 平面当量比φ 的二维分布Fig.2 Two-dimensional distributions of the equivalence ratioφ in an x-y plane

图3 上分支火焰的热释率关于过程变量的条件平均Fig.3 Heat release rate conditionally averaged on the progress variable of the upper branch flame

2.2 火焰分层

图4 展示了典型流向位置和不同火焰分区的当量比的概率密度函数(PDF).其中实线表示预热区(0.01＜c＜0.75)，短划线表示反应区(0.75＜c＜0.95)，虚线表示已燃区(0.95＜c＜0.99).可以看出，图示两流向位置的当量比分布范围在0.4～0.95 之内，均为贫燃状态.在x/H＝6 和x/H＝12 处，不同反应区PDF 峰值所对应的当量比差异较大，表明这两处流向位置处存在明显的混合物分层现象.在x/H＝6 处，预热区PDF 峰值对应当量比大于反应区对应当量比，后者又大于已燃区对应当量比.这说明随着反应的进行，当量比逐渐减小，为前支撑燃烧模式.与其相反，在x/H＝12 处，火焰是后支撑的，且反应区当量比分布范围较大，分层燃烧现象显著.

图4 典型流向位置和不同反应区的当量比的PDFFig.4 PDF of the equivalence ratio at typical streamwise locations and various reaction zones

在x/H 为6 和12 处，混合物具有明显的分层现象.为了定量分析这两处流向位置不同当量比对热释率的贡献，定义由热释率ωT加权的当量比φ的PDF 为

图5 展示了不同流向位置的 PωT(φ).下游区域总体当量比更低，这是由于射流对周围空气的卷吸效应.在x/H＝6 处，大约50%的热量在φ＜0 .8的区域内产生.而在x/H＝12 处，80%以上的热量在这一区域产生.相对于上游，下游更多热量在低当量比区域中产生.造成这种现象的主要原因是后支撑燃烧相对于前支撑燃烧对燃烧的增强作用.

图5 不同流向位置的 PωT(φ)Fig.5 PDF of the equivalence ratio φ weighted by the heat release rate ωT at various streamwise locations

为了研究后支撑增强火焰燃烧强度的原因，笔者分析了不同流向位置处当量比关于主要中间组分H2的质量分数Y2H和温度的双条件平均(图6).图中实线表示无应变一维层流自由传播火焰的结果，其对应当量比参考彩图.可以发现，相较于x/H＝6，x/H＝12处的H2质量分数Y2H在高温侧更高，H2在高温侧积累，这一点与Wang 等[20]的结果一致.H2的积累使得产物侧的当量比更高，是前支撑燃烧转化为后支撑燃烧的重要原因.此外，x/H＝12 处H2出现在更低温的区域，与Barlow 等在实验中观察到的现象一致[5-7]，说明存在H2优先输运的现象.H2的优先输运使得预热区积累更多的活性基团，是导致该位置燃烧强度增强的重要原因.与一维层流自由传播火焰结果相比，在同一温度和当量比条件下，x/H＝6 处的H2质量分数普遍更低，而x/H＝12 处的H2质量分数普遍更高.这主要是因为前支撑对燃烧的相对削弱作用，后支撑对燃烧的相对增强作用.

图6 不同流向位置处当量比关于H2 质量分数和温度的双条件平均Fig.6 Equivalence ratio φ conditionally averaged on the H2 mass fraction and temperature T in various streamwise locations

2.3 位移速度统计

为了进一步研究前支撑、后支撑燃烧模式及混合物分层对燃烧强度造成的影响，笔者对分层火焰中位移速度的统计特性进行研究.在分层火焰中，位移速度由式(4)计算[15]：

式中：ωc为反应速率；D 为质量扩散系数；A 表示由于反应物不均匀而对位移速度的贡献，由式(5)计算：

其中n为火焰面法向量，指向反应物，定义为

图7 展示了不同流向位置处c＝0.85、φ=0.7时位移速度的PDF.需要注意的是位移速度由当量比0.7 对应的层流火焰速度 SL归一化.这里选取φ=0.7作为条件是为了排除当量比不同对位移速度造成的影响.且由图4 可以发现，不同流向位置反应区在φ=0.7时产生交集的概率最大.由图7 可以发现，不同位置处的归一化位移速度基本上分布于-2 至2 之间.x/H＝6 时，PDF 峰值对应位移速度小于0，大部分样本的火焰位移速度为负.而x/H＝12 时，PDF 峰值对应的位移速度为正，燃烧强度增加.

图7 不同流向位置处c＝0.85、φ 0.7=时位移速度的PDFFig.7 PDF of displacement speed at c＝0.85，φ 0.7=at various streamwise locations

为了进一步分析产生图7 中位移速度统计结果的原因，图 8 展示了不同流向位置处 c＝0.85、φ=0.7时位移速度和切向应变率的联合概率密度函数和位移速度关于切向应变率的条件平均.其中黑色短划线表示位移速度关于切向应变率的条件平均.切向应变率定义为

图8 不同流向位置处c＝0.85、φ 0.7= 时位移速度和切向应变率的联合概率密度函数Fig.8 Joint PDF of the displacement speed and the tangential strain rate at c＝0.85 and φ 0.7=at various streamwise locations

其中，u为当地流体速度.由联合概率密度函数可以发现，x/H＝6 处的切向应变率范围大于x/H＝12 处的切向应变率范围.由条件平均可以发现，在上游x/H＝6 处，当切向应变率为正时，位移速度与切向应变率呈负相关.较大的正切向应变率会导致较大的负火焰位移速度.而在下游x/H＝12 处，正切向应变率造成火焰产生负位移速度的影响减弱.这一趋势与Wang 等[21]在强湍流预混射流火焰中不同流向位置所观察到的结果一致，且在本研究中更为显著.基于以上两点，x/H＝6 处产生了较大概率的负位移速度.

3 结论

(1)本文使用直接数值模拟和详细化学反应机理研究甲烷/空气湍流分层射流火焰.通过对火焰结构和火焰分层情况的分析，发现在上游x/H＝6 处火焰为前支撑的，在下游x/H＝12 处火焰为后支撑的.

(2)对于相同的当量比，后支撑火焰的热释率大于前支撑火焰的热释率.

(3)通过分析H2基团在温度空间的分布，发现H2基团在产物侧的积累是前支撑燃烧向后支撑燃烧转变的重要原因.同时，后支撑火焰燃烧强度的增强与H2基团的优先输运有关.

(4)通过分析位移速度的统计特性，发现x/H＝6处负火焰位移速度的概率较大，这是由较大的正切向应变率造成的.