基于释热率分析的钝体贫油熄火过程分析

2022-03-29邢竞文金捷王方

北京航空航天大学学报 2022年3期

邢竞文，金捷，2，3，王方，2，3，*

（1.北京航空航天大学能源与动力工程学院，北京 100083； 2.北京航空航天大学江西研究院，南昌 330096；3.北京航空航天大学成都航空动力创新研究院，成都 611930）

贫油熄火边界是决定燃烧室稳定工作范围的关键参数，是发动机的最小稳定工作状态［1］。可通过实验、半经验、数值仿真与半经验结合及数值仿真4种方法来实现贫油熄火边界的预测。Dawson等［2］以实验的方法对锥形钝体稳定器的甲烷-空气湍流预混火焰进行了贫油熄火边界的测量，分析了火焰形态变化，得到了该实验条件下的贫油熄火边界曲线。陶焰明等［3］使用半经验与CFD相结合的方法预测了单头部燃烧室的贫油熄火边界，预测精度在±20%以内。随着计算机及流体力学等学科的发展，数值仿真逐渐成为研究贫油熄火边界的重要手段。由于数值仿真方法可以得到熄火过程中的诸多细节，促进了学者们对熄火过程的研究。

基于数值仿真方法的熄火研究有很多，但使用的湍流燃烧模型及研究的侧重方向不尽相同。Hodize等［4］使用大涡模拟（large eddy simulation，LES）结合有限化学速率的湍流燃烧模型，对稳定和接近吹熄条件下的火焰动力学进行了分析，认为熄火是由一些大尺度的拟序结构导致的，其主要侧重于研究熄火原因。Lee等［5］利用LES-EDC（eddy dissipation concept）湍流燃烧模型研究了Volvo钝体稳定的预混湍流火焰在贫油熄火时的火焰动力学特性，分析了贫油稳定极限下不同展向大小时涡量、速度、温度及湍流度，得到了熄火过程中火焰及流场的变化。Foale等［6］使用煤油旋流燃烧器模型验证了LES-CMC（conditional moment closure）湍流燃烧模型方法捕捉湍流喷射火焰局部熄火的有效性，表明了LES-CMC能够捕获Dagaut火焰中的局部熄灭和火焰脱离，分析了不同燃料及化学反应机理下的组分及温度分布。Paul和Najm［7］以实验的方法证明，依次获得的OH和CH2O的PLIF图像逐像素乘积的分布与层流预混火焰中的释热率有很好的相关性。Kariuki等［8］以实验的方法研究了近熄火工况下轴对称钝体上湍流预混甲烷-空气火焰的释热率。Hodize等［9］证明了LES-TPDF（transport equation probability density function，输运方程概率密度函数）方法可以很好地分辨和捕捉不同工况下的燃烧变化，使用甲烷-空气DRM19（19组分84步）反应机理初步研究了单一进口速度下、不同当量比时的速度、组分、温度及释热率的变化。

对于熄火过程的研究大多是基于速度、组分及温度等变量来进行，基于释热率的研究需要进一步发展。释热率是非稳定反应流研究中最引人注意的实验观察之一［7］。实验中直接测量释热率作为流场变量是很困难的，许多实验研究都依赖于从被认为是相关的间接测量中推断出释热率。而数值模拟可以得到释热率的分布，可以更详细地得到贫油熄火过程中释热率的变化规律。释热率可表征燃烧反应，是熄火研究中的关键因素之一，对熄火过程中的释热率进行分析有助于进一步掌握熄火规律和探索准确的湍流燃烧模型。文献［9］中仅对一种进口速度下熄火过程的释热率进行了初步研究，数据不足。本文对多种进口速度下熄火过程的释热率进行了详细的分析，并得到了以释热率为依据的熄火判据。

在诸多湍流燃烧模型中，LES能够捕捉湍流的大尺度运动，TPDF能够较好地处理湍流和化学反应之间的关系，二者结合使用可以得到湍流燃烧过程中较高精度的模拟结果。该观点也在参考文献［9］中被证实。使用基于LES-TPDF方法的AECSC程序，选用甲烷15步19组分机理。首先，分别对一个冷态工况和一个热态工况进行模拟，对比结果与实验数据，进行准确性验证。然后，对3种来流速度下远离熄火、近熄火及熄火工况进行了数值模拟，分析了组分与释热率的关系和熄火过程中释热率的变化。在本文条件下，得到了一种以释热率作为熄火判据的贫油熄火边界预测方法，为后续数值仿真方法预测贫油熄火边界的研究奠定了基础。

