重型燃气轮机燃烧室展现出更高参数、更低排放的发展趋势：H/J级重型燃机的工作压力已达20~27个大气压(1个大气压= 1.013 MPa)，燃烧室出口温度高达1 700 ℃，NO_x排放已达到不超过10 mg/m³水平^[1-2]。在燃烧室研发过程中，通过计算流体力学(computational fluid dynamics, CFD)准确预测NO_x十分重要。由于计算速度快、计算资源消耗少等优点，Reynolds时均N-S方程(Reynolds-averaged Navier-Stokes equations, RANS) 被广泛应用于燃烧室模拟中。目前在商用软件平台上，求解NO_x生成多采用后验方法，即首先求解出流场、温度场、组分场，再单独解NO_x的生成反应方程^[3-4]。实践表明，这种后验方法对于燃烧室的复杂结构、不同负荷和工作模式缺乏足够的普适性，计算结果还会随不同的湍流燃烧模型而不同，给准确预测NO_x排放带来了很大的不确定性。

为解决上述问题，本文采用化学反应器网络(chemical reactor network, CRN)方法取代后验方法，进而发展成一个基于CFD与CRN耦合的完整、准确描述湍流、燃烧和污染物生成的数值模拟方法。传统的CRN方法主要用于对NO_x排放的快速准确预测，其原理是将整个燃烧室空间划分成预混区、火焰锋面区、主火焰区、掺混气流区等有限个小区域，并处理成不同类型的反应器，再根据流量分配将其连成一个网络。本文提出的方法将CRN作为数值模拟的一环，根据CFD结果统计出反应器体积、初始温度、燃料和空气流量分配等关键参数后，再根据详细的化学反应机理计算NO_x^[5-8]。这与在百千万个网格上迭代求解守恒方程组的CFD做法不同，具有快速、准确的优点。

对于CFD-CRN耦合方法，对不同结构的燃烧室进行高效合理分区以及确定分区数量，是确保NO_x预测精度的关键。为此，需要研究CRN燃烧空间自动分区及求解算法、分区策略和通用准则。

现有CRN自动分区/求解算法主要有两类：

1) 将流场中温度或组分相近的网格合并成均匀搅拌反应器(perfectly-stirred reactor, PSR)，导致反应器数量较多。然后，用反应器平均温度和详细反应机理求解化学反应过程^[9-11]。

例如，有研究用温度和化学当量比进行初步分区，再用混合不均匀指数对初步分区进行合并和重组，以减少反应器数量，再用软件DSMOKE求解锅炉炉膛的污染物排放，但反应器数量仍达到400个^[9]；还有以温度和轴向速度为分区标准，用C语言将自动构建算法编写成FLUENT用户定义函数(user define function, UDF)，再用反应动力学后处理软件KPPSMOKE求解，获得旋流非预混火焰的污染物排放值，但PSR反应器数量更高达5 000个^[11]；还有用Python语言编写算法，根据预设的相似性标准，自动聚类CFD计算所用的体网格成PSR反应器，用开源软件Cantera预测实验室尺度的燃烧器污染物生成水平，反应器数量也在1 000个以上^[10]。对反应器数量庞大的CRN进行整体求解，需要占用计算机大量的内存，故只能使用单反应器分步求解迭代算法。

2) 按一定的分区准则来划分流场，可大幅减少反应器数量，再根据详细反应机理求解化学反应过程和能量方程，获得反应器的组分和温度^{[8, 12]}。

例如，有研究根据预先设定的O₂/CH₄摩尔分数、轴向速度和混合分数的阈值进行分区，再用非线性温度间隔对高温区分区进行优化，然后用商业软件ENERGICO自动构建贫预混燃烧器的CRN，再导入商业软件CHEMKIN算出NO_x排放值，反应器数量只有84个^[12]；还有研究用Python语言开发了燃烧室CRN自动生成程序，根据燃烧室的几何结构、流场、温度场特征，预先定义分区标准，再基于网格聚类算法生成结构清晰的CRN模板，最终用软件Cantera整体耦合求解出所有反应器的平均组分和温度。由此形成的反应器数量在几十个量级，边界连接清晰，如果用于预测相近工况的NO_x排放，只需调整边界条件即可^[8]。

