Generalization Capability of Convolutional Neural Networks for Progress Variable Variance and Reaction Rate Subgrid-Scale Modeling

Abstract

Deep learning has recently emerged as a successful approach to produce accurate subgrid-scale (SGS) models for Large Eddy Simulations (LES) in combustion. However, the ability of these models to generalize to configurations far from their training distribution is still mainly unexplored, thus impeding their application to practical configurations. In this work, a convolutional neural network (CNN) model for the progress-variable SGS variance field is trained on a canonical premixed turbulent flame and evaluated a priori on a significantly more complex slot burner jet flame. Despite the extensive differences between the two configurations, the CNN generalizes well and outperforms existing algebraic models. Conditions for this successful generalization are discussed, including the effect of the filter size and flame–turbulence interaction parameters. The CNN is then integrated into an analytical reaction rate closure relying on a single-step chemical source term formulation and a presumed beta PDF (probability density function) approach. The proposed closure is able to accurately recover filtered reaction rate values on both training and generalization flames.

1. Introduction

In LES of turbulent premixed combustion, the large-scale disparity between the domain size and the flame thickness typically does not allow adequate resolution of the internal flame front structure. As a consequence, one of the main challenges for LES is the modeling of the SGS reaction source term. To this end, recent studies have successfully used machine learning to train artificial neural networks as SGS models. They have been applied to the modeling of scalar dissipation rates [1], filtered density functions [2,3], wrinkling functions [4,5], and have also been used to predict LES filtered reaction rates directly [6,7] or by deconvolution [8]. These studies have consistently shown that neural networks can outperform physical algebraic models on test cases identical or similar to their training configuration.

A major unknown for all machine learning models is their ability to generalize to unseen configurations that can greatly deviate from the data they were trained on. Existing works in combustion operate in a supervised learning framework [9], where Direct Numerical Simulations (DNS) are used to produce paired input/output data for the training procedure. If machine learning models cannot guarantee some form of generalizability, their predictive ability will be limited to applications for which DNS data is available, which typically require petascale levels of computing power for realistic configurations [10].

This paper addresses this issue by exploring under what conditions a model trained on a canonical turbulent flame performs on a very different complex configuration. A CNN model is built for the SGS variance of the progress variable, a key quantity in many turbulent combustion closures [11,12,13,14]. Nikolaou et al. [8] had previously developed a CNN model for the variance, but it relied on the deconvolution of filtered quantities instead of direct computation of the variance, and was tested on cases similar to the training configuration.

The potential of the CNN model is illustrated by its integration in a new reaction rate closure. In a simplified single-step chemistry context, the mean LES reaction rate is analytically computed by a combination of a reaction source term formulation recently proposed by Pfitzner et al. [15,16] and a presumed beta probability density function (PDF) model, forming the so-called Pfitzner beta PDF CNN (PB-CNN) reaction rate closure. The CNN guarantees accurate predictions of the progress-variable variance, which parametrizes the beta PDF. Crucially, the analytical nature of the closure is expected to induce very fast evaluations, making this formulation a compelling alternative to tabulated chemistry methods that incur storage and computational costs [17,18,19].

The paper is organized as follows: the Pfitzner formulation for the reaction rate of a premixed flame is described in Section 2.1 while its coupling with a beta PDF model is presented in Section 2.2. Existing algebraic models for the progress variable variance are discussed in Section 2.3 before the rationale for using a CNN is discussed in Section 2.4. The training and generalization configurations are described in Section 3. The training procedure of the CNN is then detailed in Section 4. Results of the a priori evaluation of the CNN model and the PB-CNN closure are shown in Section 5 and Section 6.

2. Combustion Modeling Framework

2.1. Pfitzner Source Term

The derivation of the chemical source term formulation proposed by Pfitzner et al. [15,16] is summarized here. Assuming adiabaticity, constant pressure, unity Lewis number, and single-step chemistry, the evolution of a steady, one-dimensional (1D) fully premixed flame is fully determined by a progress variable $c=1-\frac{Y_F}{Y_F^f}$, where $Y_F$ is the fuel mass fraction and $Y_F^f$ is its value in the fresh gases. The transport equation for c reads [11]

$$\rho u\frac{\partial c}{\partial x}=\frac{\partial }{\partial x}\left(\rho D\frac{\partial c}{\partial x}\right)-\frac{{\dot{\omega }}_F}{Y_F^f},$$

(1)

where x, $\rho$, u, D, ${\dot{\omega }}_F$ are the spatial coordinate, density, velocity, diffusivity, and fuel reaction rate, respectively. It is convenient to introduce a reduced spatial variable

$$\xi \left(x\right)={\int }_0^x\frac{{\rho }_f{s}_L}{\rho \left(y\right)D\left(y\right)}\mathrm{d}y,$$

(2)

with ${\rho }_f$ the fresh gas density and $s_L$ the laminar flame speed, in order to convert Equation (1) to a nondimensional form

$$\frac{\partial c}{\partial \xi }=\frac{{\partial }^{2}c}{\partial {\xi }^2}-\frac{\rho D}{{\left({\rho }_f{s}_L\right)}^2{Y}_F^f}{\dot{\omega }}_F=\frac{{\partial }^{2}c}{\partial {\xi }^2}+\omega \left(c\right),$$

(3)

where the nondimensional source term $\omega \left(c\right)$ has been introduced. Pfitzner et al. [15,16] showed that $c_m\left(\xi \right)={[1+e^{-m\xi }]}^{-1/m}$ is a solution of Equation (3) parametrized by a model coefficient m. It is associated with the nondimensional source term

$${\omega }_m\left(c\right)=(m+1)(1-c^m)c^{m+1},$$

(4)