1 数值仿真程序与模型法

1.1 AECSC程序及LES-TPDF模型

AECSC是基于LES-TPDF方法的三维亚声速湍流燃烧数值模拟程序，取得了很多可靠的成果［10-12］。目前由笔者课题组使用并发展了气相和两相版本。AECSC程序简介详见文献［13］。本文模拟使用分块结构化网格并行求解，LES使用动态Smagorinsky-Lilly亚网格模型，燃烧模型使用PDF输运模型，在小尺度混合上采用IEM模型，其中PDF输运方程采用随机场解法。采用基于SIMPLE的预测-校正步方法进行压速耦合求解。为保证LES求解精度，空间离散大都采用二阶中心差分，时间离散大部分采用Crack-Nicholson格式，随机场的维纳项采用Euler-Maruyama显式格式离散。

使用LES方法对Navier-Stokes方程进行滤波，大于网格尺度的涡直接求解，小于网格尺度的涡具有各向同性，使用亚网格模型进行模化，因此，LES方法可以很好地模拟复杂的湍流流动。TPDF模型可以耦合详细化学机理，并且在化学反应机理方面可以精确求解且不需要模化［14］。LES与TPDF两者结合可以很好地模拟湍流-化学相互作用的复杂燃烧流动问题，已有大量学者对其详细数学公式及推导进行了研究，详细内容参见文献［15-18］。

1.2 几何模型

本文选择的实验模型为Dawson等［2］在2011年对锥形钝体短火焰进行熄火研究时的模型。原始实验装置由Balachandran等［19］发展制造，针对该实验装置已有大量学者对其进行了实验研究，数据丰富，结果可靠。实验中不但有稳定火焰与接近熄火时火焰的连续拍摄图像，而且有贫油熄火边界的数据点，适合与LES-TPDF模拟得到的结果做定性或定量对比。图1为实验装置示意图［2，20］。将半顶角为45°的锥形钝体固定在一个细杆上，安装在进气道内，进气道的出口轮廓为刀口状，钝体直径d=25 mm，进气道直径D=35 mm，细杆直径为6.35 mm，进气道壁厚为1.8 mm。常温常压下，均匀的甲烷-空气预混燃气由进口流入，点火后预混燃气在钝体后方燃烧。

图1 钝体燃烧器示意图［2，20］Fig.1 Schematic diagram of bluff-body combustor［2，20］

本次模拟中的流向沿X轴方向，计算域是一个具有矩形截面扩张型的箱体。燃烧器周围设置进口速度为0.2 m/s的伴流，用以模拟环境空气的卷吸。使用商业软件ICEM对实验模型进行建模和网格划分，网格为结构化网格。文献［9］中使用冷态算例不同流向位置的平均速度和脉动速度进行网格敏感性研究和验证，测试使用的3种网格分别为56万、150万和325万，结果表明，3种网格预测的平均轴向速度分布和均方根速度分布均与实验结果吻合较好，其中150万和325万网格预测结构相近，偶尔在某一位置处，325万网格得到的结果与实验结更为接近一些。为了能在节约计算量的同时得到较为准确的结果，本文使用的总网格数为160万。计算域及网格划分示意图如图2所示，图3为图2中白色框内钝体附近局部网格放大细节图及垂直于流向平面网格图。

图2 计算域整体示意图［9］及网格划分示意图Fig.2 Schematic diagram of overall calculation domain［9］and grid division

图3 局部网格细节Fig.3 Local grid details

1.3 边界条件

计算所用的化学反应机理均为基于GRI3.0详细CH4机理简化得来的15步19组分机理［21］。各算例编号及部分边界条件设定如表1所示。其中，冷态和热态验证算例分别用A1和A2表示，远离熄火的工况使用F1、F2、F3表示，当量比φ高于熄火工况2.5%（近熄火工况）使用C1、C2、C3表示，熄火工况使用B1、B2、B3表示，φ低于熄火工况3.3%使用L1、L2、L3表示。冷态和热态验证算例基于文献［20］，由远离熄火到熄火的热态算例基于文献［2］。本次计算进口平均速度范围在19～44 m/s之间，基于管道进口水力直径的Re在3.4×104～8.0×104之间。得到贫油熄火工况及其附近工况的方法为：在远离熄火工况F1、F2和F3的基础上，保持其他条件不变，逐渐降低φ。其他边界条件设定如表2所示。