上述CRN自动分区/求解算法存在的问题是缺乏通用性，即对于不同的CFD工况，需要重新定义或调整分区准则和CRN的初始参数，并修改程序^[8]。因此，很有必要开发通用算法和程序。

在构建燃气轮机燃烧室化学反应器网络时，现有研究多针对个案应用，未从方法层面提出合理的分区策略和通用准则，因此无法确保不同燃烧室结构或工况下的NO_x预测精度。现有研究所采用的分区准则参数很多，但对于究竟哪些参数最为合适并无明确和统一的结论^[13-16]。例如，有研究用局部当量比和温度构建的5 MW和10 MW天然气玻璃窑炉化学反应器网络，对满负荷的NO_x预测误差约为3%^[17]。通过进一步优化，改用混合分数和过程变量后，预测误差降至2.7%，但在75%负荷时误差又升至63.5%，其反应器数量已多达369个^[18-19]。可见，反应器数量并非越多越好，否则会增加计算量，削弱CRN方法快速预测的优越性。对于燃气轮机燃烧室，针对主火焰区有多种分区方法。例如，有研究用单个变量——混合分数Z，以间隔值0.01将火焰锋面分成上、下游两部分，对扩散燃烧的NO_x预测误差小于6%^[20]；还有用CO质量分数最大值的20%等值线、1 800 K等温线、OH和CH₂O自由基的分布等^{[3, 6, 14, 21]}进行分区。总之，这些方法仅适用于确定结构燃烧室的单一工况。

在计算效率方面，CRN方法简化了对燃烧室湍流流动的模拟，将燃烧室划分成数个区域并使用反应器代替，而反应器的数量远小于CFD的计算网格数，因此计算时间较短。CRN方法计算时间通常比CFD低数个数量级。文[18]对CRN的计算时间的统计结果表明，对于一个包含369个反应器的CRN，不同工况下的计算时间为887~992 s，在并行计算条件下，计算时间可以降低至1 min内。

为了提高CFD-CRN耦合方法燃烧室数值模拟的NO_x预测精度，本文以某国产重型燃气轮机全尺寸全参数天然气干式低污染燃烧室(dry low NO_x, DLN)为对象，通过实验和数值计算，研究CRN自动分区/求解方法及其程序化，以及适用于多工况、不同燃烧模式、不同空气流量分配的分区策略和通用准则。

1 数值模拟的数据准备

按照燃烧室数值模拟的基本流程，需先完成流场、温度场、燃烧组分场的迭代求解，获得用于构建CRN的基础数据。

本文所采用的燃烧室由导流衬套、头部喷嘴组件、火焰管和过渡段组成，如图 1所示。头部喷嘴组件按燃料径向分级模式，分为中心喷嘴和环形区喷嘴两类。空气沿导流衬套和火焰管之间的环形通道流向头部，与燃料天然气(以甲烷CH₄代表)混合后进入燃烧区燃烧。部分空气经火焰管下游的掺混孔和头部环形区外侧的主燃孔进入燃烧区。

中心区始终按照贫预混燃烧方式工作，环形区则有两种工作模式：当相对负荷K_L＜0.8时始终按部分预混燃烧方式工作，在0.8以上直至满负荷均保持熄火状态，环形区作为预混通道使用。

实验和数值模拟的工况参数和NO_x排放如表 1所示，涉及不同负荷、不同掺混孔径和环形区的两种工作模式。其中，掺混孔径D_jet用于调节头部燃烧区当量比及火焰温度，以控制NO_x排放。模拟中用到的燃料天然气理论燃烧空气量L₀ (质量比) 取值17.11。其他有关计算网格、湍流燃烧模型、边界条件的信息详见文[22]。

图 2为数值模拟求出的燃烧室温度场。此时环形区按熄火模式工作，故区内的气流温度低，喷入的燃料与空气混合后进入主燃烧区燃烧，与中心喷嘴形成双层高温火焰锋面。掺混射流也因掺混孔径大而十分强劲，几乎贯穿整个火焰筒，下游的整体燃气温度迅速变得相当均匀。当然，此时的NO_x排放明显比小掺混孔径时高出数倍(见表 1)。