which is a polynomial alternative to the exponential Arrhenius formulation typically used to compute the reaction rate. A constant prefactor is defined as $\Lambda ={\left({\rho }_f{s}_L\right)}^2{Y}_F^f$ and the reaction rate ${\dot{\omega }}_F$ is linked to the nondimensional source term ${\omega }_m$ by

$${\dot{\omega }}_F=-\frac{\Lambda }{\rho D}{\omega }_m.$$

(5)

Equation (5) describes what will be referred to as the Pfitzner reaction rate term. Pfitzner et al. [15,16] chose the model parameter m to match the profile of an Arrhenius reaction rate. For all m values, the laminar flame speed obtained from the total fuel consumption rate is equal to the input $s_L$ value used in the prefactor $\Lambda$:

$$-\frac{1}{{\rho }_f{Y}_F^f}{\int }_{-\infty }^{\infty }{\dot{\omega }}_{F,m}\left(x\right)\mathrm{d}x=s_L.$$

(6)

On the other hand, the thermal flame thickness ${\delta }_L^0$ is inversely proportional to m and has the following analytical expression in $\xi$ space

$${\delta }_L^0=\frac{{(m+1)}^{\frac{m+1}{m}}}{m}.$$

(7)

The Pfitzner formulation can therefore be a substitute to any flame computed with Arrhenius chemistry, by matching independently the flame speed with the $\Lambda$ parameter and the flame thickness with the m parameter.

2.2. Pfitzner Beta PDF Closure for the Filtered Reaction Rate

In LES, any quantity of interest $\varphi$ is split into the sum of a mean and fluctuating component: $\varphi =\overline{\varphi }+{\varphi }^{\prime }$. Mean components are resolved by the numerical simulation on the LES mesh, while fluctuating components correspond to unresolved SGS activity that requires modeling. Mean components are related to the original quantity through a spatial filtering operation

$$\overline{\varphi }\left(x\right)={\int }_{V}G(x-x^{\prime })\varphi \left(x^{\prime }\right)\mathrm{d}{x}^{\prime },$$

(8)

where G is a low-pass filter with a cutoff scale equal to the LES mesh size. The filtered balance equations for reacting flows naturally introduce Favre-filtered (density-weighted) quantities $\tilde{\varphi }=\overline{\rho \varphi }/\overline{\varphi }$.

The evolution of all species mass fractions can be recovered from the transport equation of the Favre-filtered fuel mass fraction ${\tilde{Y}}_F$

$$\frac{\partial \overline{\rho }{\tilde{Y}}_F}{\partial t}+\nabla ·\left(\overline{\rho }\tilde{\mathit{u}}{\tilde{Y}}_F\right)=\nabla ·\left(\overline{\rho D\nabla Y_F}\right)-\nabla ·(\overline{\rho }\widetilde{\mathit{u}{Y}_F}-\overline{\rho }\tilde{\mathit{u}}{\tilde{Y}}_F)+{\overline{\dot{\omega }}}_F,$$

(9)

where ${\overline{\dot{\omega }}}_F$ is the filtered fuel reaction rate, which needs to be modeled to account for the unresolved SGS scales. To this end, PDF models rely on modeling the SGS distribution of c, noted as $p\left(c\right)$, to express the filtered reaction rate as the expected value of the DNS reaction rate

$${\overline{\dot{\omega }}}_F={\int }_0^{1}p\left(c\right){\dot{\omega }}_F\left(c\right)\mathrm{d}c.$$

(10)

A common approach [20,21,22,23,24] is to presume that $p\left(c\right)$ belongs to the family of beta PDFs that have the expression

$$p_{\beta }\left(c\right)=\frac{\Gamma (a+b)}{\Gamma \left(a\right)\Gamma \left(b\right)}{c}^{a-1}{(1-c)}^{b-1},$$

(11)

where $\Gamma$ is the gamma function, and a and b are two parameters linked to the mean $\overline{c}$ and variance $\overline{c^{\prime }{}^2}=\overline{c^2}-{\overline{c}}^2$ of the distribution by the relations

$$\begin{array}{cccc}\hfill a& =\overline{c}\left[\frac{\overline{c}(1-\overline{c})}{\overline{c^{\prime }{}^2}}-1\right],\hfill & \hfill b& =(1-\overline{c})\left[\frac{\overline{c}(1-\overline{c})}{\overline{c^{\prime }{}^2}}-1\right].\hfill \end{array}$$

(12)

In combination with the Pfitzner source term, the beta PDF approach is particularly interesting because the LES nondimensional reaction rate $\overline{{\left({\omega }_m\right)}_{\beta }}$ has an analytical expression [16]

$$\begin{array}{cc}\hfill \overline{{\left({\omega }_m\right)}_{\beta }}& ={\int }_0^1{p}_{\beta }\left(c\right){\omega }_m\left(c\right)\hfill \end{array}$$

(13)

$$\begin{array}{cc}& =(m+1)\frac{\Gamma (a+b)}{\Gamma \left(a\right)}\left(\frac{\Gamma (a+m+1)}{\Gamma (a+b+m+1)}-\frac{\Gamma (a+2m+1)}{\Gamma (a+b+2m+1)}\right).\hfill \end{array}$$

(14)

This is not the case with Arrhenius chemistry, and this property paves the way for computationally efficient evaluations of the filtered reaction rate. Using Equation (5), a similar formula can be derived for the fuel reaction rate under the simplifying assumption $\rho D\approx \overline{\rho D}=\overline{\rho }\tilde{D}$:

$$\begin{array}{cc}\hfill {\overline{\dot{\omega }}}_F={\int }_0^1{p}_{\beta }\left(c\right){\dot{\omega }}_F\mathrm{d}c& =\Lambda {\int }_0^1{p}_{\beta }\left(c\right)\frac{{\omega }_m\left(c\right)}{\rho \left(c\right)D\left(c\right)}\mathrm{d}c\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & =\frac{\Lambda }{\overline{\rho }\tilde{D}}\overline{{\left({\omega }_m\right)}_{\beta }}.\hfill \end{array}$$

