对燃气轮机燃烧室进行数值模拟是设计过程中的重要环节。然而，全尺寸燃烧室的数值模拟具有空间尺度跨度大，以及流动、传质传热、燃烧和污染物生成等多个复杂过程的耦合特性，需要使用多种数学模型，并且具有网格复杂、收敛难度大、耗时长的特点^[1]。例如，对于一个约600万网格的算例，在960 CPU(单核主频2.5 GHz)的条件下，需要运行3 h才能获得收敛解。在燃烧室设计过程中，需要对设计和结构参数等进行调整和反复迭代，因此需要进行大规模数值验算，会导致大量的CFD算例。目前，国际上发展出一种代理模型(surrogate model, SM)方法，采用在有限数量的计算流体动力学(computational fluid dynamics, CFD)算例中的插值算法，可迅速、完整地得出多个中间工况的流场参数，取代原本需要的三维CFD算例^[2-4]。对于同样的算例，SM在数秒内即可完成，与CFD相比，显著提高了工作效率。

SM是一种基于一定量样本数据来描述输入变量和输出数据之间关系且计算代价较小的数学模型^[5-6]。SM的构建方法包含3个基本步骤：样本选取、数据降维和工况插值^[2]。

样本选取即形成一个算例集，如同插值所需的一套离散数据点集。样本的选取质量对SM的精度有较大影响，因此需要在变量空间中均匀地分布样本点，以提高精度。目前有以下几种选取方法，例如正交法、均匀设计法、中心复合法，以及最常用的Latin超立方抽样方法(Latin hypercube sampling, LHS)等。其中，LHS是一种在多维变量空间中的分层抽样方法，随机性强，可以根据需要任意选择样本点的数量^[7-8]。

在完成样本算例的CFD计算后，如果对每一个网格点上的流场参数进行插值，计算量会非常大，因此需要对流场数据进行降维。本征正交分解(proper orthogonal decomposition, POD)是一种常用的方法，通过POD提取出一组代表流场主要特征的基函数矩阵(基函数个数不超过样本数)，使用这组基函数以级数叠加的形式来计算流场数据，从而将计算维数降低到与基函数的个数一致，以大幅提高计算效率^{[3-4, 9]}。

插值方法包括多项式法、径向基函数法、神经网络法和Kriging插值法等。其中Kriging插值法是最常用的最优线性无偏估计方法之一，被广泛应用于地质学、环境科学、空间遥感和计算统计等领域^[10-13]。

目前，SM已广泛应用于多个领域。例如，国内外分别构建了飞行器气动外形参数与气动特性之间的SM，包括升阻比、升力系数等。SM也应用于火电厂空冷凝汽器内部流场、温度场与环境风速的关系建模，以及发动机旋流喷嘴冷态流场的建模等^{[2, 7, 14]}。然而，文献研究中只涉及单一流体流动，并且未见有关燃烧室SM的研究，而燃烧室燃烧是一种多过程耦合的复杂系统。例如，三维高湍流度的旋流流动、燃料射流、燃料/空气混合过程、传热、燃烧化学反应和污染物生成等。因此，从设计角度还需对SM的结构进行改变。同时，多参数和多过程的耦合可能会对数据降维、插值精度等产生影响，进而影响SM的精度。

本文先构建燃烧室冷态物理混合过程的SM，检验SM方法的可行性，为后续构建包含燃烧过程、更复杂条件下的模型奠定基础。为了避免燃烧室复杂过程带来的干扰，拟先做冷态且仅关注状态参数的变化，即不涉及燃烧反应过程，不考虑结构参数的变化。

为构建适用于燃烧室冷态流动混合过程的SM的方法，本文对构建过程中的样本选取、数据降维和工况插值等关键步骤进行了研究，并制定了应用准则，为后续研究更复杂的燃烧反应过程和结构变化的模型构建奠定基础。

1 燃烧室数值模拟方法 1.1 基本结构和工作模式

