1 引言

随着航空运输的快速发展，人们对航空发动机产生的污染排放的关注度也越来越高，例如未燃碳氢化合物（UHC），一氧化碳（CO），烟和氮氧化物（NO_x）。NO_x是NO和NO₂的统称，会破坏放平流层的臭氧，造成酸雨和光化学烟雾等。近三十年来对航空发动机的NO_x排放规则也越来越严格，降低NO_x排放成为航空发动机燃烧室设计中一个主要目标。

NO是燃烧过程主要的NO_x组成。由于NO的生成反应的活化能较高，NO的生成相比其他燃烧产物缓慢许多，因此NO会以远低于平衡值的浓度排除燃烧室，这就给NO的精确预测带来非常大的困难。此外O和OH的预测精度同样对NO的预测精度起着至关重要的作用。因此NO的生成过程不仅受到流动和输运的影响，而且主要受化学动力学因素的控制。为提高NO的模拟精度，燃烧的模拟需要考虑详细的反应动力学。然而求解详细反应动力学需要巨大的计算资源，这在目前几乎是不可能实现的，尤其在工程应用中。为应用详细反应反应动力学而发展了多种简化模拟模型。Peters提出的火焰面模型^[1]广泛应用于湍流燃烧的模拟。稳态层流火焰面模型（Steady Laminar Flamelet Model，SLFM）利用一维层流火焰面假设，将燃烧中的热力化学特性参数，包括温度、密度、组分质量分数和化学反应源项等映射到混合分数和标量耗散率上，成为这两个控制变量的函数，并结合预先设定概率密度函数（Probability Density Function，PDF），构建用于湍流燃烧的查询表。在此基础上，Pierce ^[2]发展了火焰面/进度变量模型（Flamelet/Progress Variable，FPV），考虑了不稳定燃烧分支，一定程度上能够模拟不稳定燃烧现象，得到广泛的应用^{[3, 4]}。Pitsch等^[5]则进一步发展出非稳态火焰面/进度变量模型（Unsteady Flamelet/Progress Variable，UFPV），能够有效模拟火焰的抬举和不稳定燃烧现象。而后结合火焰面方法和流形方法发展了火焰面生成流形方法（Flamelet Generated Manifolds，FGM）用于预混火焰的模拟^[6]，FGM方法进一步扩展到用于部分预混火焰的模拟^{[7, 8]}。

然而NO的特征时间与主要燃烧产物的有非常大的差别，而FPV和FGM方法通常采用主要燃烧产物作为反应进度变量，无法直接对NO进行精确预测^{[9, 10]}。因此需要对NO的预测进行特殊处理。Godel等^[11]将NO，NO₂和N₂O等含N的组分引入到绝热FGM方法的反应进度变量的定义中，扩展反应进度变量的时间尺度范围，并模拟了抬举甲烷-空气火焰的NO分布，对NO的预测精度有一定的提高，但由于这些含N组分的在反应进度变量中的比重较低，对NO预测精度的提高有限。Ihme等^[10]以及El-Asrag等^[12]基于FPV方法引入非绝热效应，考虑辐射并且对NO的化学反应源项的正逆反应速率分别处理，有效提高了NO的模拟精度，但仍有较大的误差。Ketelheun等^[13]基于绝热FGM方法对比了不同的NO反应速率的处理方式对NO模拟精度的影响，得到较为满意的结果。

本文采用采用FGM方法和NO的附加输运方程对Sydney大学CH₄/H₂钝体稳定火焰^[14]中NO的生成进行模拟，湍流燃烧相互作用采用预先设定PDF方法。NO输运方程的源项采用两种不同的处理方式，并且考虑了扩展反应进度变量对NO预测精度的影响。考虑到NO对温度的敏感性，FGM方法的查询数据库中引入焓损失作为一个控制变量形成非绝热FGM方法。辐射模型采用光学薄模型（OTM）。由于NO受反应动力学的影响较大，GRI 2.11和GRI 3.0 ^[15]两个详细反应机理用于模拟NO生成。