2 CRN方法的改进

对CRN方法进行研究的目的是提高其计算效率和对于不同燃烧室结构及工况的通用性。为此，需要研究改进CRN方法的两个关键环节：1) 自动分区/求解的算法和程序化，2) 分区策略和通用准则。

2.1 CRN自动分区/求解算法和程序化

这部分主要包括数据格式标准化、程序开发和数值模拟结果的后处理方法。

1) 数据格式标准化。使用XML格式作为统一的数据存储标准，实现程序输入输出信息的规范化。

2) 程序开发。通过Python语言完成CRN自动分区/求解算法开发：分区程序可根据所输入的分区准则完成区域划分、网格聚类和参数统计，并生成CRN网络；求解程序基于开源反应动力学模拟软件Cantera，根据各反应器参数及质量流量传递求解Gri-Mech 3.0详细反应动力学机理，计算得到所有反应器的详细组分和温度，并换算成15%O₂体积分数下的NO_x体积分数。

3) 数值模拟结果的后处理方法。在商用软件Ansys FLUENT算出的燃烧场基础上，将需要用到的体网格标识(identification, ID)及体网格的体积和坐标，以及相应的温度、压力、组分等导出为CGNS文件、HDF5格式的算例和数据文件。

程序的操作步骤为：

1) 准备自动分区模块的输入文件.xml。该文件中包含CFD数据文件的存储路径、数值模拟用到的物理量名称、流体域名称及其ID、内部面和边界面名称及其ID、分区准则关系式、NO_x生成详细反应机理文件名(例如gri30.yaml)、燃料和氧化剂信息等。

2) 运行自动分区模块。选择上述输入文件.xml，从CFD数据文件中筛选出所需的物理量数据，采用网格聚类算法将符合某一分区准则的体网格聚并到一个PSR反应器中，并顺序完成对应所有分区准则的网格聚并操作；计算反应器参数(体积V、压力p、温度T、组分等)和反应器间的质量流量。

反应器体积等于所包含的体网格i(i=1, 2, …, N)的体积之和，即V=∑V_i；反应器压力p和温度T则用体积加权平均计算得到：

$ p=\frac{\sum V_i p_i}{\sum V_i}, T=\frac{\sum V_i T_i}{\sum V_i} . $

(1)

完成上述分区及参数计算之后，将有关信息存储为CRN求解模块的输入文件.xml，此后视需要还可随时改动其中的一些参数值。

3) 运行CRN求解模块。读入上一步产生的输入文件，通过接口程序与开源反应动力学模拟软件Cantera连接并进行求解，即可获得各反应器的温度和详细组分信息以及燃烧室出口的NO_x体积分数。为保证各反应器进出口的质量流量平衡，需计算反应器所有出口的流量分配比例并保持不变，再用实际输入的质量流量对各反应器间的质量流量进行校正。

以上算法程序实现了基于数值模拟结果的CRN预测NO_x排放的整个自动求解过程。相比于求解NO_x的后验方法，其计算耗时几乎可以忽略。除了手动完成Ansys FLUENT数据输出和输入文件的编辑等预备性操作外，对于反应器数量约为30个的CRN，上述算法程序在笔记本电脑上只需3 min即可完成自动分区并输出结果，计算效率非常高。

完成上述求解过程的关键在于分区准则，而分区准则是否合理和通用则取决于正确的分区策略，由此决定了反应器数量及计算效率、NO_x预测准确性和方法的普适性。

2.2 分区策略和通用准则

正确的分区策略应体现出燃烧室的主要特征，并符合关于NO_x生成与排放的燃烧基础理论。

现代DLN燃烧室有两个主要特征：1) 采用燃料分级的多喷嘴结构，燃料分级包括轴向分级、径向分级、多点式分级等多种形式，必然形成多个火炬或多重火焰结构，导致空间温度分布不均匀，如图 2所示。显然，在进行CRN分区时，必须兼顾由燃烧引起的温度轴向(沿流动方向)变化以及多重火焰结构导致的温度径向变化。

