湍流燃烧中的某些物理量，如切向应变率、位移速度、曲率等，对研究湍流和火焰间的相互作用具有重要的意义^[1-8]. 其中切向应变率是影响湍流火焰拉伸、褶皱特性的关键因素. Steinberg等^[7]通过实验测量，研究二维（two-dimensional，2D）应变率在湍流-预混火焰相互作用中对火焰拉伸的影响. Donbar等^[8]测量了二维应变率沿火焰褶皱轮廓的分布，研究火焰对应变率的响应. 切向应变率的准确计算需要速度和标量的三维（three-dimensional，3D）信息，在湍流燃烧实验测量中广泛使用的平面成像技术只能获得在测量平面上的低维信息. 目前，虽然交叉平面测量技术可以获得湍流燃烧中某些物理量的三维值^[9-10]，但是这些技术的成本较高，因此有必要建立湍流火焰切向应变率的低维度近似模型.

当前已有学者通过理论分析、实验测量、数值模拟等方法研究湍流燃烧中相关物理量的二维值和三维值之间的关系. 在实验测量方面，Zhang等^[11]通过实验测量获得的二维火焰表面密度，利用理论模型对三维统计平均值进行预测，预测结果较好. 在数值模拟方面，Hawkes等^[12-16]利用理论分析方法，获得湍流燃烧中切向应变率和位移速度等物理量的二维与三维统计平均值之间的关系，利用直接数值模拟验证了这种关系.

先前的研究仅对切向应变率及相关物理量的二维平均值与三维平均值之间的关系进行推导与验证. 本文提出切向应变率二维与三维概率密度函数（probability density function，PDF）间关系的理论模型. 利用机器学习方法，由相关低维量直接预测切向应变率的三维值. 建立包含人工神经网络和随机森林的2种模型，将实验能够直接测量得到的低维量作为模型的输入，将三维切向应变率作为输出，从而获得该物理量的三维预测结果. 所使用的样本基于直接数值模拟，样本来自不同湍流强度的燃烧算例. 研究的主要内容包括直接数值模拟数据库、数学背景和模型的介绍、模型预测结果的展现和讨论.

使用直接数值模拟（direct numerical simulation，DNS），研究三维统计上静止的平面预混火焰. 反应物由温度为300 K、当量比为0.7的贫燃甲烷/空气混合物组成. 相应的层流火焰速度 ${S_{\rm{L}}}$=0.19 m/s，层流火焰厚度 ${\delta _{\rm{L}}}$=0.66 mm，层流火焰时间 ${\tau _{\rm{L}}}$=3.47 ms. 由于流场的平均入流速度接近湍流火焰速度，火焰位置在统计上保持静止.

考虑3个具有代表性的DNS算例. 算例L、算例M、算例H分别对应弱、中、强湍流强度. 它们的模拟参数如表1所示. 表中， $u^{\prime} $为湍流脉动速度， ${l_{\rm{t}}}$为湍流积分长度尺度， ${\tau _{\rm{e}}}$为湍流时间尺度. 湍流积分长度尺度和层流火焰厚度之比 ${{{l_{\rm{t}}}} / {{\delta _{\rm{L}}}}}$=1. 湍流雷诺数 $Re = {\rm{ }}u^{\prime} {l_{\rm{t}}}/\nu $,其中 $\nu $为运动黏度. Karlovitz数 $Ka = {\tau _{\rm{L}}}/{\tau _\eta }$，其中 ${\tau _\eta }{\rm{ = }}{\left( {{{\nu {l_{\rm{t}}}} / u}{ ^{\prime2}}} \right)^{1/2}}$为Kolmogorov时间尺度. 在Peters^[17]的状态图中，算例L位于薄反应区，算例M和算例H位于破碎反应区.

计算区域为 ${L_x} \times {L_y} \times {L_z} = 28.6{\delta _{\rm{L}}} \times 4.8{\delta _{\rm{L}}} \times 4.8{\delta _{\rm{L}}}{\rm{ = }}$ $1.884\;{\rm{cm}} \times 0.314 \;{\rm{cm}}\times 0.314\;{\rm{cm}}$. 算例L和算例M的网格数为 ${N_x} \times {N_y} \times {N_z} = 768 \times 128 \times 128$，算例H的网格数为 ${N_x} \times {N_y} \times {N_z} = 1\;344 \times 224 \times 224$，且均为均匀网格. 网格精度能够求解所有的火焰和湍流尺度. 流向为x方向，2个横向分别为y和z. 流向采用无反射边界条件，横向采用周期性边界条件. 过程变量c在x-y平面切片的云图如图1所示. 过程变量是由氧气质量分数 ${w_{{{\rm{O}}_2}}} $定义的.

