扩散燃烧是指在化学反应之前，燃料与氧化剂未预先混合，燃料与氧化剂的混合和化学反应同时进行，火焰因而产生于燃料与氧化剂的介面^[1]。由于扩散火焰的燃料和氧化剂未预先混合，故在燃料的贮存上较安全，且无回火和爆炸的危险，因此在工业上的应用颇为广泛。

LES是计算成本和准确度之间的最佳折衷方法。相比于RANS，LES可以捕捉到相对详细的湍流流场中涡团的运动特征，如燃烧室的非稳态燃烧过程和燃烧不稳定性现象^[2-5]。经过近些年的发展，LES湍流燃烧模型可以概括为以下几类：唯象的方法，即根据湍流和火焰的宏观表现来简化其关系，以EDC和PaSR为代表；火焰面假设模型，即认为湍流脉动对整个火焰面内部结构无影响，以稳态火焰面，FPV和FGM为代表；另外一种就是用概率密度的方法来考虑湍流和火焰之间的相互联系，以预设和输运型PDF为代表。国内外很多学者利用这些方法模拟扩散或者预混的湍流火焰，都得到了比较满意的结果。陈义良等^[6-8]连续多年利用不同化学反应机理(详细机理GRI3.0，简化机理等)和湍流燃烧模型(火焰面，PDF等)对甲烷射流火焰进行模拟研究，包括温度场，速度场及燃烧产物分布(NO_x)。Fureby^[9-10]用大涡模拟PaSR模拟了不同类型的燃烧室，包括SGT-800燃气轮机燃烧室和超音速发动机燃烧室，并得到了不错的计算结果，为改进设计提供了详细的动态数据。Peters等^[11-12]先后提出并改进了稳态火焰面和FPV模型，理论分析和计算结果证明了该模型满足于扩散，预混和部分预混火焰。Pope等^{[3, 13-15]}使用PDF方法成功地模拟了射流和旋流火焰，内燃机及燃气轮机燃烧室的燃烧过程。另外，周力行等^[16-18]对LES湍流燃烧模拟做了详实的分析和对比。

从现有文献来看，多数学者利用同一种LES湍流燃烧模型去分析火焰理解现象规律，而对于不同大涡模拟湍流燃烧模型性能的对比研究还不多，在此基于开源计算软件OpenFOAM，用PaSR，稳态火焰面和输运型PDF三种比较有代表性的燃烧模型对Sandia实验室Flame D进行了模拟，分析不同模型的特性。

1 湍流和燃烧模型 1.1 LES湍流模型

在LES模型中，使用滤波函数将每个变量分为两部分，大尺度的滤波后的变量，在LES模拟时直接计算；小尺度分量，需要用模型来计算表示。经过盒式滤波函数滤波后的N-S方程为

$ \frac{{\partial \bar \rho }}{{\partial t}} + \frac{{\partial \bar \rho {{\tilde u}_i}}}{{\partial {x_i}}} = \overline {{S_\rho }} $

(1)

$ \begin{array}{*{20}{c}} {\frac{{\partial \bar \rho {{\tilde u}_i}}}{{\partial t}} + \frac{{\partial \left( {\bar \rho {{\tilde u}_i}{{\tilde u}_j}} \right)}}{{\partial {x_j}}} = }\\ { - \frac{{\partial \bar \rho }}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left( {{{\bar \tau }_{ij}} - \bar \rho \widetilde {{{u''}_i}{{u''}_j}}} \right) + {{\bar S}_{{u_i}}}} \end{array} $

(2)

$ \begin{array}{*{20}{c}} {\frac{{\partial \bar \rho {{\tilde Y}_k}}}{{\partial t}} + \frac{{\partial \left( {\bar \rho {{\tilde Y}_k}{{\tilde u}_j}} \right)}}{{\partial {x_j}}} = }\\ {\frac{\partial }{{\partial {x_j}}}\left( { - \overline {{V_{k,j}}{Y_k}} - \bar \rho \widetilde {{{u''}_j}{{Y''}_k}}} \right) + \overline {{{\dot \omega }_k}} + \bar S_{{Y_k}}^e} \end{array} $

(3)