2 数学模型 2.1 燃烧模型

本文中燃烧模拟由预混FGM方法完成。预混FGM方法的火焰面数据库是利用一维稳态层流无拉伸预混火焰面生成，以反应进度变量C作为控制变量，本文采用方程（1）定义反应进度变量

$ C = {Y_{{\rm{CO}}}} + {Y_{{{\rm{H}}_2}}} + {Y_{{\rm{C}}{{\rm{O}}_2}}} + {Y_{{{\rm{H}}_2}{\rm{O}}}} $

(1)

为考虑不同的混合状态来求解扩散或部分预混火焰，求解不同混合分数Z下的预混火焰面，则将层流热力化学参数Φ，例如温度、密度、组分质量分数和反应源项，映射到Z和C两个控制变量上，得到

$ \mathit{\Phi } = \mathit{\Phi }\left( {Z,C} \right) $

(2)

形成层流火焰面数据库，由于C是Z的函数，为保证PDF积分时积分变量的统计无关性，对C进行归一化处理，保证了Z和y的统计无关性。

$ y = \frac{{C - {C_0}\left( Z \right)}}{{{C_{{\rm{max}}}}\left( Z \right) - {C_0}\left( Z \right)}} $

(3)

$ \mathit{\Phi } = \mathit{\Phi }\left( {Z,y} \right) $

(4)

一维稳态层流无拉伸预混火焰面的求解采用开源动力学和火焰计算软件Cantera完成。由于预混火焰面的求解只能在贫富油熄火极限内进行，则在可燃极限外则采用外插的方式由熄火极限和燃料或氧化剂一侧的参数共同确定。本文采用Bongers提出的外插方法^[16]。

由于NO的化学反应速率与温度之间存在非常强的关系，温度的模拟精度直接影响NO的预测精度。因此考虑非绝热效应对于NO的模拟具有非常重要的影响。本文在层流火焰面的求解过程中引入焓损失ζ ^[17]考虑非绝热效应对火焰面结构的影响，其定义为

$ \zeta = h - \left[ {Z{h_{\rm{F}}} + \left( {1 - Z} \right){h_{\rm{O}}}} \right] $

(5)

式中h_F和h_O分别为火焰面中燃料和氧化剂一侧的生成焓。通过改变预混火焰面方程的边界条件的焓值h求解不同焓损失下的绝热火焰面结构，形成非绝热的层流火焰面数据库，即

$ \mathit{\Phi } = \mathit{\Phi }\left( {Z,y,\zeta } \right) $

(6)

对于湍流火焰的模拟，采用预先设定概率密度函数来考虑，Z，C和ζ采用统计无关假设。混合分数采用β-PDF，反应进度变量和焓损失则采用δ-PDF。绝热和非绝热FGM查询表分别为

$ \mathit{\tilde \Phi } = \mathit{\tilde \Phi }\left( {\tilde Z,{{\tilde Z''}^2},\tilde y} \right) $

(7)

$ \mathit{\tilde \Phi } = \mathit{\tilde \Phi }\left( {\tilde Z,{{\tilde Z''}^2},\tilde y,\tilde \zeta } \right) $

(8)

由于β-PDF的特点，其积分较为困难，本文采用分段线性积分方法（Piecewise Integration Method，PIM）^[18]处理β-PDF的积分。由于没有确定的归一反应进度变量$\tilde y$的输运方程，没办法直接求解$\tilde y$，在湍流流场中则求解非归一化反应进度变量$\tilde C$，在对FGM查询表进行查询插值是以$\tilde C$代替$\tilde y$的重新插值过程完成的。$\tilde Z$，${{\tilde Z''}^2}$，$\tilde C$和$\tilde h$的输运方程分别为

$ \frac{{\partial \left( {\bar \rho {{\tilde u}_i}\tilde Z} \right)}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_i}}}\left[ {\left( {\bar \rho {D_Z} + \frac{{{\mu _{\rm{t}}}}}{{{\sigma _Z}}}} \right)\frac{{\partial \tilde Z}}{{\partial {x_i}}}} \right] $