(1) $c = \frac{{{w_{{\rm{O}}_2}} - {w_{{{\rm{O}}_2},{\rm{b}}}}}}{{{w_{{{\rm{O}}_2},{\rm{u}}}} - {w_{{{\rm{O}}_2},{\rm{b}}}}}}.$

式中： ${w_{{{\rm{O}}_2},{\rm{u}}}} $和 ${w_{{{\rm{O}}_2},{\rm{b}}}} $分别为反应物和产物中氧气的质量分数.

模拟采用简化的甲烷/空气燃烧化学机理，对应的骨架机理由GRI-Mech3.0发展而来，包含268步基元反应和44种组分，其中28种组分通过DNS网格进行输运，其他16种组分设置为准稳态. 简化机理的验证可以参考文献[18].

以x-y平面作为测量平面，一个火焰微元所在的坐标系如图2所示. 图中，n为火焰面单位法向量，指向反应物；n_xy为n在xy平面上的投影 ${{n}}_{xy}^{\rm{p}} $再进行归一化后的单位向量； $\phi $为n_xy和n的夹角； $\theta $为y轴正方向和n_xy的夹角.

三维切向应变率的表达式为

(2) $\nabla \cdot {{u}} - {a_{{\rm{n}},3}}.$

测量平面上的二维切向应变率为

(3) ${\nabla _{xy}} \cdot {{{u}}_{xy}} - {a_{{\rm{n}},xy}}.$

式中：u为流体速度， ${a_{\rm{n}}} = {n_i}{n_j}{{ \cdot \partial {u_i}} / {\partial {x_j}}}$为法向应变率. 通过实验测量可以获得的二维量有二维法向量的2个分量 $ {n}_{x,xy}$、 $ {n}_{y,xy}$、二维速度梯度张量的分量 ${{\partial {u_i}} / {\partial {x_j}}}$（i，j=1，2）. 通过实验测量可以获得一些标量的值，如温度和部分基团的质量分数.

模型使用的样本选自DNS中火焰充分发展时刻和过程变量c=0.8的位置，即火焰面附近的数据. 算例L、算例M和算例H的数据量分别约为20万、30万和60万. 分别随机选取其中70%的样本作为人工神经网络模型的训练集及随机森林模型的袋内样本；选取30%的样本作为人工神经网络模型的测试集及随机森林模型的袋外样本，统称为测试样本. 模型的输入特征选取2.1节中实验能够直接测量的量，分别为H基团的质量分数（索引1）、二维法向量的分量 ${n_{x,xy}}{\text{、}}{n_{y,xy}}$（索引2、3）和二维速度梯度张量的分量 ${{\partial {u_i}} / {\partial {x_j}\left( {i,j = 1,2} \right)}}$（索引4~7）. 输出量为三维切向应变率. 输入特征中除了选取与输出量直接相关的量之外，还选取了H基团的质量分数. 这是由于H基团的分布能够在一定程度上反映火焰的三维特性^[19].

因为这些输入特性具有不同的量纲和量级，在训练过程开始之前，将输入特征全部进行归一化：

(4) ${x_{i,{\rm{nor}}}} = \frac{{{x_{i,{\rm{real}}}} - x_{i,{\rm{real}}}^{\min }}}{{x_{i,{\rm{real}}}^{\max } - x_{i,{\rm{real}}}^{\min }}}.$

式中： ${x_{i,{\rm{nor}}}}$为进行归一化之后第i组的输入特征， ${x_{i,{\rm{real}}}}$为归一化之前第i组输入特征的真实值， $x_{i,{\rm{real}}}^{\max }$和 $x_{i,{\rm{real}}}^{\min }$分别为该组输入特征的最大值和最小值. 进行归一化之后，所有输入特征的范围均为0~1.0，这有利于训练过程的收敛. 对所有的输出量分别进行反归一化：

(5) ${y_{i,{\rm{pred}}}} = {y_{i,{\rm{M}}}} \left( {y_{i,{\rm{real}}}^{\max } - y_{i,{\rm{real}}}^{\min }} \right) + y_{i,{\rm{real}}}^{\min }.$