表1 算例编号及边界条件Table 1 Number of cases and conditions of boundary

表2 其他边界条件Table 2 Other boundary conditions

2 结果与讨论

2.1 计算验证

将计算结果同实验数据进行对比，可在一定程度上验证数值模拟方法及化学反应机理的准确性，为后续熄火研究的结果提供支撑。图4为冷态验证算例A1的平均流向速度与实验测量速度的对比。平均流向速度分布已使用燃烧器出口理论平均速度进行了无量纲化，无量纲化后的速度用Umean/Ubulk表示，横坐标r/d为到燃烧器轴线的无量纲化径向距离，X/d为到钝体底面处的无量纲化轴向距离。若无特别说明，后续所提及的平均量均为时间平均量，平均速度均为平均流向速度，熄火工况均指贫油熄火工况。图4中平均速度的冷态计算结果与实验数据的相对均方根误差平均在10%以内，数值模拟结果整体上和实验结果吻合良好。

图4 不同流向位置处冷态工况平均速度与实验数据对比Fig.4 Comparison of average velocity and experimental data at different flow positions under isothermal condition

图5为冷态验证算例A1的平均速度云图和等值线图。其中，剪切层已在图中标出，黑色实线表示平均速度为0 m/s的位置，其围成的区域为平均回流区，平均回流区的流向长度在黑色虚线处结束，长度约为1.12d。

图5 冷态工况平均速度云图与等值线图Fig.5 Contours of average velocity under isothermal condition

图6 不同流向位置处热态工况平均速度与实验数据对比Fig.6 Comparison of average velocity and experimental data at different flow positions under thermal condition

结合图7所示热态验证算例A2的平均速度云图与等值线图可知，算例A2下流向平均回流区在黑色虚线处结束，长度约为1.85d。相比于冷态验证算例A1，燃烧对平均速度场的大小和结构分布特征影响很小，对平均回流区的长度影响较大，因此，燃烧会使平均回流区的长度增大。

图7 热态工况平均速度云图与等值线图Fig.7 Contours of average velocity and isogram under thermal condition

图8为实验平均OH的PLIF与由热态验证算例A2计算得到的平均OH云图，将二者进行定性对比。图8（b）中，靠近钝体底面的部分区域中，平均OH质量分数含量较实验PLIF值略少，该区域不大，且周围OH质量分数含量较多，结合文献［8］中“火焰越旺盛，火焰内部OH越少”的结论，说明计算所得火焰在该区域较为旺盛，火焰内部熄火情况没有实验PLIF严重。通过计算可知，云图左上方OH含量低的区域小于实验结果，进一步表明计算所得火焰状态较实验旺盛，但是二者整体趋势相差不是很大，可认为定性一致。

图8 热态工况实验［8］PLIF和计算平均OH云图Fig.8 Contours of mean OH of experiment PLIF［8］and mean OH of numerical simulation under thermal condition

2.2 释热率与组分

释热率是单位时间、单位体积内释放热量的度量，湍流燃烧中只有发生燃烧反应的区域才会产生释热率，因此，释热率可以表征反应速率。OH和CH2O作为中间体既有生成又有消耗，与释热率有着千丝万缕的联系。为解析三者之间的相互关系，将数值模拟得到的释热率、OH自由基、CH2O和温度分布进行相关性分析。

图9为算例F1在X/d=0.8处（数据提取位置可在燃烧剧烈区域任取）横跨火焰的平均温度、释热率、OH质量分数和CH2O质量分数分布，数据提取位置为图10蓝色虚线处。由于释热率本身数值很大，在108甚至109的量级，无法与其他变量同框，现将其除以106J·m-3·s-1进行缩小和无量纲化，得到图9中黑色实线，其对应右侧纵坐标。

图9 算例F1平均温度、质量分数和无量纲化释热率分布Fig.9 Distribution of mean temperature，species mass fractions and dimensionless heat release rate for case F1

图10 算例F1平均温度云图与平均速度等值线图Fig.10 Mean temperature and average velocity contour for case F1

CH2O的产生在预热区开始，在火焰的中温区达到峰值，并在高释热率区基本完全消耗。在CH2O开始产生后，OH开始产生，并在高温火焰区域达到峰值，此时CH2O已经被完全消耗。与释热率曲线相比较，OH和CH2O分布的重叠区域恰好与释热率峰值位置相吻合，得到的规律与文献［8］一致，进一步验证了本文中湍流燃烧模型和化学反应机理的可靠性。该现象的出现表明只有存在预热区域时才有机会发生剧烈的燃烧反应并释放热量，且只有预热区域和中间产物OH同时存在时才会发生剧烈的燃烧反应。