燃烧室的基本结构如图 1所示，图中数值的单位均为mm。由导流衬套、头部喷嘴组件、火焰筒和过渡段组成。燃烧室采用径向分级的结构，分为中心喷嘴和环形区；环形区内包含5个沿周向分布的自带旋流器的喷嘴，且在环形区外壁有2排主燃孔，在火焰筒下游周向均布3个掺混孔。

中心喷嘴的工作模式为：空气从喷嘴外环进入预混通道，燃料则直接喷入，与空气相互掺混，形成预混可燃混合物，并通过轴向旋流器进入燃烧区，轴线处设有扩散燃烧的值班喷嘴。环形区有2种工作模式：当燃烧室相对负荷N_e≤0.8时，发生燃烧；当N_e>0.8时，仅进行预混。此前，该燃烧室已有详细的数值模拟和试验验证，SM在制作样本算例时，使用原来的网格^[15]。

1.2 计算域选取

本文选取从中心喷嘴入口到过渡段入口截面的三维燃烧空间(见图 1)作为计算域，主要考虑以下几方面：1) 中心喷嘴是燃烧室的典型核心部件，在气体动力学方面有双重作用，一是通过旋流器在其下游的燃烧空间形成用于稳定火焰的回流区；二是实现燃料喷射及其与空气的预混过程，可反映与高速流动、高湍流度和燃料/空气的预混均匀度等方面有关的问题。因此，SM能展示燃烧室的特殊性和典型性。2) 在过渡段入口截面之前，回流区已结束，流场较均匀，因此可作为计算域的出口。3) 暂不考虑环形区和掺混孔，简化不必要的影响因素。4) 可同时考察和评价SM再现横截面和纵剖面的CFD结果的能力。

中心喷嘴出口截面将燃烧室分为2个区域，前面为流动混合区，后面为燃烧反应区。如果需要详细模拟燃烧的过程，可通过SM提供燃烧区的入口边界条件。同时，纵剖面可以反映心回流区、角回流区等主要流场特征。因此，选取中心喷嘴出口截面和纵剖面作为SM的考察截面。

1.3 数值模拟设置

1) 计算网格。

本文采用原先工作的网格文件^[15]。在整体燃烧室网格基础上去掉中心喷嘴前的导流衬套、环形区和掺混孔的网格，仅保留计算域。网格类型为非结构四面体网格，网格数为182.3万，如图 2所示。

2) 数学模型和边界条件。

除了通用的质量、动量守恒方程和组分扩散方程，湍流模型使用Realizable k-ε模型，边界条件为：空气入口使用速度进口边界条件；截取该燃烧室中心喷嘴入口截面的速度分布，并设置进口压力；燃料入口使用质量流量入口边界条件；燃烧室出口使用自由出流边界条件；壁面为无滑移边界。

1.4 SM输入参数的选择

选取在冷态条件下参数变化较大的燃料量G_f和进口压力p₂这2个参数。燃料量代表燃烧室的负荷变化，直接体现中心喷嘴中的预混状况和出口处的燃料浓度分布。参数G_f的具体数值范围为0.040 0~0.120 0 kg/s，p₂为1 330~1 510 kPa，对应的N_e范围为0.5~1.0。

1.5 流场和混合过程的数值模拟结果

图 3展示了经CFD计算得到的燃烧室纵剖面流线图。图中的计算工况为：空气流量中G_a= 1.766 kg/s，温度T₂=665.5 K，压力p₂=1 349 kPa；燃料流量中G_f=0.074 0 kg/s，温度T_f=288.0 K。在横截面突扩处形成角回流区；由于旋流作用在喷嘴下游形成中心回流区，且回流区前端进入中心喷嘴内，在喷嘴的出口截面处中心区域的轴向速度为负值，外围区域的轴向速度为正值，如图 4所示。由于该喷嘴上侧某旋流叶片有探视孔，气流通过探视孔以直流的形式进入燃烧室，造成旋流强度下降，导致上侧的回流区长度较短，因此中心回流区上下并不对称。