(9)

$ \begin{array}{*{20}{l}} {\frac{{\partial \left( {\bar \rho {{\tilde u}_i}{{\tilde Z''}^2}} \right)}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_i}}}\left[ {\left( {\bar \rho {D_Z} + \frac{{{\mu _{\rm{t}}}}}{{{\sigma _Z}}}} \right)\frac{{\partial {{\tilde Z''}^2}}}{{\partial {x_i}}}} \right] + }\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{C_{g1}}\frac{{{\mu _{\rm{t}}}}}{{{\sigma _{\rm{Z}}}}}{{\left( {\frac{{\partial {{\tilde Z''}^2}}}{{\partial {x_i}}}} \right)}^2} - {C_{g2}}\frac{{\tilde \varepsilon }}{{\tilde k}}\bar \rho {{\tilde Z''}^2}} \end{array} $

(10)

$ \frac{\partial \left( \bar{\rho }{{{\tilde{u}}}_{i}}\tilde{C} \right)}{\partial {{x}_{i}}}=\frac{\partial }{\partial {{x}_{i}}}\left[ \left( \bar{\rho }{{D}_{C}}+\frac{{{\mu }_{\rm{t}}}}{{{\sigma }_{C}}} \right)\frac{\partial \tilde{C}}{\partial {{x}_{i}}} \right]+\bar{\rho }{{{\tilde{\dot{\omega }}}}_{C}} $

(11)

$ \frac{\partial \left( \bar{\rho }{{{\tilde{u}}}_{i}}\tilde{h} \right)}{\partial {{x}_{i}}}=\frac{\partial }{\partial {{x}_{i}}}\left[ \left( \bar{\rho }{{D}_{h}}+\frac{{{\mu }_{\rm{t}}}}{{{\sigma }_{h}}} \right)\frac{\partial \tilde{h}}{\partial {{x}_{i}}} \right]-{{S}_{h}} $

(12)

式中D_Z，D_C和D_h为分子扩散系数；σ_Z和σ_C为湍流Schmidt数，取值分别为0.7和0.4；σ_h为湍流Prandtl数，取值为0.7；C_g1和C_g2取2.0；${{{\tilde{\dot{\omega }}}}_{C}}$为反应进度变量的源项，由FGM查询表获取；S_h为辐射源项。

2.2 NO预测模型

NO的时间尺度与定义反应进度变量的组分的时间尺度差别较大，NO的生成是一个相对较长的过程，尤其是热力型。这表明NO并不满足稳态假设，稳态假设会明显地高估NO的生成^[9]。

NO缓慢的生成可以通过独立求解NO质量分数的输运方程来描述^[9~11]。NO质量分数的输运方程为

$ \frac{\partial \left( \bar{\rho }{{{\tilde{u}}}_{i}}{{{\tilde{Y}}}_{\rm{NO}}} \right)}{\partial {{x}_{i}}}=\frac{\partial }{\partial {{x}_{i}}}\left[ \left( \bar{\rho }{{D}_{\rm{NO}}}+\frac{{{\mu }_{\rm{t}}}}{{{\sigma }_{\rm{NO}}}} \right)\frac{\partial {{{\tilde{Y}}}_{\rm{NO}}}}{\partial {{x}_{i}}} \right]+\bar{\rho }{{{\tilde{\dot{\omega }}}}_{\rm{NO}}} $

(13)

式中${{\tilde Y}_{{\rm{NO}}}} $为湍流流场中NO质量分数；D_NO为NO的分子扩散系数；σ_NO为NO的湍流Schmidt数，取0.7；${{{\tilde{\dot{\omega }}}}_{\text{NO}}}$为NO的净反应速率，可以表示为

$ {{{\tilde{\dot{\omega }}}}_{\rm{NO}}}=\tilde{\dot{\omega }}_{\rm{NO}}^{+}-\tilde{\dot{\omega }}_{\rm{NO}}^{-} $

(14)

式中$\tilde{\dot{\omega }}_{\text{NO}}^{+}$为正反应速率；$\tilde{\dot{\omega }}_{\text{NO}}^{-}$为逆反应速率。