(16)

with a Pfitzner beta PDF formulation, the turbulent LES reaction rate can be computed without needing on-the-fly integration of the PDF or tabulated chemistry.

2.3. Algebraic SGS Variance Models

To fully close the reaction rate, an SGS model is required for $\overline{c^{\prime }{}^2}$. Many closures rely on algebraic expressions based on the gradient of the resolved progress variable. One such model [21,25,26] is given by

$$\overline{c^{\prime }{}^2}=C_s{\Delta }^2{|\nabla \overline{c}|}^2.$$

(17)

where $\Delta$ is the local filter size and $C_s$ is the model parameter. It is important to note that this model was originally proposed for nonreactive scalars, and for reactive scalars such as the progress variable, it should be modified to account for the correlation of the scalar with its source term [27]. Nevertheless, it is still commonly found in the literature whenever an algebraic model for the progress variable is needed [25,26,28,29,30,31]. It is therefore used here as a basis of comparison for the CNN, which plays the same role as a traditional algebraic model. The simplest way of estimating $C_s$ is to use a global constant, and values such as $1/12$ [26,28], $0.18$ [25], or $0.5$ [30] have been proposed in the literature. However, a constant value is not expected to be accurate for all flow configurations and filter sizes [26,29], and the best model constant is a priori unknown. The model of Equation (17) with a constant $C_s$ parameter will be called the CST model.

$C_s$ can also be determined by a dynamic procedure [26,27,32,33] recalled in the following. A test filter of size $\widehat{\Delta }=2\Delta$ is used to generate test-filtered quantities, noted as $\widehat{\varphi }$. It is assumed that the model coefficient varies slowly in space, and that Equation (17) holds at the filter and test filter levels with the same model coefficient. Its expression at the standard filter level is filtered to the test-filter level

$$\widehat{\overline{c^2}}-\widehat{\overline{c}\overline{c}}=C_s{\Delta }^2\widehat{|\nabla \overline{c}{|}^2},$$

(18)

and independently expressed directly at the test-filter level

$$\widehat{\overline{c^2}}-\widehat{\overline{c}}\widehat{\overline{c}}=C_s{\widehat{\Delta }}^2{|\nabla \widehat{\overline{c}}|}^2.$$

(19)

After introducing a Leonard term defined as

$$\mathcal{L}=\widehat{\overline{c}\overline{c}}-\widehat{\overline{c}}\widehat{\overline{c}},$$

(20)

and a model term

$$\mathcal{M}={\widehat{\Delta }}^2{|\nabla \widehat{\overline{c}}|}^2-{\Delta }^2\widehat{|\nabla \overline{c}{|}^2},$$

(21)

the combination of Equations (18) and (19) leads to

$$\mathcal{L}=C_s\mathcal{M}.$$

(22)

As in [34], a least-squares procedure consisting of averaging over homogeneous directions of the flow is used to remedy stability issues. The final expression of the model coefficient is

$$C_s=\frac{\langle \mathcal{L}\mathcal{M}\rangle }{\langle {\mathcal{M}}^2\rangle },$$

(23)

where the averaging operation is denoted by the brackets $\langle ·\rangle$. The resulting model for $\overline{c^{\prime }{}^2}$ is called DYN.

2.4. Convolutional Neural Networks for SGS Variance Modeling

The simple algebraic model formulation of Equation (17) is driven by sound physical arguments, but also limits the family of $\overline{c^{\prime }{}^2}$ fields that these models can fit. As a consequence, algebraic models are often inaccurate [26,31,33] and transport equations for $\overline{c^{\prime }{}^2}$ are usually preferred [35,36,37] but they induce additional computations and require supplementary modeling. To combine computational efficiency and accuracy, $\overline{c^{\prime }{}^2}$ is modeled using a CNN, which can be trained to model nearly arbitrarily complex functions due to its large representational capacity [9]. Compared with fully-connected architectures used in some studies for SGS modeling [3,7], CNNs benefit from locality and translation invariance inductive biases that are helpful to process grid elements [38], and have been empirically observed to generalize better on reconstruction tasks in computer vision [39].

A CNN model for $\overline{c^{\prime }{}^2}$ is trained on a planar flame wrinkled by homogeneous isotropic turbulence with simplified Pfitzner chemistry and evaluated a priori on the DNS of a slot burner jet flame with skeletal chemistry. Since most thermophysical and chemical parameters are purposefully different from one configuration to the other, this is a difficult generalization test for the model. It is assumed that local values of $\overline{c^{\prime }{}^2}$ can be determined from three main elements:

• The profile of the $\overline{c}$ field in the neighboring flame brush;
• The amount of unresolved SGS scales;
• The effect of local turbulence on the SGS distribution of c.

If they are known by the model or kept similar in the training and generalization configuration, the model should be able to generalize properly. The essential nature of these elements is highlighted by their presence in the algebraic model of Equation (17), respectively, in the form of the terms $|\nabla \overline{c}|$, $\Delta$, and $C_s$. They therefore guide some of the key choices in the present work:

• A CNN is used because of its ability to accurately learn spatial patterns in an extended area around a location of interest (Section 4.1);
• the flame fronts seen in the training and evaluation contexts belong to the same turbulent combustion regime (Section 3.1 and Section 3.2);
• DNS snapshots are filtered and downsampled to a coarse grid with a well-chosen resolution (Section 4.2).

3. Training and Generalization Flow Configurations

3.1. Training Configuration: Planar Flame in Homogeneous Isotropic Turbulence