中心喷嘴的燃料浓度分布，如图 5a所示。X_CH4为摩尔分数；C_CH4为甲烷浓度，kmol/m³。燃料通过中心通道进入燃烧室，其中一部分作为值班燃料直接进入燃烧室；另一部分通过4个燃料喷嘴喷出，并在预混通道内与空气混合。由于预混通道的方式是直流，且预混距离较短，喷嘴出口截面处的燃料/空气预混均匀度并不理想。图 5b和5c展示了CH₄的浓度分布，图中的外围区分别形成了4个间隔分布的高浓度区，对应燃料喷管的位置。

2 结果与分析 2.1 SM构建流程

构建流程如图 6所示，使用LHS方法确定样本算例；在完成CFD计算后，通过POD进行数据降维和Kriging模型进行工况插值，完成初步构建；最后以验证算例的CFD计算结果为基准，对精度进行检验；如果误差较大，则需重新构建。

2.2 样本选取

首先，确定样本个数N_samp，主要取决于输入变量的个数N_var。为保证SM的精度，应满足的准则可表示如下：

$ {N_{{\rm{samp }}}} \ge 5{N_{{\rm{var }}}}. $

(1)

结合1.4节中冷态燃烧室的输入变量G_f和p₂，可得N_var=2，由式(1)确定样本数N_samp=12。

样本选取的具体方法为：把每个变量的变化范围等分为N_samp个区间，在每个区间内随机生成一个数值，这样每个变量都得到N_samp个抽样值，然后将各变量的抽样值随机组合，即可得到N_samp个样本。本文生成12个样本，同时也构成了12个CFD算例，将G_f按大小进行排序，样本的参数条件如表 1所示。

在表 1的样本中选择4个CFD验证算例。算例取自该燃烧室负荷N_e范围为0.5~1.0的真实运行参数，用于检验SM的误差，如表 2所示。表 1和2中的工况只有参数G_f、p₂发生了改变，其他参数保持不变，与1.5节中的工况相同。

2.3 数据降维处理

使用POD进行降维，首先将各个样本的流场数据，包括：速度、压力、燃料浓度、湍动能等，按照计算网格节点的顺序进行排列，构成m×1的矩阵Y _t，表示如下：

$ {\mathit{\boldsymbol{Y}}_t} = \left[ {\begin{array}{*{20}{c}} {{Y_{t, 1}}}\\ {{Y_{t, 2}}}\\ \vdots \\ {{Y_{t, m}}} \end{array}} \right]. $

(2)

其中：Y_t为第t个样本的流场数据，m为样本的网格节点总数，各样本的网格节点和排列顺序应保持一致。

将各样本的流场数据进行排列构成矩阵Y，表示如下：

$ \mathit{\boldsymbol{Y}} = \left[ {{\mathit{\boldsymbol{Y}}_1}, {\mathit{\boldsymbol{Y}}_2}, \cdots , {\mathit{\boldsymbol{Y}}_{{N_{{\rm{samp }}}}}}} \right]. $

(3)

其中：Y_Nsamp为样本数据，N_samp为样本总数。

按照Galerkin加权残值法完成基函数的求解，具体算法为：计算样本矩阵Y的內积；计算矩阵Q的特征值λ_i及对应的特征向量a_i(i=1, 2, …, N_samp)，并按照特征值大小进行排列；计算各特征值λ_i所对应的函数φ_i，并进行归一化；计算基函数ϕ_i的系数β_{t, i}。计算过程分别表示如下：