逆反应速率中包括来自查询表中的NO质量分数，为避免该NO质量分数预测误差对NO生成速率的影响，可采用如下反应速率表达式

$ {{{\tilde{\dot{\omega }}}}_{\rm{NO}}}=\tilde{\dot{\omega }}_{\rm{NO}}^{+}-{{{\tilde{Y}}}_{\rm{NO}}}\frac{\tilde{\dot{\omega }}_{\rm{NO}}^{-}}{\tilde{Y}_{\rm{NO}}^{\rm{FGM}}} $

(15)

式中$\tilde Y_{{\rm{NO}}}^{{\rm{FGM}}}$为由FGM查询表得到的NO质量分数。由输运方程求解得到的NO质量分数替换逆反应速率中的NO质量分数，可排除来自查询表中的NO质量分数的不精确对NO预测精度的影响。因此方程（13）则转换为

$ \begin{array}{*{35}{l}} \frac{\partial \left( \bar{\rho }{{{\tilde{u}}}_{i}}{{{\tilde{Y}}}_{\rm{NO}}} \right)}{\partial {{x}_{i}}}=\frac{\partial }{\partial {{x}_{i}}}\left[ \left( \bar{\rho }{{D}_{\rm{NO}}}+\frac{{{\mu }_{\rm{t}}}}{{{\sigma }_{\rm{NO}}}} \right)\frac{\partial {{{\tilde{Y}}}_{\rm{NO}}}}{\partial {{x}_{i}}} \right]+ \\ \quad \quad \quad \quad \quad \bar{\rho }\tilde{\dot{\omega }}_{\rm{NO}}^{+}-\bar{\rho }{{{\tilde{Y}}}_{\rm{NO}}}\frac{\tilde{\dot{\omega }}_{\rm{NO}}^{-}}{\tilde{Y}_{\rm{NO}}^{\rm{FGM}}} \\ \end{array} $

(16)

方程（13）中的净反应速率和方程（16）中的正逆反应速率由FGM查询表获得。

由于NO相对较慢的生成过程，当反应进度变量达到最大值时，NO仍然在继续变化。这会导致在反应进度变量达到最大值附加，NO质量分数会出现极大的梯度，引起较大的插值误差。Van Oijen等^[19]认为解决这个问题的方法是扩展反应进度变量的时间尺度，在反应进度变量中引入NO质量分数，并通过与详细反应机理模拟结果的对比，验证了该方法的模拟精度，这种处理方法在低排放燃烧室中得到了应用^[20]。本文定义如下的反应进度变量

$ C={{Y}_{\rm{CO}}}+{{Y}_{{{\rm{H}}_{2}}}}+{{Y}_{\rm{C}{{\rm{O}}_{2}}}}+{{Y}_{{{\rm{H}}_{2}}\rm{O}}}+\alpha {{Y}_{\rm{NO}}} $

(17)

式中加权因子α的作用是增加NO质量分数在反应进度变量中的权重。因为NO质量分数相比其他四个组分要低几个数量级，仅添加NO质量分数并不能对反应进度变量带来有效的影响，因此通过α增加NO质量分数在反应进度变量中的权重。本文研究了采用方程（17）的定义对NO预测精度的影响。

2.3 辐射模型

方程（12）中的辐射源项S_h由光学薄模型（Optically Thin Model，OTM）^{[20, 21]}求解，可以表示为

$ {{S}_{h}}=4\sigma \left( {{{\tilde{T}}}^{4}}-\tilde{T}_{\infty }^{4} \right)\sum\limits_{i}{{{p}_{i}}{{a}_{i}}} $

(18)

式中p_i和a_i分别为组分i的分压和Planck平均吸收系数。${{\tilde T}_\infty }$为环境温度。方程（18）中考虑了CO，CO₂，H₂O和CH₄四种组分的辐射效应，其Planck平均吸收系数通过温度的拟合得到。

3 试验及数值设置 3.1 试验设置