式中： ${y_{i,{\rm{M}}}}$为模型输出的预测值， $y_{i,{\rm{real}}}^{\max }$和 $y_{i,{\rm{real}}}^{\min }$分别为训练样本中对应目标量的最大值和最小值， ${y_{i,{\rm{pred}}}}$为模型最终预测值.

在火焰切平面上定义标准正交坐标系，如图2所示，其中t₁与火焰面相切且位于测量平面上，t₂与火焰面相切且与t₁垂直. 在坐标系t₁-t₂-n中，三维法向应变率等于速度梯度张量在火焰法向n上的分量，即

(6) ${a_{{\rm{n}},3}} = {{\partial {u_{\rm{n}}}}}/{{\partial {x_{\rm{n}}}}}.$

三维切向应变率^[14]为

(7) ${a_{{\rm{t}},3}} = \nabla \cdot {{u}} - \frac{{\partial {u_{\rm{n}}}}}{{\partial {x_{\rm{n}}}}} = \frac{{\partial {u_{t_1}}}}{{\partial {x_{t_1}}}} + \frac{{\partial {u_{{{t}}_2}}}}{{\partial {x_{t_2}}}}.$

此时位于测量平面内的2D切向应变率为

(8) ${a_{t_1,xy}} = {{\partial {u_{t_1}}}}/{{\partial {x_{t_1}}}}.$

位于平面外的2D切向应变率为

(9) ${a_{t_2,xy}} = {{\partial {u_{{{t}}_2}}}}/{{\partial {x_{{{t}}_2}}}}.$

在各向同性假设的基础上，进一步提出2种假设. 第1种假设为 ${a_{{t_1},xy}} = {a_{{t_2},xy}} $；第2种假设为 ${a_{{t_1},xy}} $与 ${a_{{t_2},xy}} $统计上相互独立. 当 ${a_{{t_1},xy}} = {a_{{t_2},xy}} $时，三维切向应变率a_t,3=2a_t,xy. a_t,3的概率分布函数为

(10) $\begin{split} & {F_{{a_{{\rm{t}},3}}}}\left( y \right) = 1 - P\left( {{a_{{\rm{t}},3}} \leqslant y} \right) = \\ & 1 - P\left( {{a_{{\rm{t}},2}} \leqslant y/2} \right) = 1 - \int_{ - \infty }^{\frac{1}{2}y} {{p_{{a_{{\rm{t}},xy}}}}\left( x \right){\rm{d}}x} . \end{split} $

式中： ${p_{{a_{{\rm{t}},xy}}}}$为二维切向应变率的PDF. 对a_t,3的概率分布函数求导，得到切向应变率2D与3D的PDF的关系为

(11) ${p_{{a_{\rm{t}}},{\rm{A}}1}}\left( {{a_{\rm{t}}}} \right) = \frac{1}{2}{p_{{a_{\rm{t}}},xy}}\left( {\frac{1}{2}{a_{\rm{t}}}} \right).$

式中：下标‘A1’表示由假设1获得的PDF. 当 ${a_{{t_1},xy}} $与 ${a_{{t_2},xy}} $统计上相互独立时，它们和的PDF为它们各自PDF的卷积，即

(12) ${p_{{a_{\rm{t}}},{\rm{A}}2}}\left( {{a_{\rm{t}}}} \right) = \int_{ - \infty }^\infty {{p_{{a_{t_1}},xy}}\left( {{a_{t_1}}} \right){p_{{a_{t_2}},xy}}\left( {{a_{\rm{t}}} - {a_{t_1}}} \right){\rm{d}}{a_{t_1}}} .$

由于各向同性假设， ${a_{{t_1},xy}} $与 ${a_{{t_2},xy}} $统计上是等价的，即P( ${a_{{t_1},xy}} $)=P( ${a_{{t_2},xy}} $)，因此

(13) ${p_{{a_{\rm{t}}},{\rm{A}}2}}\left( {{a_{\rm{t}}}} \right) = \int_{ - \infty }^\infty {{p_{{a_{t_1}},xy}}\left( {{a_{t_1}}} \right){p_{{a_{t_1}},xy}}\left( {{a_{\rm{t}}} - {a_{t_1}}} \right){\rm{d}}{a_{t_1}}} .$

全连接前馈人工神经网络（artificial neural network，ANN），也称为多层感知机. 采用的ANN模型的拓扑结构如图3（a）所示，包含输入层、单个隐藏层和输出层. 该模型收到输入特征后，对权值矩阵和偏置矩阵进行随机初始化，通过正向传播，对输入进行预测. 模型会引入评估预测值与真实值差距的损失函数，即均方误差 $\sigma $，表达式为