2.3 瞬时释热率变化

为更直观地分析瞬时释热率变化，将其与瞬时温度云图结合进行讨论。远离熄火工况（F1、F2和F3）瞬态温度和释热率云图如图11所示。此时火焰旺盛，火焰形态较长，释热率主要沿内侧剪切层出现，说明燃烧反应主要发生在内侧剪切层内，该工况下回流区内的高温区由高温燃烧产物来维持。释热率平滑起皱，但没有破裂，钝体底面边缘处出现Kelvin-Helmholtz不稳定性，其产生的影响沿剪切层向下游方向传播，有时会出现漩涡状的结构，如图11中矩形框内所示。越往下游，旋涡结构越大，卷吸影响范围越大，参与反应的混合气区域也随之变大。

图11 远离熄火工况瞬时温度与释热率云图Fig.11 Contours of instantaneous temperature and heat release rate when far away from blowoff

φ高于熄火工况2.5%的近熄火工况（C1、C2和C3）的瞬时温度与释热率云图如图12所示。火焰变短，火焰内部仍保持完整，释热率在回流区下游部分横向穿过流向气流进入到回流区内。

图12 近熄火工况瞬时温度与释热率云图Fig.12 Contours of instantaneous temperature and heat release rate when close to blowoff

图13为熄火工况（B1、B2和B3）的瞬时温度与释热率云图。火焰高度破碎，其形状在空间上发生了显著变化，火焰下游的局部熄火加剧。算例B1和B3左侧剪切层的燃烧反应被破坏，算例B2左侧剪切层内燃烧反应处于即将断裂的状态，将断裂处剪切层内的释热率无向内扭曲的趋势，由此判断剪切层内释热率的破损可能是回流区内低温气体来回运动冲撞导致的。

图13 熄火工况瞬时温度与释热率云图Fig.13 Contours of instantaneous temperature and heat release rate under blowoff condition

图14为当量比低于熄火工况算例瞬时温度与释热率云图。φ低于熄火点3.3%（L1、L2和L3）时，燃烧不稳定性增强，火焰开始变得不稳定。由图可以看出，φ低于实验熄火点时，算例L1和L2沿火焰边缘的释热率出现了断裂。其相较于熄火点的火焰熄火面积更大，火焰根部遭到的破坏更严重。算例L3的火焰由于湍流脉动的增强已经基本完全熄灭。通过上述分析可以看出，熄火过程中，不同工况下的火焰形态及释热率位置具有明显的特征差异，同一工况下不同算例具有相同的规律，进一步表明了基于LES-TPDF模型的数值模拟方法可以很好地捕捉到φ变化时的火焰变化。

图14 当量比低于熄火工况瞬时温度与释热率云图Fig.14 Contours of instantaneous temperature and heat release rate when equivalent ratio is lower than blowoff condition

2.4 平均释热率变化

为得到定量的结论，对平均释热率进行探究。图15为远离熄火工况下（F1、F2和F3）平均释热率云图与平均速度等值线图，由速度等值线图可以看出剪切层所在的位置。该工况下，平均释热率主要沿内侧剪切层出现，回流区内基本没有出现释热率，说明燃烧反应主要发生在内侧剪切层，该结论与瞬时释热率一致。平均速度等值线图中黄色实线表示平均轴向速度为0 m/s的位置，由此可看出流向平均回流区长度。

图15 远离熄火工况平均释热率云图与平均速度等值线图Fig.15 Mean heat release rate and average velocity contours when far away from blowoff

图16为远离熄火工况下降低φ到高于熄火工况2.5%时的平均释热率云图与平均轴向速度等值线图。该工况下，沿剪切层内侧的平均释热率横跨过流向在X轴上发生闭合，回流区下游也出现较大的平均释热率，回流区上游仍然没有较大平均释热率出现，这表明回流区上游已经开始有冷反应物侵入，燃烧反应也开始在回流区内发生。这种规律对于该工况下的3个算例C1、C2和C3都适用。由于沿剪切层的释热率仍然是连续的，证明冷反应物是在钝体正后方的下游进入回流区的，而不是在侧面。该工况下回流区内的高温区已经开始由燃烧反应和高温燃烧产物共同维持。