本文的研究对象是Sydney大学的CH₄/H₂钝体稳定火焰HM1，火焰构型结构示意如图 1所示。该构型包括一个圆柱形钝体，直径为D=50mm。钝体中心为直径为3.6mm的燃料喷孔，燃料是体积比1:1的CH₄/H₂混合物，燃料射流速度为118m/s。钝体周围为伴流空气，速度为40m/s。环境压力为0.1MPa，温度为300K。化学当量混合分数为0.05。

3.2 数值设置

计算域取钝体一个60°的扇形段，以燃料喷孔出口中心为原点，向下游（X）延伸6D距离，向上游延伸1D距离以保证喷嘴出口截面的流动是完全发展的，沿径向（r）延伸3D距离。网格采用C型六面体网格，轴向非均匀设置201个网格，周向均匀设置22个网格，网格总数约为25万，如图 2所示。

N-S方程采用同位网格上的SIMPLE算法求解。求解器具有二阶空间精度。湍流由Realizable k-ε湍流模型^[23]考虑。进口边界除压力采用Neumann边界条件外，其他变量采用Dirichlet边界条件。出口边界上，压力采用Dirichlet边界条件，其他变量则采用Neumann边界条件。壁面采用标准壁面函数。周期面上采用周期边界条件。

4 结果与讨论 4.1 流动和主要热力化学参数

图 3和图 4分别为六个轴向位置处轴向和径向速度沿径向的分布。从图中可以看出，对于轴向和径向速度，试验和计算结果总体上吻合较好，能够较好地捕捉回流区的位置和尺寸。在回流区下游（X/D= 18）处，由于回流区下游滞止点位置和燃料射流穿透深度模拟的误差，使得轴向和径向速度在该位置出现一定的偏差。非绝热效应对速度的分布有着较为显著的影响，主要表现在X/D=0.9下游的中心轴线附近。考虑辐射的影响能够提高速度模拟的精度。此外反应机理对速度的影响也较为明显的影响，其影响区域与非绝热效应的影响区域一致，GRI 3.0较GRI 2.11的模拟结果与试验数据更为吻合。

图 5和图 6分别给出6个轴向位置混合分数和温度沿径向的分布情况。从混合分数的对比结果来看，除X/D=1.8和2.4以外，模拟与试验吻合较好。X/D=1.8和2.4处模拟结果稍低，这与回流区下游滞止点的模拟误差有关。这在数值上是非常困难的^[9]，即使大涡模拟仍无法精确预测这些位置的混合分数^[13]。从温度分布的对比可以看出，模拟和试验吻合较好，在X/D=0.26处，模拟高估了回流区边界附近的温度，该位置靠近壁面，壁面的影响较为明显，引起温度的降低，本文的模拟并未考虑壁面换热影响造成温度高估。在X/D=0.9和1.3处的回流区边界附近低估了温度。考虑辐射的影响能在一定程度上提高模拟精度，表明HM1火焰的辐射较弱。GRI 3.0能够更为准确地预测混合分数和温度。

图 7和图 8分别给出不同轴向位置H₂O和OH质量分数的分布情况。对于H₂O质量分数，数值模拟结果与试验数据吻合较好，在X/D=0.9和1.3处，与温度分布一致，在回流区边界附近模拟低估了H₂O质量分数的值。而在X/D=1.8和2.4的轴线附近，GRI 2.11的模拟精度明显高于GRI 3.0。考虑辐射有效提高H₂O的模拟精度。从OH质量分数的对比结果来看，FGM方法能够准确地扑捉到火焰位置。除X/D=0.26外，模拟结果与试验数据吻合较好。X/D=0.26处模拟结果出现了较大的偏差，这与温度的偏差一致，都是由壁面换热的影响造成的。辐射和反应机理对OH分布的影响并不明显。

4.2 NO的预测