(14) $\sigma = {{\displaystyle\sum\limits_{i = 1}^N {{{\left( {{y_{i,{\rm{M}}}} - {y_{i,{\rm{real}}}}} \right)}^2}} }}\Big/{N}.$

直至该误差达到理想值时，训练才会结束，否则将进行反向传播过程，权值矩阵和偏置矩阵通过输出层的信息进行修正，使误差不断下降^[20]. 模型引入L2正则化方法，以缓解模型过拟合的现象^[21]. 输出量与输入特征的非线性关系可由下式表示：

(15) ${y_{i,{\rm{M}}}} = {f_2}\left[ {{{{W}}_2}{f_1}\left( {{{{W}}_1}{{{x}}} + {{{B}}_1}} \right) + {{{B}}_2}} \right].$

式中：x=[x₁, x₂, ···, x_n]；f_i为各层神经元的激活函数，采用的激活函数为ReLU函数，表达式为

(16) $f\left( {{m}} \right) = \left\{ \begin{array}{l} {{m}},{{m}} > {{0}}; \\ {{0}},{{m}} \leqslant {{0}}. \end{array} \right.$

随机森林（random forest，RF）是包含多棵决策树、用于分类、回归目标的集成学习方法^[22]. 该方法的原理是在训练时建立大量不相关的决策树，各自独立地学习和预测，将每棵树的预测结果进行综合判断，因此最终结果一般优于单棵树的结果，可以有效地避免过拟合现象. 如图3（b）所示，对于含有M个特征的N个样本，从原始样本中随机抽样生成与决策树个数相等的样本集，生成的样本集分别用于每棵树的学习过程. 对于每个样本集，将样本分为袋内（in bag）和袋外（out of bag，OOB）2种类型，袋内样本进行训练，袋外样本对模型进行测试. 最终的预测结果由每棵树进行加权后平均得到.

影响ANN模型和RF模型预测精度的主要模型超参数是隐藏层神经元个数和决策树个数. 采用试错法对它们开展超参数调优，调优的评估参数分别是测试集和袋外样本的均方误差，具体的调优结果在3章中展示和讨论.

理论模型预测的PDF结果如图4所示. 图中，‘A1’表示假设1，‘A2’表示假设2. 切向应变速率的范围随着Ka的增加而变大，这是由强湍流火焰中产生的大应变率引起的. A2预测的PDF与DNS实际的PDF非常吻合，相比之下，A1的预测精度低于A2.

ANN模型隐藏层神经元个数和RF模型决策树个数的选择采用试错法进行确定. 根据相关研究的经验^[23-24]，试错法选取神经元个数的测试范围为2~99，决策树个数的测试范围为2~49，分别选择测试集和袋外样本的均方误差作为试错法的评价指标. 如图5所示为隐藏层神经元个数N_neu、决策树个数N_tree与对应 $\sigma $的关系. 可以看出，随着隐藏层神经元个数或决策树个数的增加，均方误差的总体趋势是先下降，后虽稍有波动但整体平稳. 考虑到模型精度和计算成本，算例L、算例M和算例H最终选择的隐藏层神经元个数分别为88、61、51，决策树个数分别为22、34、33.

为了比较2种机器学习模型的性能，选择一种更好的模型. 为了分析模型的预测值与真实值之间的误差，判断模型是否能够满足要求，定义以下3种评估参数进行衡量：Pearson相关系数R，用来衡量预测值与真实值的相关性；决定系数R²，称作模型优度，在机器学习中，多采用该系数来评价模型的性能；平均绝对百分误差MAPE，能够反映预测值与真实值之间的相对误差. 它们的表达式如下：

(17) $ \begin{split} &R =\\ &\frac{{N\displaystyle\sum\limits_{i = 1}^N {{y_{i,{\rm{pred}}}}{y_{i,{\rm{real}}}}} \!-\! \displaystyle\sum\limits_{i = 1}^N {{y_{i,{\rm{pred}}}}} \displaystyle\sum\limits_{i = 1}^N {{y_{i.{\rm{real}}}}} }}{{\sqrt {N\displaystyle\sum\limits_{i = 1}^N {{y_{i,{\rm{pred}}}}^2} \!-\! {{\left( {\!\displaystyle\sum\limits_{i = 1}^N {{y_{i,{\rm{pred}}}}} } \right)}^2}} \sqrt {\!N\displaystyle\sum\limits_{i = 1}^N {{y_{i,{\rm{real}}}}^2} \!-\! {{\left( {\displaystyle\sum\limits_{i = 1}^N {{y_{i,{\rm{real}}}}} } \right)}^2}} }}, \end{split} $

(18) ${R^2} = 1{\rm{ - }}\frac{{\displaystyle\sum\limits_{i = 1}^N {{{\left( {{y_{i,{\rm{pred}}}} - {y_{i,{\rm{real}}}}} \right)}^2}} }}{{\displaystyle\sum\limits_{i = 1}^N {{{\left( {{y_{i,{\rm{real}}}} - {{\overline y }_{i,{\rm{real}}}}} \right)}^2}} }},$

(19) ${\rm{MAPE}} =\frac{1}{N}\sum\limits_{i = 1}^N {\frac{{\left| {{y_{i,{\rm{pred}}}} - {y_{i,{\rm{real}}}}} \right|}}{{\left| {{{\overline y }_{_{i,{\rm{real}}}}}} \right|}}} \times 100{\text{%}} .$

式中： ${\overline y _{i,{\rm{real}}}}$为真实值的平均值.

利用2种模型测试样本切向应变率三维预测值与真实值之间的散点图，如图6所示. 图中， ${a_{\rm{t}}}{\tau _{\rm{L}}}$为切向应变率的三维真实值， ${a_{{\rm{t,p}}}}{\tau _{\rm{L}}}$为模型预测值. 如表2所示为模型的R、R²和MAPE. 从图6、表2可以看出，3种算例的相关系数R分别从ANN模型的0.908、0.818、0.792增大到了RF模型的0.997、0.984、0.981. RF模型的相关系数已经非常接近1，这在湍流研究中是可以接受的误差. R²分别从ANN模型的0.822、0.669、0.623增大到了RF模型的0.994、0.968、0.960，增幅随着湍流强度的增强而增大. 这是因为随着湍流强度的增强，样本的无序性和随机性增大，使得ANN模型产生了很强的过拟合. RF模型的学习过程具有更强的稳定性，减弱了噪声数据对结果的影响，能够预测强湍流条件下的样本. 2个模型的平均绝对百分误差MAPE是类似的情况. RF模型的性能优于ANN模型，这在许多机器学习的研究中是类似的^[25-26].

由于湍流运动的无序性与随机性，在统计学上对预测结果进行分析. RF模型的性能更好，仅基于该模型进行分析. 利用RF模型预测得到的PDF、二维和三维真实PDF如图7所示. 可以看出，利用RF模型预测的PDF与三维真实的PDF非常接近. 与图4相比，RF模型预测得到的PDF比理论模型预测得到的PDF与真实PDF吻合更好，说明在统计层面，RF模型的精度达到了非常高的程度.

综合多种评估参数可以发现，RF模型的性能达到了非常高的水平，通过RF模型由相关低维量来预测火焰三维切向应变率是可行的.

在湍流燃烧中，二维火焰面法向量、二维流体速度梯度张量以及组分质量分数与切向应变率有着直接或间接的联系，影响切向应变率二维和三维之间的关系. 分析相关低维量对切向应变率的具体影响程度，即特征重要性，为某一特征在所有特征中贡献的百分占比. 从3.2节的分析可知，RF模型的性能优于ANN模型，因此基于RF模型进行分析. 在RF模型中，特征重要性是通过判断每个特征在每棵树上的贡献，再取加权平均值.

各输入特征对切向应变率的重要性J如图8所示. 图中，横坐标为各输入特征的索引. 可以看出，除二维法向量及二维速度梯度张量这些与切向应变率直接相关的量之外，H基团质量分数对三维切向应变率有较大影响，且随着湍流强度的增强，这种影响逐渐增强. 这说明在较强湍流条件下某些基团的质量分数和切向应变率有深刻的关联，能够在一定程度上反映它在第3个方向上的特性.

本文提出切向应变率2D与3D的PDF之间的代数关系，其中假设各主应变率相互独立时模型预测的PDF与DNS的真实PDF吻合较好. 构建ANN模型和RF模型，基于直接数值模拟数据库，由相关低维量对三维切向应变率进行预测. 通过对不同模型预测结果的分析和对比发现，RF模型预测值与真实值的R和R²均更大，其中R² > 0.96；平均绝对百分误差MAPE更小，MAPE < 17%. RF模型的性能很好，能够通过相关的低维量预测三维切向应变率.