$ \begin{array}{c} \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}} = {\mathit{\boldsymbol{Y}}^{\rm{T}}}\mathit{\boldsymbol{Y}}, }\\ {{\lambda _1} > {\lambda _2} > \cdots > {\lambda _{{N_{{\rm{samp}}}}}}, }\\ {{\mathit{\boldsymbol{\varphi }}_i} = \mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{a}}_i}, }\\ {{\mathit{\boldsymbol{\phi}} _i} = \frac{{{\mathit{\boldsymbol{\varphi }}_i}}}{{\left| {{\mathit{\boldsymbol{\varphi }}_i}} \right|}}, } \end{array}\\ \left[ {\begin{array}{*{20}{l}} {{\beta _{t, 1}}}\\ {{\beta _{t, 2}}}\\ \vdots \\ {{\beta _{t, {N_{{\rm{samp}}}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\phi}} _1^{\rm{T}}}\\ {\mathit{\boldsymbol{\phi}} _2^{\rm{T}}}\\ \vdots \\ {{\mathit{\boldsymbol{\phi}} _{N_{{\rm{samp}}}}^{\rm{T}}}} \end{array}} \right]{\mathit{\boldsymbol{Y}}_t}. \end{array} $

(4)

通过上述计算，将流场数据Y_t分解为由一组基函数ϕ_i级数相加的形式，表示如下：

$ {\mathit{\boldsymbol{Y}}_t} = \sum\limits_{i = 1}^{{N_{{\rm{samp }}}}} {{\beta _{t, i}}} {\mathit{\boldsymbol{\phi}} _i}. $

(5)

其中：基函数ϕ_i为m×1的矩阵，β_{t, i}为对应的系数，并满足以下条件。

1) 基函数是一组标准正交基，表示如下：

$ \left( {{\mathit{\boldsymbol{\phi}} _i}, {\mathit{\boldsymbol{\phi}} _j}} \right) = \left\{ {\begin{array}{*{20}{l}} {1, }&{i = j;}\\ {0, }&{i \ne j.} \end{array}} \right. $

(6)

2) 样本矩阵Y_t在基函数上的投影最大，表示如下：

$ \max \sum\limits_{t = 1}^{{N_{{\rm{samp }}}}} {\sum\limits_{i = 1}^{{N_{{\rm{samp }}}}} {\left( {{\mathit{\boldsymbol{Y}}_t}, {\mathit{\boldsymbol{\phi}} _i}} \right)} } . $

(7)

在后续的计算中，通过对系数β_{t, i}的插值完成SM的构建，计算维度由样本的网格点数降低为基函数的个数。

2.4 Kriging插值

在对流场数据进行POD后，需对系数β_{t, i}和SM的输入参数(即G_f和p₂)进行插值，从而得到非样本点上的系数值，并计算得到各参数的分布，表示如下：

$ {\mathit{\boldsymbol{Y}}_{{\rm{SM}}}} = \sum\limits_{i = 1}^{{N_{{\rm{samp }}}}} {{\beta _{{\rm{SM}}, i}}} {\mathit{\boldsymbol{\phi}} _i}. $

(8)

其中：Y_SM为SM计算得到的参数矩阵，β_{SM, i}为非样本点上的系数值。

Kriging插值法是指，假设对一个目标函数y进行插值，共有n个数据点，定义x₁、x₂、…、x_n为自变量，y₁、y₂、…、y_n分别为对应的函数值。

该方法将y(x)表述成函数平均值μ与一个随机过程ε(x)的组合，表示如下：

$ y(x) = \mu + \varepsilon (x). $

(9)

随机过程ε(x)应满足无偏性和数据点x_l、x_m之间的协方差关系，分别表示如下：

$ E[\varepsilon (x)] = 0, $

(10)

$ {\mathop{\rm corr}\nolimits} \left[ {\varepsilon \left( {{x_l}} \right), \varepsilon \left( {{x_m}} \right)} \right] = \exp \left[ { - \sum\limits_{k = 1}^d {{\theta _k}} {{\left( {{x_{kl}} - {x_{km}}} \right)}^2}} \right]. $

(11)

其中：θ_k为核函数，x_kl和x_km分别为变量x_l和x_m在第k维的数值，d为x_l和x_m变量的维变。

假设使用Kriging插值法求得在预测点x_new处的函数值为y_new，最优解是指该点处预测值的误差方差最小，即var(y_new-y₀)最小，其中y₀为真实值，使用最大似然估计方法进行推导，表示如下：