2) 采用贫预混燃烧模式。由燃烧理论可知，贫预混燃烧模式可通过调节当量比将火焰平均温度控制在一定范围内，进而实现对NO_x总体排放水平的控制。但是，燃料与空气预混的均匀程度直接关系到燃烧空间中可能存在呈扩散燃烧性质的炽热高温区，进而对NO_x排放产生关键影响。由图 3可见，尽管燃烧室设计有一定长度的预混距离，在中心喷嘴出口端面(图 2中的截面A—A)和环形区出口(文丘里喉部截面B—B)处的甲烷质量分数(w_CH4)分布仍很不均匀。故在进行CRN分区时，尤应考虑燃烧室横截面内燃料质量分数和温度分布的不均匀性。

以上认识主要针对有火焰的燃烧反应区。实际上任何一个燃烧室均可按气流流动方向顺序分成纯流动区、燃烧反应区、后燃烧区3类大区，在它们之下再细分出用于CRN计算的小区，如表 2所示。

纯流动区只涉及空气、燃料(甲烷CH₄)以及二者的混合物，并未发生燃烧，因此使用气流温度和组分的质量分数w作为分区准则较为合适，即w_O2≥ 0.229为纯空气，w_CH4≥0.999为纯甲烷燃料。T≤T₂为空气/燃料的预混区，其中T₂是燃烧室入口空气温度(由于入口燃料温度低于空气温度，在发生燃烧前，预混区温度不超过T₂)。同时，考虑预混不均匀性, 使用预混区的平均当量比对预混区进行划分。本文所用燃烧室的预混区有两个：1) 中心喷嘴的预混通道，2) K_L≥0.8的环形区(参见表 1)。

燃烧反应区无疑是大量生成NO_x的关键区域。研究表明^{[3, 6]}，该区的分区准则对CRN计算精度有重大影响。其中，回流区的划分很容易，以“轴向速度v_a=0”作为分区准则，其等值面的包络区就是回流区。除了中心喷嘴下游的回流区外，在环形区出口的文丘里结构背面也有一个小回流区。

对火焰锋面区和主火焰区的分区准则相对复杂。本文研究发现，采用“燃尽率η”和“燃烧前的区平均当量比φ_av”最为合理。从物理上讲，燃尽率体现了燃烧释热程度和燃烧效率，决定了总体控制NO_x生成的平均温度水平。由于本文所用燃烧室采用了径向燃料分级，即沿径向分成了中心区和外侧的环形区，且不同区域对NO_x的贡献不同，故要分别计算两区的燃尽率，

$ \eta_k=\left(1-\frac{w_{\mathrm{CH}_4}}{w_{\mathrm{CH}_{4, \text { in }}}}\right)_k \times 100 \%, k=\mathrm{c}, \mathrm{r} . $

(2)

式中：下标c代表中心区，r代表环形区。式(2)中的w_{CH4, in}为该区进口的甲烷平均质量分数，w_CH4为燃烧后的剩余值。相较于温度准则，使用燃尽率η能更准确地表达燃烧反应进度。

燃烧前的区平均当量比φ_av，实际上代表的是燃烧前的燃料质量分数。如果已知任一CRN轴向分区后得到的预混区、火焰锋面区和出火焰区的平均当量比，就能分辨出上述各区内部的燃料质量分数的高低差异。通过比较局部当量比φ与各区平均当量比的大小再细分子区，φ≥φ_av为燃料浓区，φ＜φ_av为燃料稀区，就能计算燃料/空气预混的不均匀度，进而得到预混不均匀度对NO_x的关键影响。φ可用数值模拟算出的燃烧产物质量分数来计算，

$ \varphi=L_0 \frac{0.233}{0.75} \frac{\frac{3}{4} w_{\mathrm{CH}_4}+\frac{3}{11} w_{\mathrm{CO}_2}+\frac{3}{7} w_{\mathrm{CO}}}{w_{\mathrm{O}_2}+\frac{8}{11} w_{\mathrm{CO}_2}+\frac{4}{7} w_{\mathrm{CO}}+\frac{8}{9} w_{\mathrm{H}_2 \mathrm{O}}} . $

(3)