The training database is built from the DNS of a freely-propagating, initially planar flame immersed in decaying homogeneous isotropic turbulence (HIT), similar to [31,40,41,42]. The cubic domain is uniformly meshed with ${384}^3$ hexahedral elements and a mesh size $\Delta x=36 \mathsf{\mu }\mathrm{m}$. In the x direction, NSCBC inlet and outlet boundary conditions are prescribed [43], while the other boundary conditions are periodic (Figure 1). The flame front is initialized along the x direction using a DNS solution of a 1D laminar propane–air flame at stoichiometry, temperature $T=300 \mathrm{K}$, and pressure $P=1 \mathrm{bar}$. Chemistry is represented by a single irreversible reaction $C_3{H}_8+5O_2\to 3C{O}_2+4H_{2}O$ and the fuel reaction rate is computed using the Pfitzner DNS formulation of Equation (5). As key inputs needed in the formulation, the laminar flame speed $s_L=0.383 \mathrm{m}/\mathrm{s}$ and the model parameter $m=3.8$ are chosen to match the laminar flame speed and flame thickness ${\delta }_L=352 \mathsf{\mu }\mathrm{m}$ obtained with single-step Arrhenius chemistry [44]. The flame thickness is based on the gradient of c:

$${\delta }_L=\frac{1}{max\frac{\mathrm{d}c}{\mathrm{d}x}}.$$

(24)

The initial laminar flame front in the HIT is therefore resolved on approximately 10 mesh points.

In the fresh gases, turbulent velocity fluctuations are superimposed on the velocity field of the laminar flame. They are generated by a von Karman–Pao spectrum [45]:

$$E\left(k\right)=\alpha \frac{u^{\prime }{}^2}{k_e}\frac{{(k/k_e)}^4}{{[1+{(k/k_e)}^2]}^{17/6}}exp\left[-2{\left(\frac{k}{k_{\eta }}\right)}^2\right],$$

(25)

where $E\left(k\right)$, $u^{\prime }$, $k_e$, $k_{\eta }$ are the turbulent kinetic energy at the wavenumber k, turbulent fluctuation, wavenumber of the most energetic eddies, and wavenumber of the Kolmogorov scale $\eta$, respectively. $\alpha =1.453$ is a numerical constant. The integral length scale $l_t$ is related to $k_e$ through $l_t\approx 0.747/k_e$ [45]. $u^{\prime }$ and $k_e$ are chosen so that

$$\begin{array}{cccc}\hfill u^{\prime }/s_L& =10,\hfill & \hfill l_t/{\delta }_L& =4.8,\hfill \end{array}$$

(26)

leading to $u^{\prime }=3.83 \mathrm{m}/{\mathrm{s}}^{-1}$ and $k_e=423/{\mathrm{m}}^{-1}$. $k_{\eta }$ is obtained via the dissipation rate $ϵ\approx u^{\prime }{}^3/l_t$, which leads to $k_{\eta }=2\pi /\eta =2\pi {\nu }^{-3/4}{ϵ}^{1/4}=3.49\times 5 {\mathrm{m}}^{-1}$. The ratio of the mesh size to the Kolmogorov scale is $\Delta x/\eta =2$, so that the mesh is fine enough to resolve all the eddies of the initial spectrum.

Table 1. HIT initial parameters.

${\mathit{s}}_{\mathit{L}}$	${\mathit{\delta }}_{\mathit{L}}$	${\mathit{u}}^{\prime }$	${\mathit{l}}_{\mathit{t}}$	$\mathit{\eta }$	$\mathit{Da}$	$\mathit{Ka}$	${\mathit{Re}}_{\mathit{t}}$	${\mathit{Re}}_{\mathit{\lambda }}$	$\mathit{\tau }$
$0.383$$\mathrm{m}$/$\mathrm{s}$	352 $\mathsf{\mu }$$\mathrm{m}$	$3.83$$\mathrm{m}$/$\mathrm{s}$	$1.77$$\mathrm{m}$$\mathrm{m}$	18 $\mathsf{\mu }$$\mathrm{m}$	0.48	14	450	82	$0.46$$\mathrm{m}$$\mathrm{s}$

summarizes the main flame and turbulence parameters at initialization, including the Damköhler number $Da=(l_t/{\delta }_L)/(u^{\prime }/s_L)$, Karlovitz number $Ka={(u^{\prime }/s_L)}^{3/2}{(l_t/{\delta }_L)}^{-1/2}$, turbulent Reynolds number $R{e}_t=u^{\prime }{l}_t/\nu$, Taylor-scale Reynolds number $R{e}_{\lambda }=u^{\prime }{\lambda }_g/\nu$ based on the Taylor microscale ${\lambda }_g$, and eddy turnover time $\tau =l_t/u^{\prime }$.

The DNS is run using the fully compressible explicit code AVBP [46,47] with the TTGC Taylor–Galerkin finite element scheme of third-order accuracy in time and space [48]. The simulation is run for $1.84 \mathrm{m}\mathrm{s}$, corresponding to 4-eddy turnover times.

3.2. Generalization Configuration: R2 Slot Burner Jet Flame

The configuration used to assess the generalizability of the CNN is the DNS of the R2 slot jet flame studied by Luca et al. [49,50,51]. A fully premixed methane–air mixture at temperature $T=800 \mathrm{K}$, pressure $P=4 \mathrm{bar}$, equivalence ratio $\varphi =0.7$, and bulk velocity $U=100 \mathrm{m}/{\mathrm{s}}^{-1}$ is injected through a slot of width $H=1.2 \mathrm{m}\mathrm{m}$ and surrounded by a coflow of burnt gases. Chemistry is described by a skeletal chemical mechanism containing 16 species and 72 reactions [50]. The corresponding laminar flame has a speed of $s_L=1 \mathrm{m}/{\mathrm{s}}^{-1}$ and a thermal thickness ${\delta }_L=110 \mathrm{m}/{\mathrm{s}}^{-1}$. The domain dimensions are $28.8 \mathrm{m}\mathrm{m}$, $19.2 \mathrm{m}\mathrm{m}$, and $5.16 \mathrm{m}\mathrm{m}$ in the streamwise (x), crosswise (y), and spanwise (z) directions, respectively. It is uniformly meshed with a resolution $\Delta x=20 \mathsf{\mu }\mathrm{m}$. The progress-variable field of the fully developed jet flame is shown in Figure 2. 5 uncorrelated snapshots from instants where the jet is fully developed are used in this study, and form what will be called the generalization set.