图16 近熄火工况平均释热率云图与平均速度等值线图Fig.16 Mean heat release rate and average velocity contours when close to blowoff

图17为模拟实验熄火工况得到的平均释热率云图与平均速度等值线图。实验中该工况下火焰完全熄灭，本次模拟中仍有残存火焰。但不同学者对熄火的定义不同，一些学者将熄火定义为远离稳定器火焰的熄灭，此时稳定器后方仍有残存火焰，而另一些学者则认为熄火是火焰的完全消失［2］。无论使用哪一种定义，只要能有合理的熄火判据适用于多个实验熄火工况下的数值模拟结果即可。适用于本次数值模拟的熄火判据将在后文进行讨论。

图17 熄火工况平均释热率云图与平均速度等值线图Fig.17 Mean heat release rate and average velocity contours under blowoff condition

图17为熄火工况平均释热率云图与平均速度等值线图，该图表明熄火工况下回流区上游和下游均有较大平均释热率出现，说明该工况下回流区上游也发生明显的燃烧反应，未燃混合物已经进入到回流区上游，并到达了火焰根部的位置，未燃气体对火焰根部的入侵使得火焰难以继续维持，认为这种状态下就是熄火。此时，火焰根部内侧剪切层的释热率与回流区的释热率差异减小。

图18为熄火工况基础上φ降低3.3%所得到的平均释热率云图与平均速度等值线图。可以看出，全场平均释热率在回流区的分布更加均匀，火焰根部内侧剪切层的平均释热率与回流区内平均释热率的差异更小。算例L3的平均释热率云图下翼较明亮，平均释热率大于上翼，结合图14（c）可以发现，下翼所处的位置存在600 K左右的温度区，对周围未燃燃气持续加热，产生化学反应。虽然有少量化学放热反应发生，但无法维持燃烧。

图18 当量比低于熄火工况平均释热率云图与平均速度等值线Fig.18 Mean heat release rate and average velocity contours when equivalent ratio is lower than blowoff condition

2.5 熄火判据

由上述对平均释热率的分析可以发现，在实验熄火工况中，较大的平均释热率在回流区上游和下游均有出现，且与远离熄火工况、φ高于实验熄火工况2.5%和φ低于实验熄火工况3.3%的算例相比有明显特征差异，可以此为线索来探究熄火判据。仔细观察各平均释热率云图可以发现，不同工况下平均释热率明显的差异在于靠近钝体底面附近位置处（火焰根部）内侧剪切层平均释热率与回流区平均释热率之间的差距。

提取3种进口速度下12个热态算例紧邻钝体底面X/d=0.2处（即钝体后方5 mm）垂直于X轴截面的平均释热率数据，后续称该截面位置为紧邻钝体底面处，得到平均释热率分布曲线，如图19所示。图19显示了在其他条件一定的情况下，随着φ降低，回流区的平均释热率不均匀性增强。这是由于越接近熄火湍流燃烧相互作用越强烈越复杂，放热均匀性越差。

图19 12个算例X/d=0.2处平均释热率分布Fig.19 Distribution of mean heat release rate at X/d=0.2 of twelve cases

图19（a）～（c）为远离熄火工况下紧邻钝体底面处平均释热率分布曲线，观察其横坐标，平均释热率峰值在内侧剪切层处，此处回流区平均释热率相对很低，内侧剪切层的平均释热率Hs与回流区的平均释热率Hr相差2～3个数量级，说明远离熄火时紧邻钝体底面处燃烧反应主要发生在内侧剪切层；图19（d）～（l）为熄火工况及其附近工况紧邻钝体底面处，内侧剪切层的平均释热率与回流区的平均释热率相差1个数量级左右，由此证明可通过释热率进行燃烧状态的判断。

统计紧邻钝体底面处（钝体后方5 mm）二维截面中环状内侧剪切层内平均释热率峰值的平均值Hs和圆状回流区内平均释热率的平均值Hr，得到12个热态算例的Hs/Hr，令ΔH=Hs/Hr。通过比较12个算例的ΔH值来探究熄火及其附近工况ΔH值的详细细节，并由此获得熄火判据。

各算例的ΔH值如表3所示。熄火工况及当量比低于熄火工况3.3%算例的ΔH值小于4，当量比大于熄火工况2.5%算例的ΔH值在10左右，稳定火焰算例的ΔH值更大。因此，可将ΔH值作为本次计算条件下的熄火判据，结合熄火工况下ΔH值的大小，在本次计算条件下，当ΔH＜4时发生熄火。