1 引言

经过几十年的发展，大涡模拟（LES）方法已经成为研究湍流燃烧强有力的工具。使用LES方法研究燃烧的难点在于：燃烧反应的尺度通常比计算网格尺寸更小，薄层中的燃烧必须由适当的亚网格模型来处理。目前文献中的模型有很多种^[1]，火焰增厚（TF）模型是其中之一。深入分析TF模型，可以发现它具有一些与众不同的优点。首先，它不是针对某一个燃烧模式提出来的，而是一个通用模型，已被成功用于预混火焰^{[2, 3]}、部分预混火焰^[4]和非预混火焰^[5]的计算。其次该模型的计算结果对网格精度的影响不甚敏感，计算得到的平均温度、速度场甚至是声场特性均表现出一定的网格独立性^[6]；而且由于火焰很薄，其它亚格子燃烧模型为了保证LES计算稳定需要在算法中加入数值扩散，使用多个网格来求解该火焰，这可称为“数值”增厚，而TF模型采用的是人为操纵的“物理”增厚，显著降低了所需的数值扩散，消除了一些数值不确定性^[7]。大量的计算研究表明，TF模型在预测复杂湍流与燃烧相互作用方面具有良好的表现。

在TF模型中，扩散系数及Arrhenius公式中的指前常数都需要进行人工修正，从而产生能够运用LES网格计算的增厚火焰。然而对火焰的人为加厚使燃烧模式由受输运控制变为受化学反应控制，与此同时对应于该湍流流动的火焰也发生了变化。为了考虑上述变化，Colin等^[8]引入了褶皱函数E。

研究人员做了大量的工作来改善TF模型，使之能够精准的预测火焰。但为了保证计算数据的可信度，多数运用TF模型对湍流火焰进行大涡模拟的研究中还需要一定的验证工作，找出效果最佳的参数组合。

为了使利用TF模型计算火焰时可快速确定应用参数的范围，且所得结果更加精确，本研究应用LES方法和TF模型对著名的PRECCINSTA燃气轮机模型燃烧器进行计算。对比实验与TF模型计算结果，同时对比分析该模型在不同褶皱函数E，不同火焰厚度网格数n，有无动态加厚过程等一系列条件下的计算结果：瞬时反应速率、速度、温度以及组分等数据，以开展TF模型的参数敏感性分析。最终得到利用该模型计算火焰时最佳的参数组合。

2 数值方法

在TF模型中引入褶皱函数E和增厚系数F之后，第k种组分质量分数的输运方程可表示成如下形式