3.3. Comments on the Differences and Similarities between the Two Configurations

The HIT and R2 flames feature very different domain geometries, flow dynamics, and simulation parameters (

Table 2. Main differences in the simulation parameters of the HIT and R2 DNS.

	HIT	R2
Fuel	$C_3{H}_8$	$C{H}_4$
Reactions	1	72
Species	5	16
$\varphi$	1	0.7
T	300 $\mathrm{K}$	800 $\mathrm{K}$
P	1 bar	4 bar
$s_L$	0.383 $\mathrm{m}/{\mathrm{s}}^{-1}$	1 $\mathrm{m}$/${\mathrm{s}}^{-1}$
${\delta }_L$	352 $\mathsf{\mu }$$\mathrm{m}$	85 $\mathsf{\mu }$$\mathrm{m}$
$\Delta x$	36 $\mathsf{\mu }$$\mathrm{m}$	20 $\mathsf{\mu }$$\mathrm{m}$

). In addition, turbulence in the flame brush is induced by homogeneous isotropic turbulence in one case and shear between the jet and the coflow in the other. Previous studies involving CNNs investigated minor parametric variations in the inlet condition [4], fuel species, and Karlovitz number [1], or turbulence intensity [5,8]. Substantial generalization was observed by Wan et al. [7] for a fully connected neural network trained as a surrogate model for the filtered reaction rate on a micromixing database.

The present work investigates the capacity of the CNN trained on the HIT to generalize well to flames in the same premixed turbulent combustion regime, characterized by the velocity and length ratios $u^{\prime }/s_L$, $l_t/{\delta }_L$ [11,52,53]. To this end, in Equation (26), the parameters of the initial turbulence spectrum of the HIT were chosen to match the $u^{\prime }/s_L$, $l_t/{\delta }_L$ values computed by Luca et al. [51] at the crosswise location of maximum heat release and $60\%$ of the flame length. They are considered to be representative of the turbulent combustion regime of the R2 flame.

Finally, three orders of magnitude separate the computational cost of the two DNS simulations. To simulate $0.9 \mathrm{m}\mathrm{s}$ of physical time, R2 requires 1 million CPU hours on an Intel Xeon Haswell-based supercomputer [49], while the HIT only needs 1600 CPU hours on a comparable Intel Xeon Skylake-based cluster. This highlights the greatly reduced computational demands of training an ML model on simple canonical configurations instead of full-scale realistic flames.

4. Machine Learning Framework

4.1. U-Net Architecture

The CNN trained in this study is a U-Net [54] with the same architecture used by Lapeyre et al. [4] to model SGS flame wrinkling. U-Nets are convolutional neural networks that are well-suited to field-to-field visual pattern recognition tasks, and are commonly used in computer vision for 2D and 3D image segmentation [54,55,56]. They rely on a succession of convolutional layers and downsampling operations to learn visual patterns in a wide area around the location of the prediction [4,9]. This is an appealing approach for combustion SGS modeling since the CNN is able to make local predictions based on the knowledge of the full spatial structure of the flame front.

The full architecture is detailed in Figure 3. It follows a fully convolutional, symmetrical encoder–decoder structure. Each stage is composed of two successions of

• A 3D convolution with a $3\times 3\times 3$ kernel with “same” padding;
• A batch normalization layer [57];
• A rectified linear unit (ReLU) nonlinear activation unit;
• A $2\times 2\times 2$ max-pooling downsampling operation in the encoder or a $2\times 2\times 2$ upsampling operation in the decoder.

The specificity of the U-Net architecture lies in the skip connections that concatenate the end of each encoder stage to the start of the corresponding decoder stage. Due to the downsampling operations, each stage of the encoder processes feature maps at decreasing resolutions. The combination of the upsampling and concatenation encourages each decoder stage to aggregate information from two separate spatial scales. Skip connections also allow the network to combine low-level features learned by the shallow layers of the encoder with the more complex, abstract features learned by the decoder, and accelerate training convergence [56]. Multiscale CNN architectures for field-to-field fluid predictions were successfully used by Lapeyre et al. [4] to model SGS flame wrinkling and Ajuria et al. [58] to solve the Poisson equation in incompressible flows. In the following, the U-Net model will simply be called the CNN. In the context of SGS modeling, this multiscale approach can be thought of as an extension of the dynamic procedure, which extracts information at a test-filtered scale to inform local predictions. Its ability to predict the entire output field at once is also attractive in terms of computational efficiency, as a single inference of the network is needed for the whole domain.

4.2. Data Preparation

The CNN is trained to process $\overline{c}$ input fields and predict associated $\overline{c^{\prime }{}^2}$ output fields. Following the LES framework, the filtered fields are represented on a grid of size $\Delta$ that is coarser than the DNS mesh. To generate these filtered coarse fields, 3D fields of c are first extracted from instantaneous snapshots of the DNS simulations. c is computed as

$$c_{HIT}=1-\frac{Y_F}{Y_F^f},$$

(27)

with the fresh gas fuel mass fraction $Y_F^f$.

The DNS progress variable fields are then filtered with a Gaussian filter

$$G\left(x\right)={\left(\frac{6}{\pi {\Delta }^2}\right)}^{3/2}exp\left[\frac{-6x^2}{{\Delta }^2}\right],$$