综合以上研究，本文提出了对燃烧反应区的分区策略和通用准则：

1) 采用“两维度”分区法。先进行轴向/径向(纵剖面内)分区，然后完成横截面内分区。

2) 轴向/径向分区准则。可根据燃尽率η₁=99%和η₂=90%，确定火焰锋面的轴向/径向边界。

3) 横截面内分区准则。根据燃烧前的区平均当量比φ_av，可将已有的轴向/径向分区进一步分成稀、浓两区。

4) 环形区处于燃烧状态时(K_L≤0.8)，只在区内按以上策略和准则分区；当处于熄灭状态时(K_L≥0.8)，则在其下游燃烧区的外围进行分区。

上述分区策略和通用准则有如下优点：

1) 体现了燃烧空间各个维度上的不均匀性特征，充分反映了燃烧室结构特征和燃料/空气预混不均匀度对NO_x生成的关键影响；

2) 在不同维度上采用不同的分区准则，以及先粗后细的分区顺序，体现了不同物理过程的特征；

3) 可区别对待环形区有无燃烧的不同工作模式，对预测部分负荷燃烧室的NO_x排放有重要意义。

下面将通过数据验算证明本文提出的分区策略和通用准则(简称“η-φ准则”)的合理性和准确性。

3 应用结果与分析

本章采用第2章发展的CRN方法对表 1的燃烧工况进行计算，作为对比，还给出了文献中常用的温度准则分区方法的预测结果。

3.1 基于温度准则的NO_x预测结果

采用文献中常用的温度准则分区方法^[8]，先按表 2分出纯流动区、燃烧反应区和后燃烧区3个大区，再用不同的温度间隔将燃烧反应区，特别是NO_x主要生成区域，进一步分出若干小区，结果如图 4所示。然后，用基于温度准则的CRN(CRN-T)预测对应表 1各工况的燃烧室出口NO_x体积分数，结果如表 3所示。

可见，基于温度准则的CRN方法完全不能适应工况的变化，除工况5外，其余工况的预测误差均超过70%，通用性很差。

参照文[18]中提高预测精度的惯常做法，本文进一步以工况2为例，研究了缩小温度间隔、增加分区即反应器数量的效果。由图 5可见，当温度间隔从初始的ΔT=200 K依次缩小到100、50、25 K，相应地，反应器数量从33增加到43、61、98个后，NO_x预测值随之增大约5倍，但仍与实验值有很大差距，且增大趋势逐渐趋缓，明显形成一个系统误差。同时，计算耗时增加很快。这表明，单纯增加反应器数量的做法缺乏实用价值，除了如加密数值计算网格的作用一样，从量上对NO_x预测分辨率有所改进外，对预测精度和计算耗时均有负面作用。从根本上讲，预测精度低是由于分区策略失当、分区准则不合理、未体现燃料/空气预混不均匀性的关键影响所致。

3.2 基于η-φ准则的NO_x预测结果

采用η-φ准则对燃烧反应区进行CRN分区的具体步骤如下：

1) 径向分区。根据燃烧室径向特征，将燃烧区沿径向用直线简单分为3个区，即值班火焰区(7区)、中心火焰区(8、9区)、环形火焰区(10、11区)。

2) 轴向分区。用η₁=99%和η₂=90%依次确定火焰锋面区的两个轴向边界，见图 6a。

3) 横截面内分区。分别用中心区和环形区的入口平均当量比φ_c和φ_r, 以及总平均当量比φ₀，将以上各区再分出浓、稀两个子区，如图 6b和6c所示。

轴向分区个数究竟多少合适，目前文献中无统一结论^[5]。通常认为，火焰锋面区内的温度梯度很大，合理的分区可使区内的参数分布适当均匀化，且有利于算准高温燃气的停留时间和NO_x体积分数。本文分别将火焰锋面区的轴向分区由1增加至4进行计算：首先以η=99%作为该区域的轴向边界，分区个数为1表示火焰锋面区在轴向不分区；分区个数为2时，以η=90%进行再次分区；分区个数为3，则是以η=90%、η=80%进行分区，依此类推；并且在轴向分区后，仍使用φ准则将各区分为燃料浓、稀两类子区。由此得到NO_x体积分数随轴向分区个数的变化，如图 7a所示。