$ \begin{array}{*{20}{c}} {\frac{{\partial \bar \rho \tilde h}}{{\partial t}} + \frac{{\partial \left( {\bar \rho \tilde h{{\tilde u}_j}} \right)}}{{\partial {x_j}}} = }\\ {\frac{{D\bar p}}{{Dt}} - \frac{\partial }{{\partial {x_j}}}\left( {{{\bar q}_j} + \bar \rho \widetilde {{{u''}_j}h''}} \right) + \overline {{\tau _{ij}}\frac{{\partial {u_i}}}{{\partial {x_j}}}} + \bar S + \bar S_H^e} \end{array} $

(4)

式中：ρ为流体密度，u为流体速度，Y为组分；$ {\bar w_k}$为组分生成速率，h为组分热焓，S为源项，p为压力，τ、V和q分别为动量方程、组分方程和能量方程的应力项、扩散项和热传导项。可以看出，LES经过滤波后的N-S方程与RANS雷诺时均后的N-S方程形式完全一致，但是代表的意义却完全不同，区别在于这里的变量是滤波后的值，仍为瞬时值，而非时均值。组分和能量方程中的二阶矩项按下式进行封闭，

$ \bar \rho \widetilde {{{Y''}_k}{{u''}_j}} = - \frac{{{\mu _t}}}{{S{c_t}}}\frac{{\partial {{\tilde Y}_k}}}{{\partial {x_i}}}or\bar p\left( {\widetilde {{Y_k}{u_j}} - {{\tilde Y}_k}{{\tilde u}_j}} \right) = - \frac{{{\mu _{{\rm{sgs}}}}}}{{S{c_{{\rm{sgs}}}}}}\frac{{\partial {{\tilde Y}_k}}}{{\partial {x_i}}} $

(5)

$ \begin{array}{*{20}{c}} {\bar \rho \widetilde {h''{{u''}_j}} = - \frac{{{\mu _t}}}{{\mathit{P}{\mathit{r}_t}}}\frac{{\partial \tilde h}}{{\partial {x_i}}}}\\ {\bar \rho \left( {\widetilde {h{u_j}} - \tilde h{{\tilde u}_j}} \right) = - \frac{{{\mu _{{\rm{sgs}}}}}}{{\mathit{P}{\mathit{r}_{{\rm{sgs}}}}}}\frac{{\partial \tilde h}}{{\partial {x_i}}}} \end{array} $

(6)

和RANS封闭方法类似，SGS应力项也可以用涡粘假设描述，

$ {R_{ij}} = \bar \rho \left( {\widetilde {{u_i}{u_j}} - \widetilde {{u_i}{u_j}}} \right) = - 2{\mu _{{\rm{sgs}}}}\tilde S_{ij}^D + \frac{2}{3}\bar \rho k{\delta _{ij}} $

(7)

为了求解u_sgs，提出了Smagrinsky，动态Smagrinsky和一方程等模型。

1.2 燃烧模型 1.2.1 部分搅拌反应模型

PaSR模型最早由Golovitche等^[19]在WSR模型的基础上提出的，其认为燃烧场中某一点的状态变化可以由两部分组成：一部分是基于良搅拌的化学反应变化；另一部分是由于湍流的混合影响造成的变化，即

$ \dot \omega _r^{{\rm{WSR}}}{\mathit{\boldsymbol{\tau }}_{{\rm{chem}}}} = {{\dot \omega }_r}\left( {{\mathit{\boldsymbol{\tau }}_{{\rm{chem}}}} + {\mathit{\boldsymbol{\tau }}_{{\rm{mix}}}}} \right) $

(8)

式中：τ_chem为化学反应时间尺度，τ_mix为混合时间尺度，$\dot \omega _r^{{\rm{WSR}}} $为WSR模型反应速率，${{\dot \omega }_r} $为化学反应速率。

1.2.2 火焰面/反应进度变量模型

当湍流流场中的Kolmogorov涡旋尺度大于火焰面的厚度时，湍流脉动对整个火焰面内部结构无影响，即可以将火焰面内部结构与湍流对火焰面的作用分开处理，从而将湍流流动与化学反应解耦。通过一维火焰计算湍流火焰中的化学热力学状态值${\tilde \varphi } $，并建成以Z，Z″，C为参变量的三维数据表；求解LES湍流输运方程的连续和动量方程，求解Z、Z″、C的输运方程：