(28)

and downsampled to a coarse-mesh-structured grid with a resolution $\Delta$.

The filtering and coarsening procedure for HIT and R2 is illustrated in Figure 4. The value of the coarse mesh resolution for the HIT and R2 is chosen to ensure that the ratio $\Delta /{\delta }_L=0.81$ is the same for both cases. This is key to the ability of the CNN to generalize from one configuration to the other despite their different flame thickness and DNS mesh sizes. Since $\Delta /{\delta }_L$ is kept constant, the ratio of filtered laminar flame thickness to filter size ${\delta }_L^{filt}/\Delta$ is also approximately equal. It can be interpreted as the number of coarse mesh points on which the filtered flame brushes are resolved. This therefore guarantees that similar flame front patterns and gradients are found in the coarse HIT and R2, allowing the CNN to transfer from one case to the other. Additionally, the conservation of $\Delta /{\delta }_L$ establishes that these similar flame front gradients are unequivocably related to equivalent SGS activity, and equivalent values of $\overline{c^{\prime }{}^2}$ predicted by the CNN. In this regard, the filter sizes are chosen as equal to $288 \mathsf{\mu }\mathrm{m}$ and $70 \mathsf{\mu }\mathrm{m}$ for the HIT and R2, respectively.

4.3. Training Procedure

The training database is built from 46 instantaneous snapshots of the HIT DNS, extracted at regular intervals of $0.04 \mathrm{m}\mathrm{s}$. The first 38 constitute the training set for the CNN, while the next 4 form the validation set and the final 4 are kept as a hold-out test set for the results shown below. Instead of directly predicting $\overline{c^{\prime }{}^2}$, the CNN is trained to predict the unmixedness factor $\overline{c^{\prime }{}^2}/\left(\overline{c}(1-\overline{c})\right)$, a normalized version of the variance that was empirically found to lead to slightly better results. As a preprocessing step, the input and output fields are periodically padded in the y and z directions to enforce the periodicity boundary conditions of the DNS.

At each training iteration, 16 snapshots from the training set are sampled and, for each snapshot, 4 randomly cropped ${32}^3$ cubes are passed through random ${90}^{\circ }$ rotations and mirror operations before being added to the training batch. Gradient descent optimization is performed by an Adam optimizer [59] with an initial learning rate of $0.01$. A mean squared error loss function encourages the CNN outputs to match the ground truth $\overline{c^{\prime }{}^2}$ fields. Training was performed with the Tensorflow 2 Python library and stopped after 10,000 iterations when convergence was well-established. On an NVIDIA Tesla V100 GPU, this equates to one hour of training time. This training procedure only needs to be done once, after which the weights of the CNN are frozen and the trained model can predict $\overline{c^{\prime }{}^2}$ from a given $\overline{c}$ field in a few milliseconds.

5. A Priori Evaluation of the $\overline{c^{\prime }{}^2}$ Model

5.1. Evaluation on the HIT Test Set

The a priori performance of the CNN on the test set of the HIT is assessed to ensure that it has successfully learned to represent $\overline{c^{\prime }{}^2}$ on samples from the same distribution as the training database. The CNN model predictions are compared to ground truth $\overline{c^{\prime }{}^2}$ values computed from the filtered downsampled DNS, called true values, and to the CST model predictions. The constant coefficient $C_s=0.18$ [25] was chosen by picking the value in the literature [25,26,28,30] that best matched the true $\overline{c^{\prime }{}^2}$ values. For the DYN model, since there are no fully homogeneous directions, the averaging procedure is conducted on isocontours of $\overline{c}$. Another possibility is to perform averaging over the whole model, leading to a dynamically determined constant $C_s=0.17$, which is very close to the value chosen for the CST model.

Figure 5 shows sample slices of $\overline{c}$, true $\overline{c^{\prime }{}^2}$, and modeled $\overline{c^{\prime }{}^2}$ fields from the test set. The CNN matches the ground truth faithfully in all regions of the flow, while the gradient model does not seem as accurate.

Model predictions are aggregated on the whole test set and plotted against true values in Figure 6. The excellent performance of the CNN is confirmed by the tight spread of its predictions around the perfect model line, and the lack of bias of their linear trend. In comparison, CST and DYN model predictions display significantly higher variance and nonzero bias. Note that with prior knowledge of the ground truth, a $C_s$ value that leads to an unbiased CST model can be chosen but this would have little effect on the variance. The dynamic procedure proves to be effective at estimating a local model coefficient without any prior knowledge and leads to a slightly smaller error than the CST model.

Most of the data points plotted in Figure 6 are in regions of the flow where no flame front is present and using $\overline{c^{\prime }{}^2}$ for SGS closure is not important. The critical zones for $\overline{c^{\prime }{}^2}$ models are located in the reaction zone of the flame, identified as the regions where the laminar Pfitzner reaction rate is greater than 10% of its maximum value:

$$\mathcal{Z}=\{\overline{c}:{\dot{\omega }}_F\left(\overline{c}\right)\ge 0.1\times \underset{0\le c\le 1}{max}\left({\dot{\omega }}_F\left(c\right)\right)\}.$$

(29)

Hexplots conditioned on $\overline{c}\in \mathcal{Z}$ are plotted in Figure 7, and error values are also reported in

Table 3. HIT test set error values for $\overline{c^{\prime }{}^2}$.

Model	$\mathit{ϵ}$	$\mathit{ϵ},\overline{\mathit{c}}\in \mathcal{Z}$
CNN	0.060	0.054
CST	0.417	0.733
DYN	0.398	0.730

. It is interesting to note that the CNN performs even better in these critical regions, and is very close to a perfect model. On the other hand, the performances of the gradient models are significantly worse. Error values computed using a normalized mean square metric