$ {y_{{\rm{new}}}} = \hat \mu + \mathit{\boldsymbol{r}}_{{\rm{new}}}^{\rm{T}}{\mathit{\boldsymbol{R}}^{ - 1}}\left( {\mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{A}}_n}\hat \mu } \right), $

(12)

$ \hat \mu = \frac{{\mathit{\boldsymbol{A}}_n^{\rm{T}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{y}}}}{{\mathit{\boldsymbol{A}}_n^{\rm{T}}{\mathit{\boldsymbol{R}}^{ - 1}}{\mathit{\boldsymbol{A}}_n}}}, $

(13)

$ \mathit{\boldsymbol{y}} = {\left[ {{y_1}, {y_2}, \cdots , {y_n}} \right]^{\rm{T}}}{\rm{. }} $

(14)

其中：A_n为n×1的矩阵，矩阵每一个元素均是1；R为n×n的协方差矩阵，元素R (i, j)=corr[ε(x_l), ε(x_m)]，如式(11)所示；r_new为n×1的矩阵，元素分别是预测点x_new与其他数据点的协方差，表示如下：

$ \mathit{\boldsymbol{r}}(l, 1) = {\mathop{\rm corr}\nolimits} \left[ {\varepsilon \left( {{x_{{\rm{new}}}}} \right), \varepsilon \left( {{x_l}} \right)} \right]. $

(15)

Kriging插值方法不仅可以输出预测值，还可以给出预测精度，即预测值方差s²，可用于对不确定度的量化^[16]，表示如下：

$ {s^2} = {\hat \sigma ^2}\left[ {1 - \mathit{\boldsymbol{r}}_{{\rm{new}}}^{\rm{T}}{\mathit{\boldsymbol{R}}^{ - 1}}{\mathit{\boldsymbol{r}}_{{\rm{new}}}} + \frac{{{{\left( {1 - \mathit{\boldsymbol{A}}_n^{\rm{T}}{\mathit{\boldsymbol{R}}^{ - 1}}{\mathit{\boldsymbol{A}}_n}} \right)}^2}}}{{\mathit{\boldsymbol{A}}_n^{\rm{T}}{\mathit{\boldsymbol{R}}^{ - 1}}{\mathit{\boldsymbol{A}}_n}}}} \right], $

(16)

$ {\hat \sigma ^2} = \frac{{{{\left( {\mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{A}}_n}\hat \mu } \right)}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\left( {\mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{A}}_n}\hat \mu } \right)}}{n}. $

(17)

其中${\hat \sigma ^2}$为Kriging过程方差。

2.5 SM计算结果

采用上述方法构建SM，并计算出各验证算例(见表 2)在喷嘴出口横截面和纵剖面上各参数的分布。验证算例1(YZ₁)的计算结果如图 7—11所示，v_a为轴向速度，v_r为切向速度。由图可知，SM完全再现了CFD计算得到的横截面和纵剖面的轴向速度和切向速度分布。此外，其他所有参数，如k、p、CH₄浓度场等的分布均与CFD计算结果相似。

2.6 误差分析

为了综合评价经SM计算的变量大小和空间分布与CFD结果之间的误差δ，需要采用向量运算进行表征^[2]，表示如下：

$ \delta = \frac{{\left| {{\mathit{\boldsymbol{Y}}_{{\rm{SM}}}} - {\mathit{\boldsymbol{Y}}_{{\rm{CFD}}}}} \right|}}{{\left| {{\mathit{\boldsymbol{Y}}_{{\rm{CFD}}}}} \right|}} \times 100\% . $

(18)

其中：Y_SM和Y_CFD分别为SM和CFD结果，二者均为N_samp×1的矩阵，格式与式(2)相同。

