重型燃气轮机燃烧室展现出更高参数、更低排放的发展趋势：H/J级重型燃机的工作压力已达20~27个大气压(1个大气压= 1.013 MPa)，燃烧室出口温度高达1 700 ℃，NOx排放已达到不超过10 mg/m3水平[1-2]。在燃烧室研发过程中，通过计算流体力学(computational fluid dynamics, CFD)准确预测NOx十分重要。由于计算速度快、计算资源消耗少等优点，Reynolds时均N-S方程(Reynolds-averaged Navier-Stokes equations, RANS) 被广泛应用于燃烧室模拟中。目前在商用软件平台上，求解NOx生成多采用后验方法，即首先求解出流场、温度场、组分场，再单独解NOx的生成反应方程[3-4]。实践表明，这种后验方法对于燃烧室的复杂结构、不同负荷和工作模式缺乏足够的普适性，计算结果还会随不同的湍流燃烧模型而不同，给准确预测NOx排放带来了很大的不确定性。
为解决上述问题，本文采用化学反应器网络(chemical reactor network, CRN)方法取代后验方法，进而发展成一个基于CFD与CRN耦合的完整、准确描述湍流、燃烧和污染物生成的数值模拟方法。传统的CRN方法主要用于对NOx排放的快速准确预测，其原理是将整个燃烧室空间划分成预混区、火焰锋面区、主火焰区、掺混气流区等有限个小区域，并处理成不同类型的反应器，再根据流量分配将其连成一个网络。本文提出的方法将CRN作为数值模拟的一环，根据CFD结果统计出反应器体积、初始温度、燃料和空气流量分配等关键参数后，再根据详细的化学反应机理计算NOx[5-8]。这与在百千万个网格上迭代求解守恒方程组的CFD做法不同，具有快速、准确的优点。
对于CFD-CRN耦合方法，对不同结构的燃烧室进行高效合理分区以及确定分区数量，是确保NOx预测精度的关键。为此，需要研究CRN燃烧空间自动分区及求解算法、分区策略和通用准则。
现有CRN自动分区/求解算法主要有两类：
例如，有研究用温度和化学当量比进行初步分区，再用混合不均匀指数对初步分区进行合并和重组，以减少反应器数量，再用软件DSMOKE求解锅炉炉膛的污染物排放，但反应器数量仍达到400个[9]；还有以温度和轴向速度为分区标准，用C语言将自动构建算法编写成FLUENT用户定义函数(user define function, UDF)，再用反应动力学后处理软件KPPSMOKE求解，获得旋流非预混火焰的污染物排放值，但PSR反应器数量更高达5 000个[11]；还有用Python语言编写算法，根据预设的相似性标准，自动聚类CFD计算所用的体网格成PSR反应器，用开源软件Cantera预测实验室尺度的燃烧器污染物生成水平，反应器数量也在1 000个以上[10]。对反应器数量庞大的CRN进行整体求解，需要占用计算机大量的内存，故只能使用单反应器分步求解迭代算法。
上述CRN自动分区/求解算法存在的问题是缺乏通用性，即对于不同的CFD工况，需要重新定义或调整分区准则和CRN的初始参数，并修改程序[8]。因此，很有必要开发通用算法和程序。
在构建燃气轮机燃烧室化学反应器网络时，现有研究多针对个案应用，未从方法层面提出合理的分区策略和通用准则，因此无法确保不同燃烧室结构或工况下的NOx预测精度。现有研究所采用的分区准则参数很多，但对于究竟哪些参数最为合适并无明确和统一的结论[13-16]。例如，有研究用局部当量比和温度构建的5 MW和10 MW天然气玻璃窑炉化学反应器网络，对满负荷的NOx预测误差约为3%[17]。通过进一步优化，改用混合分数和过程变量后，预测误差降至2.7%，但在75%负荷时误差又升至63.5%，其反应器数量已多达369个[18-19]。可见，反应器数量并非越多越好，否则会增加计算量，削弱CRN方法快速预测的优越性。对于燃气轮机燃烧室，针对主火焰区有多种分区方法。例如，有研究用单个变量——混合分数Z，以间隔值0.01将火焰锋面分成上、下游两部分，对扩散燃烧的NOx预测误差小于6%[20]；还有用CO质量分数最大值的20%等值线、1 800 K等温线、OH和CH2O自由基的分布等[3, 6, 14, 21]进行分区。总之，这些方法仅适用于确定结构燃烧室的单一工况。
在计算效率方面，CRN方法简化了对燃烧室湍流流动的模拟，将燃烧室划分成数个区域并使用反应器代替，而反应器的数量远小于CFD的计算网格数，因此计算时间较短。CRN方法计算时间通常比CFD低数个数量级。文[18]对CRN的计算时间的统计结果表明，对于一个包含369个反应器的CRN，不同工况下的计算时间为887~992 s，在并行计算条件下，计算时间可以降低至1 min内。
为了提高CFD-CRN耦合方法燃烧室数值模拟的NOx预测精度，本文以某国产重型燃气轮机全尺寸全参数天然气干式低污染燃烧室(dry low NOx, DLN)为对象，通过实验和数值计算，研究CRN自动分区/求解方法及其程序化，以及适用于多工况、不同燃烧模式、不同空气流量分配的分区策略和通用准则。
1 数值模拟的数据准备
按照燃烧室数值模拟的基本流程，需先完成流场、温度场、燃烧组分场的迭代求解，获得用于构建CRN的基础数据。
本文所采用的燃烧室由导流衬套、头部喷嘴组件、火焰管和过渡段组成，如图 1所示。头部喷嘴组件按燃料径向分级模式，分为中心喷嘴和环形区喷嘴两类。空气沿导流衬套和火焰管之间的环形通道流向头部，与燃料天然气(以甲烷CH4代表)混合后进入燃烧区燃烧。部分空气经火焰管下游的掺混孔和头部环形区外侧的主燃孔进入燃烧区。
中心区始终按照贫预混燃烧方式工作，环形区则有两种工作模式：当相对负荷KL＜0.8时始终按部分预混燃烧方式工作，在0.8以上直至满负荷均保持熄火状态，环形区作为预混通道使用。
表 1 工况参数和NOx排放实测值及数值模拟预测值 |
| 参数 | 工况编号 | |||||
| 1 | 2 | 3 | 4 | 5 | ||
| 相对负荷KL | 1.0 | 0.5 | 0.8 | 0.8 | 1.0 | |
| 环形区工作模式 | 熄火 | 燃烧 | 燃烧 | 熄火 | 熄火 | |
| 掺混孔径Djet/mm | 50.8 | 69.85 | 69.85 | 69.85 | 69.85 | |
| 空气 | 流量Gair/(kg·s-1) | 15.64 | 15.7 | 15.7 | 15.7 | 15.4 |
| 压力p2/kPa | 1 486 | 1 349 | 1 440 | 1 440 | 1 491 | |
| 温度T2/K | 676 | 656 | 661 | 661 | 675 | |
| 燃料 | 总流量Gf/(kg·s-1) | 0.316 | 0.206 | 0.270 | 0.277 | 0.310 |
| 环形区Gf, r/(kg·s-1) | 0.262 | 0.132 | 0.162 | 0.215 | 0.240 | |
| 中心区Gf, c/(kg·s-1) | 0.053 7 | 0.074 | 0.108 | 0.058 | 0.07 | |
| 温度Tf/K | 288 | 288 | 288 | 288 | 288 | |
| 燃烧室出口平均温度T3 | 实测值/K | 1 483.0 | 1 197.8 | 1 365.9 | 1 365.9 | 1 483.0 |
| CFD计算值/K | 1 463.9 | 1 180.6 | 1 329.8 | 1 342.3 | 1 447.0 | |
| 误差δT3/% | 1.3 | 1.4 | 2.6 | 1.7 | 2.4 | |
| NOx体积分数(@15%O2) | 实测值/10-6 | 29.1 | 87.0 | 65.2 | 33.0 | 76.0 |
| CFD计算值/10-6 | 3.4 | 32.2 | 2.6 | 15.6 | ||
| 误差δCFD/% | 88.3 | 63.0 | 92.1 | 82.7 | ||
2 CRN方法的改进
对CRN方法进行研究的目的是提高其计算效率和对于不同燃烧室结构及工况的通用性。为此，需要研究改进CRN方法的两个关键环节：1) 自动分区/求解的算法和程序化，2) 分区策略和通用准则。
2.1 CRN自动分区/求解算法和程序化
这部分主要包括数据格式标准化、程序开发和数值模拟结果的后处理方法。
1) 数据格式标准化。使用XML格式作为统一的数据存储标准，实现程序输入输出信息的规范化。
2) 程序开发。通过Python语言完成CRN自动分区/求解算法开发：分区程序可根据所输入的分区准则完成区域划分、网格聚类和参数统计，并生成CRN网络；求解程序基于开源反应动力学模拟软件Cantera，根据各反应器参数及质量流量传递求解Gri-Mech 3.0详细反应动力学机理，计算得到所有反应器的详细组分和温度，并换算成15%O2体积分数下的NOx体积分数。
3) 数值模拟结果的后处理方法。在商用软件Ansys FLUENT算出的燃烧场基础上，将需要用到的体网格标识(identification, ID)及体网格的体积和坐标，以及相应的温度、压力、组分等导出为CGNS文件、HDF5格式的算例和数据文件。
程序的操作步骤为：
1) 准备自动分区模块的输入文件.xml。该文件中包含CFD数据文件的存储路径、数值模拟用到的物理量名称、流体域名称及其ID、内部面和边界面名称及其ID、分区准则关系式、NOx生成详细反应机理文件名(例如gri30.yaml)、燃料和氧化剂信息等。
2) 运行自动分区模块。选择上述输入文件.xml，从CFD数据文件中筛选出所需的物理量数据，采用网格聚类算法将符合某一分区准则的体网格聚并到一个PSR反应器中，并顺序完成对应所有分区准则的网格聚并操作；计算反应器参数(体积V、压力p、温度T、组分等)和反应器间的质量流量。
反应器体积等于所包含的体网格i(i=1, 2, …, N)的体积之和，即V=∑Vi；反应器压力p和温度T则用体积加权平均计算得到：
完成上述分区及参数计算之后，将有关信息存储为CRN求解模块的输入文件.xml，此后视需要还可随时改动其中的一些参数值。
3) 运行CRN求解模块。读入上一步产生的输入文件，通过接口程序与开源反应动力学模拟软件Cantera连接并进行求解，即可获得各反应器的温度和详细组分信息以及燃烧室出口的NOx体积分数。为保证各反应器进出口的质量流量平衡，需计算反应器所有出口的流量分配比例并保持不变，再用实际输入的质量流量对各反应器间的质量流量进行校正。
以上算法程序实现了基于数值模拟结果的CRN预测NOx排放的整个自动求解过程。相比于求解NOx的后验方法，其计算耗时几乎可以忽略。除了手动完成Ansys FLUENT数据输出和输入文件的编辑等预备性操作外，对于反应器数量约为30个的CRN，上述算法程序在笔记本电脑上只需3 min即可完成自动分区并输出结果，计算效率非常高。
完成上述求解过程的关键在于分区准则，而分区准则是否合理和通用则取决于正确的分区策略，由此决定了反应器数量及计算效率、NOx预测准确性和方法的普适性。
2.2 分区策略和通用准则
正确的分区策略应体现出燃烧室的主要特征，并符合关于NOx生成与排放的燃烧基础理论。
现代DLN燃烧室有两个主要特征：1) 采用燃料分级的多喷嘴结构，燃料分级包括轴向分级、径向分级、多点式分级等多种形式，必然形成多个火炬或多重火焰结构，导致空间温度分布不均匀，如图 2所示。显然，在进行CRN分区时，必须兼顾由燃烧引起的温度轴向(沿流动方向)变化以及多重火焰结构导致的温度径向变化。
以上认识主要针对有火焰的燃烧反应区。实际上任何一个燃烧室均可按气流流动方向顺序分成纯流动区、燃烧反应区、后燃烧区3类大区，在它们之下再细分出用于CRN计算的小区，如表 2所示。
表 2 燃烧室整体分区 |
| 分区 | 大区类别 | ||
| 纯流动区 | 燃烧反应区 | 后燃烧区 | |
| CRN小区类别 | · 空气进口区 | · 回流区 | · 掺混区 |
| · 燃料进口区 | · 火焰锋面区 | · 过渡段区 | |
| · 预混区 | · 主火焰区 | ||
纯流动区只涉及空气、燃料(甲烷CH4)以及二者的混合物，并未发生燃烧，因此使用气流温度和组分的质量分数w作为分区准则较为合适，即wO2≥ 0.229为纯空气，wCH4≥0.999为纯甲烷燃料。T≤T2为空气/燃料的预混区，其中T2是燃烧室入口空气温度(由于入口燃料温度低于空气温度，在发生燃烧前，预混区温度不超过T2)。同时，考虑预混不均匀性, 使用预混区的平均当量比对预混区进行划分。本文所用燃烧室的预混区有两个：1) 中心喷嘴的预混通道，2) KL≥0.8的环形区(参见表 1)。
对火焰锋面区和主火焰区的分区准则相对复杂。本文研究发现，采用“燃尽率η”和“燃烧前的区平均当量比φav”最为合理。从物理上讲，燃尽率体现了燃烧释热程度和燃烧效率，决定了总体控制NOx生成的平均温度水平。由于本文所用燃烧室采用了径向燃料分级，即沿径向分成了中心区和外侧的环形区，且不同区域对NOx的贡献不同，故要分别计算两区的燃尽率，
式中：下标c代表中心区，r代表环形区。式(2)中的wCH4, in为该区进口的甲烷平均质量分数，wCH4为燃烧后的剩余值。相较于温度准则，使用燃尽率η能更准确地表达燃烧反应进度。
燃烧前的区平均当量比φav，实际上代表的是燃烧前的燃料质量分数。如果已知任一CRN轴向分区后得到的预混区、火焰锋面区和出火焰区的平均当量比，就能分辨出上述各区内部的燃料质量分数的高低差异。通过比较局部当量比φ与各区平均当量比的大小再细分子区，φ≥φav为燃料浓区，φ＜φav为燃料稀区，就能计算燃料/空气预混的不均匀度，进而得到预混不均匀度对NOx的关键影响。φ可用数值模拟算出的燃烧产物质量分数来计算，
综合以上研究，本文提出了对燃烧反应区的分区策略和通用准则：
1) 采用“两维度”分区法。先进行轴向/径向(纵剖面内)分区，然后完成横截面内分区。
2) 轴向/径向分区准则。可根据燃尽率η1=99%和η2=90%，确定火焰锋面的轴向/径向边界。
3) 横截面内分区准则。根据燃烧前的区平均当量比φav，可将已有的轴向/径向分区进一步分成稀、浓两区。
4) 环形区处于燃烧状态时(KL≤0.8)，只在区内按以上策略和准则分区；当处于熄灭状态时(KL≥0.8)，则在其下游燃烧区的外围进行分区。
上述分区策略和通用准则有如下优点：
1) 体现了燃烧空间各个维度上的不均匀性特征，充分反映了燃烧室结构特征和燃料/空气预混不均匀度对NOx生成的关键影响；
2) 在不同维度上采用不同的分区准则，以及先粗后细的分区顺序，体现了不同物理过程的特征；
3) 可区别对待环形区有无燃烧的不同工作模式，对预测部分负荷燃烧室的NOx排放有重要意义。
下面将通过数据验算证明本文提出的分区策略和通用准则(简称“η-φ准则”)的合理性和准确性。
3 应用结果与分析
本章采用第2章发展的CRN方法对表 1的燃烧工况进行计算，作为对比，还给出了文献中常用的温度准则分区方法的预测结果。
3.1 基于温度准则的NOx预测结果
表 3 基于温度准则与η-φ准则的CRN方法预测结果及对比 |
| 参数 | 工况编号 | |||||
| 1 | 2 | 3 | 4 | 5 | ||
| NOx体积分数(@15%O2) | 实测值/10-6 | 29.1 | 87.0 | 33.0 | 76.0 | |
| CFD计算值/10-6 | 3.4 | 32.2 | 65.2 | 2.6 | 15.6 | |
| CRN-T计算值/10-6 | 6.7 | 9.4 | 215.3 | 9.3 | 73.6 | |
| CRN (η-φ) 计算值/10-6 | 28.4 | 92.6 | 489.3 | 33.8 | 75.4 | |
| 误差 | δCFD/% | 88.3 | 63.0 | 92.1 | 82.7 | |
| δCRN-T/% | 77.1 | 89.1 | 71.7 | 3.2 | ||
| δCRN(η-φ)/% | 2.4 | 6.4 | 2.4 | 0.8 | ||
可见，基于温度准则的CRN方法完全不能适应工况的变化，除工况5外，其余工况的预测误差均超过70%，通用性很差。
参照文[18]中提高预测精度的惯常做法，本文进一步以工况2为例，研究了缩小温度间隔、增加分区即反应器数量的效果。由图 5可见，当温度间隔从初始的ΔT=200 K依次缩小到100、50、25 K，相应地，反应器数量从33增加到43、61、98个后，NOx预测值随之增大约5倍，但仍与实验值有很大差距，且增大趋势逐渐趋缓，明显形成一个系统误差。同时，计算耗时增加很快。这表明，单纯增加反应器数量的做法缺乏实用价值，除了如加密数值计算网格的作用一样，从量上对NOx预测分辨率有所改进外，对预测精度和计算耗时均有负面作用。从根本上讲，预测精度低是由于分区策略失当、分区准则不合理、未体现燃料/空气预混不均匀性的关键影响所致。
3.2 基于η-φ准则的NOx预测结果
采用η-φ准则对燃烧反应区进行CRN分区的具体步骤如下：
1) 径向分区。根据燃烧室径向特征，将燃烧区沿径向用直线简单分为3个区，即值班火焰区(7区)、中心火焰区(8、9区)、环形火焰区(10、11区)。
2) 轴向分区。用η1=99%和η2=90%依次确定火焰锋面区的两个轴向边界，见图 6a。
由图 7a可见，在采用η-φ准则下，轴向分区只要超过2个，NOx体积分数预测值就与实验值几乎完全吻合。如果不分区(Nzone=1)，预测值高达299.5×10-6，是实验值76×10-6的4倍。这是因为火焰温度变化剧烈，而CRN却采用的是完全均匀反应器(PSR)，二者反应状态差异甚大，使高温燃气停留时间的计算误差很大。由3.1节已知，分区个数并非越多越好。因此，本文用两个η值分出3个区就显得十分合理，既保证了预测精度，又留有足够的“安全裕度”。
由图 7b可见，在用η1=99%分出整个火焰锋面区的前提下，用于再分区的η2值就会使NOx预测值呈不规律的变化。这意味着，对火焰锋面区的再分区不能随意进行，而要使分出的两区体积(或区内的高温燃气停留时间)形成最佳匹配关系。η2=90% 正是最佳点，且所有其他工况亦如此，因此η2的取值也是合理的。
使用φ准则在横截面对燃烧空间进行分区，与轴向分区的情形类似，横截面内的分区个数和分区准则同样直接影响NOx的预测精度，因为分区必须体现出燃料/空气预混的不均匀特征(见图 6b和6c)。由图 8a可见，根据本文提出的η-φ准则，以中心区No.8区为研究对象，将其细分为多个子区，对分区个数进行检验，具体的划分方法为：分区个数为2时，以整区的平均当量比进行划分；分区个数为3时，则在2个分区的基础上，以浓区的平均当量比进一步划分，依此类推；将序号为No.8的小区再细分成多个子区后，本研究也得到与轴向分区类似的结论，即分区个数超过2的NOx预测值均与实验值吻合。可见，横截面内必须分区，且只需2~3个分区就完全满足精度要求。
本文进一步针对两个分区的情况，研究了φ取不同的分界值对NOx预测值的影响。由图 8b可见，当分界值取到待分区(即中心区和环形区)的平均当量比，再细分出稀(φ<φav)、浓(φ≥φav)两类子区时，预测值最接近实验值。除此之外的其他任何φ值都会产生更大误差。偏离φav越远，误差越大，且误差分布相对于φav还呈一定的对称性。例如，φ=0.8和0.7均偏离φav=0.75约6.7%，会使NOx的误差达到~10%。
η-φ准则同样适用于环形区有燃烧的情形(工况3)，其分区方案如图 9所示。
最终，基于η-φ准则的CRN方法预测结果如表 3所示。可见，预测精度有实质性的改善，最大误差仅为6.4%。由此证明本文提出的分区策略和通用准则是合理可行的。
4 结论
本文针对化学反应器网络法在燃烧室数值模拟中的应用开展研究，开发了CRN分区的自动生成和求解程序，提出了适用于某国产DLN燃烧室的分区策略和通用准则，并通过实验值进行检验，得到的主要结论如下：
1) 开发了CRN自动分区/求解算法和程序，采用XML数据格式对输入输出数据的格式进行标准化处理，提高了方法在不同燃烧室类型和工况条件下的通用性；
2) 提出了适用于燃气轮机燃烧室的CRN分区策略和通用准则，即η-φ准则，体现了燃烧室结构特征，符合基础燃烧理论的要求，物理意义明确，参数取值合理；
3) 基于η-φ分区准则的CRN方法适用于多负荷、不同燃烧工作模式以及燃烧室局部结构的变化等情况，从根本上提升了NOx排放的预测精度，最大误差不超过6.4%；
4) 发展了基于CFD-CRN耦合的完整准确描述湍流、燃烧和污染物生成的数值模拟方法，可以完全取代现有商用软件对NOx排放的后验式计算模型，大幅提高了计算效率和预测准确性。