$$ϵ=\frac{\sum {(y_{model}-y_{true})}^2}{\sum y_{true}^2}$$

(30)

are reported in Table 3 and confirm the excellent a priori performance of the CNN on the HIT test set.

5.2. Evaluation on the R2 Generalization Set

The various models are now assessed on the 5 snapshots of R2 composing the generalization set. For the DYN model, averaging is performed in the spanwise direction of the burner.

Figure 8 shows a sample slice of the $\tilde{c}$ field in a filtered snapshot of R2 and the predictions of the CNN, CST, and DYN model for the corresponding SGS variance. Once again, the value $C_s=0.18$ is retained for the CST model. All models reproduce qualitatively the evolution of the filtered DNS $\overline{c^{\prime }{}^2}$. The hexplots of Figure 9 and the error values reported in

Table 4. R2 error values for $\overline{c^{\prime }{}^2}$.

Model	$\mathit{ϵ}$	$\mathit{ϵ},\overline{\mathit{c}}\in \mathcal{Z}$
CNN	0.241	0.214
CST	0.371	0.331
DYN	0.358	0.322

reveal that the CNN performs the best out of all the models, which is a demonstration of its ability to transfer its predictive power from the HIT flame to the R2 flame. However, it is naturally less accurate than on the test of the HIT, which contained samples from the same distribution as its training set. In comparison, the CST and DYN models show zero bias but a higher variance. Unlike in the HIT, all the models perform similarly when evaluated on the whole domain or only the reaction zone $\mathcal{Z}$, as indicated by the error values in Table 4.

For a sample snapshot, the spatial distribution of model predictions is investigated by plotting $\overline{c^{\prime }{}^2}$ values averaged on transverse slices in Figure 10. The CNN recovers the filtered DNS values nearly perfectly, while the algebraic models consistently overestimate the variance.

This a priori evaluation of $\overline{c^{\prime }{}^2}$ models has shown that the CNN trained on the HIT configuration was able to generalize accurately to the R2 flame and outperform constant coefficient and dynamic algebraic models.

5.3. Discussion on the Conditions for Generalization

Between the HIT training/test sets and the R2 generalization set, three key ratios are conserved: $u^{\prime }/s_L$, $l_t/{\delta }_L$, and $\Delta /{\delta }_L$. It is important to understand the sensitivity of the generalization of the CNN to these ratios, as they could place strict limits on the practical applicability of the model.

The influence of $u^{\prime }/s_L$ and $l_t/{\delta }_L$ is analyzed by training a second CNN on the same HIT configuration with halved initial values of $u^{\prime }/s_L$, $l_t/{\delta }_L$, leading to lesser wrinkling in the turbulent flame front. When evaluated on R2, this second CNN had an error $ϵ=0.301$, which is $30\%$ higher than the reference CNN. Figure 11 shows that the effect of training on a weaker HIT is noticeable in the first $4 \mathrm{m}\mathrm{m}$ of the flame near the inlet. The second CNN underpredicts the SGS variance compared with the reference CNN and the filtered DNS. This is coherent with the fact that it is trained on a weaker HIT, which contains lower $\overline{c^{\prime }{}^2}$ values than the original configuration. This brief study seems to indicate that the choice of $u^{\prime }/s_L$, $l_t/{\delta }_L$ is indeed impactful on the generalization accuracy of the CNN, but more extensive work should be performed to fully understand the effect of these ratios.

The restriction to relying on a single value of $\Delta /{\delta }_L$ can be relaxed by training the CNN on a range of filter size values. An instance of the CNN was trained on a dataset comprising the collection of training HIT snapshots, which were filtered at $\Delta /{\delta }_L=0.8$, $1.2$ and $1.6$, for a total of 126 snapshots with 3 separate filter sizes. It was then tested on R2 for the same filter sizes, as well as unseen values $\Delta /{\delta }_L=1$, $1.4$, and 2. Figure 12 shows the evolution of the error made by the CNN and DYN models with $\Delta /{\delta }_L$. When $\Delta /{\delta }_L$ increases, the performance of the DYN model deteriorates as the SGS modeling task becomes more difficult. In contrast, the error made by the CNN is stable across all filter sizes, including unseen ones. This suggests that with proper training, the CNN can be accurate on a range of $\Delta /{\delta }_L$ values instead of a single one. Note that the error made by this CNN on $\Delta /{\delta }_L=0.8$ is higher than the one reported in Table 4 for the CNN solely trained on this filter size. This implies that a balance must be found between accuracy at a single filter size, and validity across a range of values. Interestingly, providing the $\Delta /{\delta }_L$ value as an additional input channel for the CNN did not improve the results. The rest of this work is presented with the original CNN trained with $\Delta /{\delta }_L=0.8$.

6. A Priori Evaluation of the PB-CNN Model for ${\overline{\dot{\omega }}}_F$

6.1. Evaluation on the HIT Test Set

In the previous section, the CNN was shown to be able to learn a model for $\overline{c^{\prime }{}^2}$, which is accurate on the test and generalization configurations. In the following, modeled values of $\overline{c^{\prime }{}^2}$ are incorporated in the presumed beta PDF approach detailed in Section 2.2 to form the PB-CNN model for ${\overline{\dot{\omega }}}_F$. Beta PDFs are known to sometimes be inaccurate models for the SGS distribution of c in premixed combustion [22,26,35]. Therefore, the presumed beta PDF approach needs to be justified independently from the model for $\overline{c^{\prime }{}^2}$. For the HIT, Figure 13 shows the conditional median, 5% and 95% quantiles plotted against $\overline{c}$ for the true filtered reaction rate, and the Pfitzner beta PDF reaction rate with true $\overline{c^{\prime }{}^2}$ values. The beta PDF appears to be a good approximation for the true SGS PDF, showing that the presumed beta PDF assumption is sound in the case of the HIT.