SM计算结果的误差，如表 3所示，表中v_t为径向速度，可知SM在横截面和纵剖面内的计算结果的误差均≤4%，远低于国际上10%的平均水平^[2-4]，表明尽管燃烧室内的气流速度和湍流度较高，存在旋流和回流区，且存在不同气体组分的混合，在这种情况下，本文构建的SM仍有较强的处理能力和较高的准确性。

2.7 SM构建方法的应用准则

1) 采用内插值的方式。

在SM中，外插和内插是指预测点工况与样本工况区间的关系。如果预测点工况在样本工况区间之外，称为外插；如果预测点工况在样本工况区间之内，称为内插。

研究表明SM是一种数据驱动(data driven，DD)方法，由于样本集不包含参数范围外的信息，采用外插会显著降低SM的精度。本文的算例也证明了这一点，外插使SM算出的CH₄浓度场的误差成倍增加，如图 12所示。图中以表 1中算例11的CH₄浓度场作为验证算例，从表 1和2的算例中选取8个不同的组合构建SM，并且这些组合的下边界值为G_{f, down}=0.043 3 kg/s，上边界值G_{f, up}由0.117 0 kg/s变化至0.074 0 kg/s，可表示如下：

$ \Delta {G_{\rm{f}}} = {G_{{\rm{f}}, 11}} - {G_{{\rm{f}}, {\rm{up}}}}. $

(19)

其中ΔG_f为样本算例11的边界值与样本组合上边界值的差值。ΔG_f为负值时，表明插值方式为内插，此时以表 1中的算例12为上边界，其余所有组合均使算例11构成了外插；ΔG_f越大，表明验证算例外插得越远。

2) 样本数目足够多。

从数学原理上讲，用于插值的样本数量越多，SM计算结果的精度就越高。使用验证算例1的喷嘴出口截面的CH₄场进行验证，当样本数量N_samp由12降至7时，误差约增大2倍，如图 13所示。

3) 基函数个数选取合理。

在POD分解中，得到一组基函数ϕ_i，分别对应样本內积矩阵的特征值λ_i(按从大到小排列)。通常前几个特征值较大，其对应的基函数也代表了样本分布的主要特征。因此，在实际应用中，通常不会使用所有的基函数，而只选用前几阶基函数就可以达到较高的精度，以提高计算效率。定义前q个基函数的能量系数ξ^[2]，表示如下：

$ \xi = \frac{{\sum\limits_{j = 1}^{{N_{{\rm{samp }}}}} {\sum\limits_{i = 1}^q {{{\left( {{\beta _{j, i}}{\mathit{\boldsymbol{\phi}} _i}} \right)}^2}} } }}{{\sum\limits_{j = 1}^{{N_{{\rm{samp }}}}} {\sum\limits_{i = 1}^{{N_{{\rm{samp }}}}} {{{\left( {{\beta _{j, i}}{\mathit{\boldsymbol{\phi}} _i}} \right)}^2}} } }} \times 100\% . $

(20)

表 4以喷嘴出口界面的CH₄浓度场为例，展示了各验证算例的误差、能量系数与基函数个数的关系，可知：1) 第一阶基函数的能量占比极高，说明第一阶基函数的分布代表了各样本的主要特征结构；2) 随着基函数个数的增加，能量系数随之增加，SM的精度也逐渐提高。然而，后面几个基函数的能量占比较低，即使增加基函数的个数，误差下降的速度仍越来越慢。由表 4可知，当能量系数达到99.99%时，相较于使用全部基函数，仅使用3个基函数的误差增加不超过0.05%。

3 结论

本文提出了一种构建燃烧室冷态流动混合过程的SM的方法，使用了LHS方法、本征正交分解、Kriging插值等算法，并证明了其合理可行。在只改变状态参数的情况下，SM具有处理冷态高速、高湍流、强旋流动和燃料/空气掺混的能力，其精度高于国际平均水平。本文还制定了SM构建方法的应用准则，包括：1) 使用内插值的方式；2) 遵循样本个数和基函数个数选取的约束准则。为后续更复杂、包含燃烧反应过程以及结构变化的燃烧室SM构建方法研究奠定基础。