由图 7a可见，在采用η-φ准则下，轴向分区只要超过2个，NO_x体积分数预测值就与实验值几乎完全吻合。如果不分区(N_zone=1)，预测值高达299.5×10^-6，是实验值76×10^-6的4倍。这是因为火焰温度变化剧烈，而CRN却采用的是完全均匀反应器(PSR)，二者反应状态差异甚大，使高温燃气停留时间的计算误差很大。由3.1节已知，分区个数并非越多越好。因此，本文用两个η值分出3个区就显得十分合理，既保证了预测精度，又留有足够的“安全裕度”。

由图 7b可见，在用η₁=99%分出整个火焰锋面区的前提下，用于再分区的η₂值就会使NO_x预测值呈不规律的变化。这意味着，对火焰锋面区的再分区不能随意进行，而要使分出的两区体积(或区内的高温燃气停留时间)形成最佳匹配关系。η₂=90% 正是最佳点，且所有其他工况亦如此，因此η₂的取值也是合理的。

使用φ准则在横截面对燃烧空间进行分区，与轴向分区的情形类似，横截面内的分区个数和分区准则同样直接影响NO_x的预测精度，因为分区必须体现出燃料/空气预混的不均匀特征(见图 6b和6c)。由图 8a可见，根据本文提出的η-φ准则，以中心区No.8区为研究对象，将其细分为多个子区，对分区个数进行检验，具体的划分方法为：分区个数为2时，以整区的平均当量比进行划分；分区个数为3时，则在2个分区的基础上，以浓区的平均当量比进一步划分，依此类推；将序号为No.8的小区再细分成多个子区后，本研究也得到与轴向分区类似的结论，即分区个数超过2的NO_x预测值均与实验值吻合。可见，横截面内必须分区，且只需2~3个分区就完全满足精度要求。

本文进一步针对两个分区的情况，研究了φ取不同的分界值对NO_x预测值的影响。由图 8b可见，当分界值取到待分区(即中心区和环形区)的平均当量比，再细分出稀(φ<φ_av)、浓(φ≥φ_av)两类子区时，预测值最接近实验值。除此之外的其他任何φ值都会产生更大误差。偏离φ_av越远，误差越大，且误差分布相对于φ_av还呈一定的对称性。例如，φ=0.8和0.7均偏离φ_av=0.75约6.7%，会使NO_x的误差达到~10%。

总之，横截面内合理分区和当量比φ的合理取值是确保NO_x预测精度的“充分必要条件”，其中，前者是必要条件。采用φ_av作为分区准则时，横截面内分2个区足以保证预测精度，无需更多，如图 6b和6c所示。

η-φ准则同样适用于环形区有燃烧的情形(工况3)，其分区方案如图 9所示。

最终，基于η-φ准则的CRN方法预测结果如表 3所示。可见，预测精度有实质性的改善，最大误差仅为6.4%。由此证明本文提出的分区策略和通用准则是合理可行的。

4 结论

本文针对化学反应器网络法在燃烧室数值模拟中的应用开展研究，开发了CRN分区的自动生成和求解程序，提出了适用于某国产DLN燃烧室的分区策略和通用准则，并通过实验值进行检验，得到的主要结论如下：

1) 开发了CRN自动分区/求解算法和程序，采用XML数据格式对输入输出数据的格式进行标准化处理，提高了方法在不同燃烧室类型和工况条件下的通用性；

2) 提出了适用于燃气轮机燃烧室的CRN分区策略和通用准则，即η-φ准则，体现了燃烧室结构特征，符合基础燃烧理论的要求，物理意义明确，参数取值合理；

3) 基于η-φ分区准则的CRN方法适用于多负荷、不同燃烧工作模式以及燃烧室局部结构的变化等情况，从根本上提升了NO_x排放的预测精度，最大误差不超过6.4%；

4) 发展了基于CFD-CRN耦合的完整准确描述湍流、燃烧和污染物生成的数值模拟方法，可以完全取代现有商用软件对NO_x排放的后验式计算模型，大幅提高了计算效率和预测准确性。