Given the accuracy of the CNN model, replacing the true $\overline{c^{\prime }{}^2}$ values with the CNN predictions does not noticeably affect the shape of the beta PDF. The isolated effect of the CNN model on the beta PDF reaction rate appears minimal in Figure 14, where the beta PDF reaction rates with CNN values of $\overline{c^{\prime }{}^2}$ are compared with beta PDF reaction rates with true values of $\overline{c^{\prime }{}^2}$.

Finally, the combined effect of the beta PDF assumption and the $\overline{c^{\prime }{}^2}$ model is shown in Figure 15. The high variance values predicted by the DYN model lead to underpredictions of the reaction rate when combined with the Pfitzner beta PDF closure. In contrast, using CNN predictions leads to a much closer fit to the true filtered reaction rate. Figure 16, where true and modeled reaction rates averaged on $yz$-planes are plotted against the streamwise x coordinate, illustrates how the quality of the CNN model for $\overline{c^{\prime }{}^2}$ leads to a good very estimate of fuel consumption rates using a beta PDF assumption.

6.2. Evaluation on the R2 Generalization Set

The PB-CNN closure is finally evaluated on R2. From the R2 coarse fields of $\overline{c}$, beta PDF reaction rates are computed according to Equation (15) and using either true $\overline{c^{\prime }{}^2}$ or predictions from the CNN. A reference 1D laminar flame matching the conditions of the R2 simulation is used to compute the prefactor $\Lambda ={\left({\rho }_f{s}_L\right)}^2{Y}_F^f$. The value for the Pfitzner model parameter $m=4.7$ was chosen so that the total fuel consumption rate of the reference 1D flame is matched by the Pfitzner reaction rate based on the c field of the reference 1D flame. As a result, by construction, both reaction rates would lead to the same total fuel consumption rate on a planar laminar flame. Figure 17 shows the close match of the resulting Pfitzner reaction rate to the fuel reaction rate of the skeletal mechanism.

Figure 18 compares the PB-CNN reaction rates with the true filtered values of R2. For clarity, only the points where the filtered DNS reaction rate is greater than 10% of its maximum value are plotted. This corresponds to the same condition as Equation (29). Compared with the HIT flame, the scatter plot still exhibits low bias, with a large portion of the predictions concentrated around the reference linear trend, but it also has a higher variance, which is emphasized by the logarithmic scale of the color bar. This indicates that some local values of the reaction rate are not well predicted by the PB-CNN model.

One particular group of mispredicted points is located on the x-axis of the plot and corresponds to regions in the reaction zone where the PB-CNN model erroneously predicts a reaction rate of exactly 0. Analyzing the reason for this misprediction is insightful to understand the behavior of the model in difficult edge cases. This group of points is characterized in Figure 19 by the distribution of their values of the true filtered reaction rate, filtered progress variable, true SGS variance, and CNN predictions for the variance. It appears that these points are located in the postflame region, where $c\approx 1$ and the skeletal mechanism still predicts some chemical activity. On the other hand, the SGS variance is extremely small and estimated by the CNN to be 0, causing the beta PDF-based reaction rate to also be exactly 0. This highlights a limitation of comparing a simplified reaction rate formulation relying on the SGS variance against a skeletal mechanism. However, these mispredictions are a rare occurrence since this group of points only constitutes 0.1% of the points in the reaction zone.

Despite discrepancies in local predictions of the reaction rate, the PB-CNN model provides correct estimates of spatially averaged reaction rates. This is important as it shows that the model would be useful in LES to recover key integral quantities of the burner such as the mean length of the flame. This is shown in Figure 20, which compares the streamwise evolution of the filtered DNS and PB-CNN reaction rates averaged on transverse slices. The PB-CNN closure recovers the correct averaged reaction rate in all regions of the flame. This leads to a fuel consumption rate in the whole domain within $5\%$ of the reference value. A posteriori simulations could be carried out in the future to assess if the PB-CNN reaction rate leads to an LES flame with the same length and total burning rate as the R2 DNS.

7. Conclusions

A deep convolutional neural network was trained to model the progress variable SGS variance $\overline{c^{\prime }{}^2}$ on a planar flame in homogeneous isotropic turbulence, and was shown to generalize accurately to the R2 complex slot burner jet flame in an a priori evaluation. Despite very different flow dynamics, thermophysical parameters, and chemistry, the conservation of the normalized filter size $\Delta /{\delta }_L$ and the turbulence–flame interaction ratios $u^{\prime }/s_L$, $l_t/{\delta }_L$ were assumed to be key to the generalizability of the CNN, and their influence on the results was discussed. This work indicates that CNNs can be trained on canonical simple cases and used in practical configurations. Future studies could focus on applying this methodology to different combustion models. The universal nature of the conditions for generalization that are exhibited here should also be investigated, as well as ways to integrate them into the machine learning model.

In an effort to combine computational efficiency and modeling accuracy, the PB-CNN closure for the mean LES fuel reaction rate was proposed. It combines the reaction source term formulation of Pfitzner et al. [15,16], a presumed beta PDF assumption, and the CNN model for $\overline{c^{\prime }{}^2}$ in an analytical formula for the reaction rate. It showed excellent a priori results on the R2 flame when compared with the burning rate produced by the original skeletal chemistry. In the future, an a posteriori evaluation would be a challenging test of the accuracy of the $\overline{c^{\prime }{}^2}$ model given by the CNN and the PB-CNN closure for ${\overline{\dot{\omega }}}_F$. This formulation would be especially well suited to simulating large-scale premixed combustion configurations, such as explosions in venting chambers. The resulting overpressure can be accurately captured using simplified chemistry, but remains a challenge at scale due to insufficient mesh resolution, prohibitive computational costs, and strong sensitivity to turbulent reaction rates [44,60].