图 9给出了采用GRI 3.0反应机理模拟得到的六个轴向位置处NO质量分数沿径向分布。图中“table”表示NO来自FGM查询表直接查询插值得到的NO质量分数，“NET”和“PC”分别表示由方程（13）和方程（16）得到的NO质量分数。从图中可以看出，由FGM查询表得到的NO质量分数过低预测NO质量分数，而采用方程（14）定义的反应速率则过高地预测了NO质量分数。采用方程（15）定义的反应速率排除了来自FGM查询表的NO质量分数对NO反应速率的影响，提高了NO的逆反应的反应速率，进而减小净反应速率，避免由FGM查询表造成的NO质量分数的低估和方程（13）引起的NO质量分数的高估，较为准确地反映了真实的NO反应速率，因此在所有位置与试验数据都吻合较好。从图中还可以看出，考虑辐射能提高对NO质量分数的预测精度，这对采用方程（13）得到的NO质量分数的影响尤为明显。

图 10表示基于非绝热FGM方法，分别采用GRI 2.11和GRI 3.0两种CH₄详细反应机理，由方程（13）和方程（16）得到的六个轴向位置处NO质量分数沿径向分布。

从GRI 2.11和GRI 3.0模拟结果的对比来看，相比温度和其他主要组分，反应机理对NO的模拟精度有更大的影响。对方程（14）定义的反应速率，由GRI 2.11得到的NO质量分数都明显低于GRI 3.0的；对方程（15）定义的反应速率，由GRI 2.11得到NO质量分数稍低于GRI 3.0的。采用GRI 3.0和方程（16）得到的模拟结果精度最高，与试验结果吻合较好，而采用GRI 2.11和方程（16）得到的结果精度稍低，但并没有采用方程（13）得到的结果的差别大，表明数据库中的NO质量分数对NO的预测精度有着非常大影响，排除逆反应速率中的NO质量分数能够显著提高NO的预测精度，减小了对反应机理的敏感性。

图 11~13分别为基于非绝热FGM方法，采用方程（17）定义C得到的由FGM查询表、方程（13）和方程（16）得到的六个轴向位置处NO质量分数分布，反应机理采用GRI 3.0。α的取值分别为0，10，100和1000，该取值范围保证了C的单调递增。从图 11~13可以看出，在α=1000时才对NO的模拟有非常明显的影响，由于NO质量分数相比其他定义C的组分小几个数量级，过小的α使得NO质量分数在C的定义中权重过小，无法对FGM查询表带来有效的影响，而过大的α则会造成其他定义C的组分在C的定义中权重过小，影响主要热力化学参数的预测精度，并且过大的α会导致某些混合分数下C单调性发生变化，根据文献[19]，1000是一个合适的值。

在C的定义中添加NO质量分数对由FGM查询表和方程（13）得到的NO质量分数有非常显著的影响，而对于由方程（16）得到的NO质量分数的影响则较小，这与方程（16）中的化学反应速率排除了来自FGM查询表中的NO的影响有关，这也说明在C的定义中加入NO质量分数主要影响FGM查询表中NO质量分数的模拟精度，而对NO的正反应速率的影响则非常弱，而对于逆反应速率的影响则是通过逆反应速率中的NO施加的，对其他组分几乎没有影响。

5 结论

（1）考虑辐射对钝体稳定火焰HM1的速度、混合分数、温度和主要组分的模拟精度有一定的提高，GRI 3.0比GRI 2.11在速度、混合分数、温度和主要组分的模拟上具有更高的精度。

（2）考虑辐射可有效提高NO的预测精度，GRI 3.0比GRI 2.11在NO的预测上具有更高的精度。

（3）采用输运方程求解的NO质量分数替换逆反应速率中的NO质量分数可实现高精度的NO预测，而直接采用FGM查询表中的反应速率会高估NO，特别是在高温区对NO质量分数的预测至少会高出一个数量级。

（4）在反应进度变量添加NO，只有在NO的加权因子达到1000时，能够适当提高由FGM查询表和直接采用FGM查询表中的反应速率得到的NO预测结果的精度，但仍有较大的偏差，对由NO输运方程求解得到的NO质量分数替换逆反应速率中的NO质量分数得到的NO预测结果的精度并没有明显的影响。