$ \frac{{\partial \bar \rho \tilde Z}}{{\partial t}} + \frac{{\partial \bar \rho {{\tilde u}_j}\tilde Z}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left( {\bar \rho \tilde D\frac{{\partial \tilde Z}}{{\partial {x_j}}} - {{\rm{\Phi }}_{Z,j}}} \right) $

(9)

$ \begin{array}{*{20}{c}} {\frac{{\partial \bar \rho {{\tilde Z''}^2}}}{{\partial t}} + \frac{{\partial \bar \rho {{\tilde u}_j}{{\tilde Z''}^2}}}{{\partial {x_j}}} = }\\ {\frac{\partial }{{\partial {x_j}}}\left( {\bar \rho \tilde D\frac{{\partial {{\tilde Z''}^2}}}{{\partial {x_j}}} - {{\rm{\Phi }}_{{Z_v},j}}} \right) - 2\bar \rho {{\tilde s}_{{\chi ^z}}} - 2{\mathit{\Phi }_{Z,j}}\frac{{\partial \tilde Z}}{{\partial {x_j}}}} \end{array} $

(10)

$ \frac{{\partial \bar \rho \tilde C}}{{\partial t}} + \frac{{\partial \bar \rho {{\tilde u}_j}\tilde C}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left( {\bar \rho \tilde D\frac{{\partial \tilde C}}{{\partial {x_j}}} - {\mathit{\Phi }_{C,j}}} \right) + {{\bar{\dot{\omega }}}_{C}} $

(11)

式中：Φ_{Z, j}、Φ_{C, j}为混合分数和反应进度变量的湍流通量，根据标准梯度扩散假设来进行模型模型封闭。

1.2.3 联合标量概率密度函数

定义标量ψ=[Y₁, Y₂, …, Y_Ns, h]的单点单时间联合概率密度函数$\tilde P\left( {\psi ;x, t} \right) $，经过计算，可以得到${\tilde P} $的输运方程，

$ \begin{array}{*{20}{c}} {\bar \rho \frac{{\partial \tilde P\left( \psi \right)}}{{\partial t}} + \bar \rho {{\tilde u}_j}\frac{{\partial \tilde P\left( \psi \right)}}{{\partial {x_j}}} + \frac{\partial }{{\partial {\psi _\alpha }}}\left[ {\bar \rho {{\dot \omega }_\alpha }\left( \psi \right)\tilde P\left( \psi \right)} \right] = }\\ { - \frac{\partial }{{\partial {x_i}}}\left[ {\left\langle {{{u''}_i}\left| \psi \right.} \right\rangle \bar \rho \tilde P\left( \psi \right)} \right] + \frac{\partial }{{\partial {\psi _\alpha }}}\left[ {\left\langle {\frac{{\partial J_i^\alpha }}{{\partial {x_i}}}\left| \psi \right.} \right\rangle \tilde P\left( \psi \right)} \right]} \end{array} $

(12)

式中：〈·|ψ〉表示条件平均。方程左侧是封闭的，第一和第二项表示$\tilde P\left( {\psi ;x, t} \right) $沿平均流线的输运，很显然化学反应源项${\dot \omega } $在TPDF输运方程中以封闭的形式出现，不需要任何近似，这也是TPDF模型的优势所在。

本文使用非线性小尺度混合模型，表示为

$ \left\langle {\frac{{\partial J_i^\alpha }}{{\partial {x_i}}}\left| \psi \right.} \right\rangle = \frac{1}{2}{C_\varphi }{\omega _T}\left( {{\psi _\alpha } - {{\tilde \psi }_\alpha }} \right) $

(13)

式中：ω_T为湍流频率；C_φ是模型系数，对于不同火焰可能不同，一般取2.0。

由于PDF输运方程的高维度，本文采用欧拉随机场的Monte Carlo方法求解PDF方程。认为燃烧场的N_s个标量；ξ_αⁿ(x, t)，1≤n≤N_F, 1≤α≤N_s，由N_F个随机场来表示。颗粒系综的标量演化则遵循随机微分方程(stochastic differential equation，SDE)，

$ \begin{array}{l} \bar \rho {\rm{d}}\xi _\alpha ^n = - \bar \rho {{\tilde u}_i}\frac{{\partial \xi _\alpha ^n}}{{\partial {x_i}}}{\rm{d}}t + \bar \rho S_\alpha ^r\left( {{\psi ^n}} \right){\rm{d}}t + \frac{\partial }{{\partial {x_i}}}\left( {{\mathit{\Gamma }_t}\frac{{\partial \xi _\alpha ^n}}{{\partial {x_i}}}} \right)dt - \\ \;\;\;\;\;\;\;\;\;\;\frac{1}{2}\bar \rho {C_\varphi }{\omega _T}\left( {\xi _\alpha ^n - {{\tilde \varphi }_\alpha }} \right){\rm{d}}t + \bar \rho \sqrt {2\frac{{{\mathit{\Gamma }_t}}}{{\bar \rho }}} \frac{{\partial \xi _\alpha ^n}}{{\partial {x_i}}}{\rm{d}}W_i^n \end{array} $

(14)

式中：$ \bar \rho S_\alpha ^r({\psi ^{(n)}}){\rm{d}}t$为化学反应源项增量，dWⁿ是一个向量Wiener过程。

1.3 燃烧室几何模型

Flame D为美国sandia实验室测量的带有值班火焰的甲烷/空气射流火焰。该系列火焰包含了火焰稳定，局部熄火再燃，污染物NO_x排放等实验数据，为模型验证提供了很好的参考依据。Flame D的燃料进口直径D=7.2 mm，平均速度为49.6 m/s，值班火焰射流进口直径为18.2 mm，温度1 880 K左右。计算域为15D×15D×75D，最小网格尺度为0.1 mm，网格数量为2.3 M，如图 1所示。考虑到PDF和PaSR的计算量，本文分别使用GRI 3.0详细反应机理(53组分，325步反应)^[20]和SKEL简化机理(16组分，41步反应)^[14]。

2 湍流燃烧室模型对比分析

图 2为不同模型计算的温度和速度沿中轴线的分布。FPV模型温度上升较快，但整体符合实验趋势。PaSR模拟的温度最高，因为系数是根据湍流雷诺数估算的，虽然是可以模拟多数火焰的经验值，但可能对Flame D来说不太准确，还需要后续调整。TPDF的温度场是吻合最好的。因为射流出口的剪切层影响燃料和氧化剂的掺混，因此速度入口边界条件对模拟结果的影响非常大。速度分布显示，采用出口回收方式形成充分发展的湍流入口的方法是比较准确的。

图 3为温度和速度场的中截面分布。从温度场可以看出，射流出口形成的剪切层脉动很强烈，强化燃料和空气的掺混，有利于燃烧发生。速度场中很明显的形成了充分发展的湍流，在壁面有很薄的边界层，与实验和DNS模拟速度分布一致。火焰指数FI(takeno flame index)，可以用来表示火焰趋向预混还是扩散的程度，FI= $\nabla {Y_{\rm{F}}} \cdot \nabla {Y_{\rm{O}}} $，扩散火焰部分如图 3(a)的封闭曲线区域。Flame D为典型的扩散主导的火焰，火焰高温区分布在恰当化学当量比的混合分数(Z_st)周围，图 3(b)的实线为Z_st，甲烷和氧气被快速的消耗掉，形成燃烧产物并放出大量的燃烧热。

对于中轴线上燃料消耗和燃烧产物的分布，如图 4所示，相比于实验值，FPV在开始阶段反应较慢，因为忽略了湍流强化燃烧的过程，随着反应过程的发展，CH₄消耗变快，刚好对应于温度快速上升阶段，弱化了湍流对温度场带来的扩散和热传导作用。燃烧产物与燃料变化的过程刚好相反。

对于中间产物的模拟，如图 5所示，FPV模型的模拟都相对弱一些，对于OH和CO都有明显的高估。另外，在创建火焰面库的时候认为刘易斯数为动态变化可能会稍微改善一下模拟准确性。对比3种湍流燃烧模型的准确性，TPDF的表现最好，与实验值对比最好，PaSR和FPV次之。

图 6是中间产物的中截面分布。在燃烧火焰的实验中，认为CH₂O代表火焰的预热区，HCO代表火焰的强烈反应区内层，OH代表燃料发生氧化区。在预热区，化学反应速度较小，一般认为该区域没有化学反应发生，反应区释放的热量通过扩散或热传导使反应物逐渐预热；在进入到反应区之前，燃料开始分解成为CH₂O，随后OH等活化基的形成使CH₂O迅速消耗，并放出热量(HRR)，反应式为CH₂O+OH→HCO+H₂O，因此CH₂O和OH区域中间便是反应区，可以用HRR或HCO分布区域表示。在化学反应区，化学反应较为剧烈，燃料被转换成其它组分，反应速度、温度和活化中心的浓度都达到了最大值；在氧化区，燃烧生成的中间产物最终被转化为生成物。

经过对扩展的泽尔多维奇NO_x机理和敏感性分析来看，以下是NO直接生成率较高的基元反应，

$ {\rm{N}} + {\rm{OH}} \to {\rm{NO}} + {\rm{H}} $

(15)

$ {\rm{HNO}} + {\rm{OH}} \to {\rm{NO}} + {{\rm{H}}_2}{\rm{O}} $

(16)

说明中间组分，如OH模拟的准确性，将影响到NO的准确程度。因为使用NO_x机理计算量比较大，因此本文借助Wang的计算结果(RANS+PDF+GRI2.11)来对比FPV的NO_x结果^[21]，如图 7所示。结果很明显，FPV模型模拟出来的NO要比实验值高出几倍，虽然趋势是和实验值相符。而TPDF模型则展现出了对NO模拟的优势，结果与实验值非常接近。

对于FPV模型，由于NO形成具有很强的温度敏感性，因此需要在考虑辐射热损失的建库方法才能准确地预测火焰温度，得到更准确的NO质量分数。与其他组分不同，NO的形成是一个缓慢的过程，受动力学反应速率限制，NO的平均值不能准确的从火焰面库插值得到。因此加入一个新的NO输运方程来改善NO_x的模拟结果^[22]，表达式为

$ \begin{array}{*{20}{c}} {\frac{{\partial \bar \rho \widetilde {{Y_{{\rm{NO}}}}}}}{{\partial t}} + \frac{{\partial \bar \rho {{\tilde u}_j}\widetilde {{Y_{{\rm{NO}}}}}}}{{\partial {x_j}}} = }\\ {\frac{\partial }{{\partial {x_j}}}\left( {\bar \rho \tilde D\frac{{\partial \widetilde {{Y_{{\rm{NO}}}}}}}{{\partial {x_j}}} - {\mathit{\Phi }_{{\rm{NO}},j}}} \right) + \bar \rho \widetilde {{{\dot \omega }_{{\rm{NO}}}}}} \end{array} $

(17)

$ \widetilde {{{\dot \omega }_{{\rm{NO}}}}} = \int\limits_0^1 {\int\limits_0^\infty {{{\dot \omega }_{{\rm{NO}}}}\left( {Z,{\chi _{st}}} \right)\tilde P\left( {Z,{\chi _{st}}} \right){\rm{d}}Z{\rm{d}}{\chi _{st}}} } $

(18)

图 8为x=15D和30D截面上的温度和OH组分以混合分数为横坐标的散点图。可以看出，温度最高值刚好对应于Z_st，而且30D截面的最高温度比15D略高一些，因为火焰在x=15D以后会存在再燃现象，在Flame E和F上会更明显。数值模拟的散点分布在平均值曲线周围，吻合良好，但是与实验测量的散点分布有略微区别，主要是PDF输运方程小尺度混合项的模型准确度带来的误差。

各个模型所使用的CPU时间如表 1所示。理论上分析，FPV(建库需要1 d时间左右)只求解很少的输运方程，而PaSR需要求解多个组分输运方程，可能是因为化学反应刚性较小，所以求解化学反应占用的时间并不是很长。如果遇到更大的且刚性比较强的反应机理，可以用ISAT或者CCM来加速化学反应。本文使用的TPDF模型用8个随机场模拟求解PDF输运方程，所以耗时差不多是PaSR的8倍左右。尽管TPDF模型的模拟结果最好，但是对于计算能力要求较高，而FPV的计算速度是工程应用上最大的优势所在。

3 结论

1) 三种模型均能较好的模拟扩散火焰的温度速度和组分分布趋势，TPDF的计算精度最高，PaSR次之，FPV对于中间组分的模拟与实验值差别较大。

2) 对NO_x机理分析表明，中间组分OH和温度对NO_x结果影响最大，因此对于OH组分模拟较差的FPV模型需要加入NO输运方程来改善NO_x的模拟结果。

3) 相比TPDF和PaSR模型，FPV的方法在计算量上大大缩减，在工程应用中具有很大价值，但需要进一步的研究修正。