$ \frac{{\partial {\rho ^{{\rm{th}}}}Y_k^{{\rm{th}}}}}{{\partial t}} + \frac{{\partial {\rho ^{{\rm{th}}}}Y_k^{{\rm{th}}}u_j^{{\rm{th}}}}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left( {{\rho ^{{\rm{th}}}}{D_k}EF\frac{{\partial Y_k^{{\rm{th}}}}}{{\partial {x_j}}}} \right) + \frac{E}{F}{\dot \omega _k}\left( {Y_k^{{\rm{th}}}, {T^{{\rm{th}}}}} \right) $

(1)

式中D，${\dot \omega _k} $，Y和T分别为扩散系数，生成速率，组分浓度和温度，上标‘th’表示加厚火焰中的物理量。通常，增厚系数F的值与网格尺寸Δ_x通过下式相联系

$ F\left( n \right) = n{\Delta _x}/\delta _{\rm{L}}^0 $

(2)

式中δ_L⁰为层流火焰的火焰厚度。通过式（2）进行加厚的火焰就可以由n个网格来进行计算，通常系数n的值在5~10 ^[9]。

文献[10]中提出了一种方法来动态控制F值，使得F能够在远离火焰处降低其值为1，亦即加厚效应完全取消，从而在此区域提供一个正确的扩散系数。本研究运用的动态过程建立在反应式参见式（10））的反应速率Q_a的基础之上。F值表示如下

$ {F_{{\rm{dyn}}}} = 1 + \left( {F\left( n \right)-1} \right)\tan h\left( {10\frac{{{Q_{\rm{a}}}}}{{{Q_{{\rm{a}}, \max }}}}} \right) $

(3)

式中Q_{a, max}为模拟一维层流火焰时Q_a的最大值^[11]。在此动态过程中，F和E是依赖于时间和空间的量，在每一个时间步得以重新计算；而在非动态过程中，仅有E依赖于时间和空间，F则是通过式（2）计算得到并且在之后的计算过程中保持不变。

Colin等基于对火焰与涡旋相互作用的直接数值模拟，引入一个亚格子尺度褶皱函数E用来考虑亚格子尺度模型（SGS）的增厚效果^[8]。该褶皱函数反映了当地过滤尺度Δ，当地亚格子尺度湍流速度u′_Δ，层流火焰流速s_L⁰，层流火焰厚度δ_L⁰和加厚火焰厚度δ_L¹ = Fδ_L⁰间的关系。函数E表示如下

$ E = \frac{{\Xi \left( {\delta _{\rm{L}}^0} \right)}}{{\Xi \left( {\delta _{\rm{L}}^1} \right)}} $

(4)

$ \Xi \left( \cdot \right) = 1 + \alpha \mathit{\Gamma} \left( {\frac{\Delta }{ \cdot }, \frac{{{{u'}_\Delta }}}{{s_{\rm{L}}^0}}} \right)\frac{{{{u'}_\Delta }}}{{s_{\rm{L}}^0}} $

(5)

$ \mathit{\Gamma} \left( {p, q} \right) = 0.75\exp \left( {\frac{{1.2}}{{{q^{0.3}}}}} \right){p^{2/3}} $

(6)

$ {u'_\Delta } = 2\Delta _x^3\left| {{\nabla ^2}\left( {\nabla \times \bar u} \right)} \right| $

(7)

式中α为模型中的常数，大小约与Re^-1/2成正比，Ξ(∙)为给定长度尺度的褶皱因数。采用式（7）能够估算由湍流现象而不是热扩散现象引起的湍流速度。

通过假设火焰面形成与破坏效应达到平衡，Charlette等提出了一种幂率褶皱函数E^[12]

$ E = {\left( {1 + \min \left[{\frac{\Delta }{{\delta _{\rm{L}}^0}}, {\mathit{\Gamma} _\Delta }\frac{{{{u'}_\Delta }}}{{s_{\rm{L}}^0}}} \right]} \right)^\beta } $

(8)

$ {\mathit{\Gamma} _\Delta }\left( {\frac{\Delta }{{\delta _{\rm{L}}^0}}, \frac{{{{u'}_\Delta }}}{{s_{\rm{L}}^0}}, {{Re }_\Delta }} \right) = {\left[{\left( {{{\left( {f_u^{-a} + f_\Delta ^{-a}} \right)}^{-1/a}} + f_{Re }^{ - b}} \right)} \right]^{ -1/b}} $

(9)

式中Γ是一个应变函数，f_u，f_Δ，f_Re是通过渐近分析得到的三个有效的模型方程^[12]。指数a和b控制着渐近线之间相互转换时的斜率。本研究将式（8）中的指数β取为0.5。

上面所讨论的模型及后续的LES计算，均通过自己的有限体积程序LESOCC2C^{[11, 13]}来实施，在分块结构贴体曲线网格上计算低马赫数的可压缩NS方程。计算中运用简化的两步化学反应机理2sCM2^[14]来描述试验中贫燃预混甲烷-空气火焰，简化反应如下

$ {\rm{C}}{{\rm{H}}_4} + 1.5{{\rm{O}}_2} \to {\rm{CO}} + 2{{\rm{H}}_2}{\rm{O}} $

(10)

$ {\rm{CO + 0}}{\rm{.5}}{{\rm{O}}_2} \leftrightarrow {\rm{C}}{{\rm{O}}_2} $

(11)

3 计算设置

本文所研究的火焰是一个大气压下的甲烷-空气火焰。燃烧器的结构参见文献[15]。它包含底部稳压室、十二通道的旋流器、燃烧室和尾气出口各一个。甲烷从旋流通道底部的小孔进入旋流器，在进入燃烧室之前与空气进行快速混合，在燃烧室内形成了“碗”形的湍流火焰。本文研究了当量比分别为0.75和0.83的两种火焰，实验测量发现这两种火焰均没有观察到强烈的热声不稳定现象^[16]，因此本研究采用低马赫近似方法是可行的。

在入口处，压力和温度分别设置为0.1MPa和320K，在当量比ϕ=0.75时旋流通道的出口处冷态体积流速设置为U_b=20m/s，对当量比ϕ=0.83时其流速为U_b=20.2m/s。因此，在已设定的喷管出口直径D，出口速度和冷流体粘度等条件下，雷诺数为3.1 × 10⁴。这两种火焰符合Borghi-Peters图的薄反应区特性，因此TF模型适应于该两种火焰。

本文把尾气出口处设置为对流出口边界条件，而在入口处采用无湍流的均匀速度分布，并假设空气和燃料已均匀混合并达到全局当量比。在所有的实体壁面处，运用绝热边界条件和Werner-Wengle壁面函数^[15]。

本文进行了不同褶皱函数E（式（4）、（8））、不同火焰厚度网格数n（2.5，5，10）和不同加厚过程（动态、静态）等条件下的数值模拟，通过对比分析不同的LES计算结果来分析TF模型对各种条件的敏感程度。表 1展示了将要讨论的十个算例的信息。其中，第一个项目中的F75表示当量比ϕ=0.75，对应的F83表示ϕ=0.83；第二个项目中的D表示具有动态加厚过程，而S表示无动态加厚过程；第三个项目中的2，5，10分别表示火焰厚度网格数n=2.5，5，10；最后一个项目中的C表示模型运用Colin褶皱方法，而P表示运用幂率函数褶皱方法。

本文采用的分块结构网格共有网格单元380万个，由ICEM-CFD软件生成，整个计算域被分割成423块。此网格曾被用于其他研究，有关计算网络划分等详情见文献[15]。

4 结果与讨论

一般的，较强的旋转入流流场内都会包含一个中心回流区（CRZ）^[17]。如文献[15]中所述，当n=5时可以观察到饱满的心形CRZ，其底部与燃烧器中心体接触，被两个小尺寸的外部回流区（ORZ）所包围。同时流场中存在3~6个螺旋型涡核^[18]，它们以同样的旋转方式规律地围绕着燃烧器的中心轴线。

4.1 火焰厚度网格数n的影响

为了在LES的结果上比较火焰厚度网格数n的影响，本研究采用了三个不同的n值，分别为2.5，5和10。图 1定性地比较了以上三种情况的瞬时反应速率Q_a。随着n值的增加，火焰高度随之单调下降，火焰褶皱减少从而导致火焰表面积减小。证明这符合人工把火焰加得越厚，小尺度涡旋造成火焰锋面褶皱就越困难的事实。预测得到的火焰表面积随着n值下降，火焰厚度却在增加，因此，总的燃料消耗率可以保持不变。如表 1所示：在相同的系数n下，平均废气温度只有1K之差，这说明整体的燃料消耗率在不同的系数n下是相同的。

LES和实验所得的流向速度的均值和均方根值展示在图 2中。F75_D_2_C和F75_D_5_C两个算例的LES结果和实验结果整体吻合良好。在x=1.5mm处存在两个剪切层，一个位于环形射流和其周围流体中间位置，另一个位于射流和中心回流区的中间区域。往下游，剪切层强度减弱，与此同时，反应物也逐渐被消耗殆尽。喷口边缘的流向流速均方根分布的两个峰值比x=1.5mm处的实验值稍小一点（如图 2（b）所示）。然而在更为下游的区域，这两个峰值的高度和位置都被很好地捕捉到了。

从图 2中可以发现：F75_D_2_C和F75_D_5_C两个算例的计算结果的差别可以忽略不计，而算例F75_D_10_C的结果与实验结果差别明显，其预测的CRZ形状和位置均不能令人满意（x≥25 mm）。结合图 1和图 2可以发现，火焰面过度的加厚导致LES较难预测湍流和火焰的相互作用，从而导致结果欠佳。另一方面，取n=10时，会导致其扩散系数加大很多，降低了可允许的最大时间步长，从而增加总体的计算工作量。综合以上讨论，当n=10时，LES计算效果不好，本文推荐n值的取值范围在2.5~5为佳。

4.2 褶皱函数E的影响

为了比较不同褶皱函数E的影响，本节对比了采用Colin褶皱函数（式（4））和幂率褶皱函数（式（8））计算得到的结果，图 3展示三种状态（F75_S_5_C，F75_D_5_C，F75_D_5_P）的瞬时火焰锋面图像（当T= 1440K时），图中的三种火焰锋面均发生了强烈的褶皱现象。但是F75_D_5_C算例下的火焰锋面皱褶最明显。另外从动画中看出三种情况中均可观察到脱离了主反应区的可燃混合气气团，但这些气团均很快被消耗，火焰顶端很不规则且存在很多鼓包，这些鼓包甚至会脱离主火焰，在F75_D_5_C状态下尤为明显。这些鼓包会间歇性高频率的产生和消失。

图 4显示的是上述三种算例下流向速度的平均值和均方根值。其中存在动态加厚过程的两种算例，其计算结果的差别几乎可以忽略不计。然而在x≥25mm范围内，无动态加厚过程算例（F75_S_5_C）的计算结果与存在动态加厚过程算例的结果及实验值都有明显偏差。如图 4（a）所示，F75_S_5_C算例下在x=15mm和25mm处的中心回流速度远大于实验值，而CRZ在x=35mm处的径向尺寸小于实验值。

经过分析，F75_S_5_C算例中不理想的结果是火焰末端不切实际的扩散率造成的。式（1）中E，F值在带有动态加厚过程的算例中是受限的，仅在火焰面处有较大的值，而在F75_S_5_C算例中E，F值在整个流场中均很大。F75_S_5_C算例中沿燃烧器中心轴在x/D=1，2和3位置处的E，F值分别为12，20和42，而其它两个算例中其值均为1。越大的E，F值所造成的扩散作用就越强烈（包括分子尺度和亚格子尺度）。靠近火焰末端（x≥35mm），没有动态加厚过程算例（F75_S_5_C）的LES结果因强烈的扩散作用改变了CRZ的行为，同时也改变了CRZ的平均尺寸。由此可见动态加厚过程对提高TF模型的性能有非常重要的作用。

4.3 温度和CO₂组分的统计信息

对于当量比ϕ=0.83的火焰，图 5和图 6分别展示了温度和CO₂的质量分数的统计数据在四个不同轴向距离沿径向的分布。从图中可以看出LES结果和实验结果吻合良好。向下游移动，更多的反应物被消耗，因此平均温度的最小值会随着x的增大而增大。在x>30mm的空间内燃烧已经基本结束，因此径向的平均温度分布图变得更加平缓。CO₂的质量分数分布也同样如此。

图 5（b）展示出了温度的均方根分布，图中对应两个平均火焰锋面位置的两个峰值（T=1440K）由火焰锋面的波动造成。在径向分布图x=30mm处平均温度的最大值和最小值之差大约为600K，然而在对应的均方根图像中其峰值为600K，这并不能完全由火焰波动来解释。如图 3所示，瞬时火焰锋面在周向（不同方位角）差别很大。在火焰顶端能发现许多的孤立的凸起，这种现象在周向具有很强的间歇性（火焰是旋流火焰），从而导致了很大的T_rms值。

如图 5所示，LES预测在x < 15mm范围内的温度要比外部的实验值高出400K（例如在r/D>0.8，x= 10mm处）。这是由本模拟中运用的绝热壁面条件忽略向环境中的热释放造成的。燃烧器底部（x=0mm）由金属制造，在测量过程中测得其平均温度为600K。

由图 5，图 6可见，F83_D_2_C和F83_D_5_C算例的结果和实验结果吻合度较高，而另外三个算例的计算结果欠佳。另一方面，F83_D_5_P中温度和CO₂质量分数均方根值在下游处（图 5（b）、6（b）中x≥30mm处）的预测值较小是火焰锋面过度的平滑造成的。如图 3所示，从F75_D_5_P算例中看到脱落的可燃混合气气团更少，这使得火焰顶端比另两种情况更加的平滑，使得间歇性以及温度和CO₂质量分数的波动性更弱。因此幂率褶皱函数替代传统的Colin褶皱函数并没有带来明显的性能改善，这可能是由于在本文的计算中将β取为常值所导致的。如果采用动态过程来确定β ^[19]，幂率褶皱函数的性能或许能有所提高。

5 结论

本文采用火焰增厚模型对著名的PRECCINSTA燃气轮机模型燃烧室进行了LES计算，同时开展了一系列的模型参数敏感性对比计算，并与实验结果进行了对比分析，得到的主要结论如下：

（1）计算获得的速度、温度以及主要组分的统计结果与实验值吻合良好，TF模型具有良好的性能及精准度。

（2）计算表明很多因素都能影响TF模型的效果。对于所研究的火焰，较小的火焰厚度网格数（n= 2.5~5）配合以动态加厚过程，加上Colin褶皱函数得到的结果与实验结果有最好的吻合度。

（3）对于较大的火焰厚度网格数（n=10）或者没有使用动态加厚过程的算例，其预测得到的CRZ与实验结果有较明显的差别。

（4）采用固定指数的幂律褶皱函数，其计算结果相比于采用传统的Colin褶皱函数并没有明显的改善。