2.1. Problem statement

The simulation of turbulent flows at high Reynolds numbers cannot generally be performed with sufficient resolution to represent all the scales of the flow. LES aims to simulate only a part of the turbulent spectrum to achieve a reasonable computational cost. Scales are cut by spatial filtering, and unresolved scales are modeled using so-called subfilter-scale models. In the LES framework, the Favre-filtered transport equation for the mass fraction of the $ {k}^{\mathrm{th}} $ species $ {Y}_k $ follows Eq. (1.2). As discussed in the introduction, modeling the filtered chemical source terms in LES of reacting flows is extremely challenging because chemistry largely occurs at the subfilter-scale level in most cases. This is further complicated in lean premixed $ {\mathrm{H}}_2 $ -air mixtures because of thermodiffusive effects, causing additional flame wrinkling and pronounced fluctuations in reaction rates along the flame front. A precise description of the preferential and differential diffusion, turbulence wrinkling, and their interactions (Berger et al., Reference Berger, Attili and Pitsch2022c) at the subfilter-scale level is required.

In this work, we focus on the development of a data-driven model for the filtered chemical source term of $ {\mathrm{H}}_2 $ in lean premixed turbulent combustion. ML is used to discover structural relations between resolved scalar quantities and subfilter-scale turbulence-chemistry interactions. We develop this method for the $ {\mathrm{H}}_2 $ filtered “burning rate” $ \overline{\dot{\omega}} $ , defined as

(2.1) $$ \overline{\dot{\omega}}=-{\overline{\dot{\omega}}}_{{\mathrm{H}}_2}, $$

with the aim of establishing an ML-based modeling framework and a proof of concept for $ {\mathrm{H}}_2 $ turbulent combustion. The final implementation will depend on the given LES modeling framework. For example, the ML model could replace the tabulation method proposed in Lapenna et al. (Reference Lapenna, Berger, Attili, Lamioni, Fogla, Pitsch and Creta2021, Reference Lapenna, Remiddi, Molinaro, Indelicato and Creta2024) for the chemical source term of the progress variable. This is what has been done in Seltz et al. (Reference Seltz, Domingo, Vervisch and Nikolaou2019) with direct mapping of the closure terms for a methane-air flame. The ML framework could also be based on a database with a one-step chemistry (Millán-Merino and Boivin, Reference Millán-Merino and Boivin2024; Schiavone et al., Reference Schiavone, Detomaso, Torresi and Laera2024) to close the LES multi-species transport equations.

The basic ingredients to describe turbulent flames remain the quantities introduced for laminar flame analysis: the progress variable $ c $ and the mixture fraction $ \xi $ . The progress variable equals zero in fresh gas, and one in burnt gas, and varies monotonically during combustion. Mixture fraction $ \xi $ or the equivalence ratio $ \phi $ is used to inform the ML-based model of the composition of the mixture, which impacts burning rates. The mixture fraction is based on Bilger’s definition (Bilger et al., Reference Bilger, Stårner and Kee1990).

(2.2) $$ \xi =\frac{\beta -{\beta}^o}{\beta^f-{\beta}^o}\hskip0.24em \mathrm{with}\hskip0.24em \;\beta =\frac{Z_{\mathrm{H}}}{2{W}_{\mathrm{H}}}-\frac{Z_{\mathrm{O}}}{W_{\mathrm{O}}}, $$

where $ {Z}_m $ and $ {W}_m $ are the elemental mass fraction and atomic mass of the element $ m $ . The superscripts $ o $ and $ f $ denote oxidizer and fuel stream. The filtered progress variable $ \tilde{c} $ is defined as

(2.3) $$ \tilde{c}=\frac{Y_{{\mathrm{H}}_2}^u\left(\tilde{\xi}\right)-{\tilde{Y}}_{{\mathrm{H}}_2}}{Y_{{\mathrm{H}}_2}^u\left(\tilde{\xi}\right)-{Y}_{{\mathrm{H}}_2}^b\left(\tilde{\xi}\right)}=\frac{\tilde{\xi}-{\tilde{Y}}_{{\mathrm{H}}_2}}{\tilde{\xi}-\max \left(0,\frac{\tilde{\xi}-{\xi}_s}{1-{\xi}_s}\right)}, $$

where superscripts $ u $ and $ b $ denote values in unburnt and burnt gas and $ {\xi}_s $ is the mixture fraction value at stoichiometry. We have arbitrarily chosen to base the progress variable on the $ {\mathrm{H}}_2 $ mass fraction. It is shown in Supplementary Material that it is also possible to base the progress variable on the $ {\mathrm{H}}_2 $ O mass fraction without any impact on the training and prediction accuracy of the model. In this framework, the challenge is now to model $ \overline{\dot{\omega}} $ as a function of $ \tilde{c} $ and $ \tilde{\xi} $ fields.

2.2. Machine learning framework

The aim of the ML-based model is to approximate the function

(2.4) $$ \overline{\dot{\omega}}=f\left(\tilde{c},\tilde{\xi}\hskip0.22em \mathrm{or}\hskip0.22em \tilde{\phi},\dots \right). $$

The evolution of the reaction rate depends on the thermochemical state described by the progress variable and mixture fraction, but also on the flow. Indeed, the level of turbulence, flame curvature, and strain, have a major influence on the reaction rate (Poinsot and Veynante, Reference Poinsot and Veynante2011). In this work, we propose to include the topological information by means of spatial convolutions using a CNN. The CNN is found to perform better when the equivalence ratio $ \tilde{\phi} $ is used instead of the mixture fraction $ \tilde{\xi} $ . The filtered equivalence ratio can be approximated as

(2.5) $$ \tilde{\phi}=\frac{\tilde{\xi}\left(1-{\xi}_s\right)}{\xi_s\left(1-\tilde{\xi}\right)}. $$

The CNN takes as input the field of $ \tilde{c} $ and $ \tilde{\phi} $ over a domain $ \Omega $ , performs successive convolutions and operations, and outputs an approximation of $ \overline{\dot{\omega}} $ over the same domain $ \Omega $ . The function $ {f}_{\mathrm{CNN}} $ therefore reads

(2.6) $$ {\overline{\dot{\boldsymbol{\omega}}}}^{\mathrm{NN}}={f}_{\mathrm{CNN}}\left(\tilde{\boldsymbol{c}},\tilde{\boldsymbol{\phi}}\right)\approx \overline{\dot{\boldsymbol{\omega}}}. $$

CNNs were designed to learn spatial hierarchies of features, from low- to high-level patterns. They are therefore well fitted for recognizing complex flame topology and learning associated burning rates. Furthermore, the convolutions enable to train models on large inputs via parameter sharing. The application of CNNs to the modeling of flame features has already been successfully performed in Lapeyre et al. (Reference Lapeyre, Misdariis, Cazard and Poinsot2018, Reference Lapeyre, Misdariis, Cazard, Veynante and Poinsot2019), Seltz et al. (Reference Seltz, Domingo, Vervisch and Nikolaou2019) and Nikolaou et al. (Reference Nikolaou, Chrysostomou, Vervisch and Cant2019) for fuel-air mixture without significant preferential diffusion effects, unlike the $ {\mathrm{H}}_2 $ -air mixtures considered in this article.

The CNN architecture used is inspired by the U-net architecture originally developed for image segmentation (Ronneberger et al., Reference Ronneberger, Fischer, Brox, Navab, Hornegger, Wells and Frangi2015) (Figure 1). However, the output activation function is replaced by a linear activation function for the regression task considered in this work. A similar architecture for the approximation of flame surface density has already been used successfully (Lapeyre et al., Reference Lapeyre, Misdariis, Cazard and Poinsot2018, Reference Lapeyre, Misdariis, Cazard, Veynante and Poinsot2019). Using $ 64 $ channels in the first multi-channel feature map instead of $ 32 $ in Lapeyre et al. (Reference Lapeyre, Misdariis, Cazard and Poinsot2018, Reference Lapeyre, Misdariis, Cazard, Veynante and Poinsot2019) increased the accuracy of the approximation of $ \overline{\dot{\omega}} $ (Eq. (2.6)) in the present work. In total, the network comprises $ \mathrm{5,424,193} $ trainable parameters.

Figure 1. U-net-type architecture is used in the present work. Each black box corresponds to a multi-channel feature map. The number of channels is denoted on top of the box. The input sample has three channels: $ \tilde{c} $ , $ \tilde{\phi} $ and $ {\delta}_L^0/{\delta}_L^1 $ . The ratio $ {\delta}_L^0/{\delta}_L^1 $ is the inverse of the laminar flame thickening due to filtering (Section 2.3.2). The original size of the cubic sample is $ {N}^3 $ . It is then reduced to $ {\left(N/2\right)}^3 $ and $ {\left(N/4\right)}^3 $ during the contracting path before going back to the original size during the expansive path. Gray boxes represent copied feature maps. The arrows denote the different operations.

2.3. Data generation and training strategy

The data generation strategy consists of filtering and downsampling data from DNS to emulate LES data with known filtered quantities. The configuration chosen to generate training data is a slot burner (Figure 2) at constant pressure $ P=1 $ atm and fresh gas temperature $ {T}^u=300 $ K, as in Lapeyre et al. (Reference Lapeyre, Misdariis, Cazard and Poinsot2018, Reference Lapeyre, Misdariis, Cazard, Veynante and Poinsot2019). The physical domain consists of a central inlet where a premixed $ {\mathrm{H}}_2 $ -air mixture flows at a bulk velocity $ {U}_b=24\;\mathrm{m}/\mathrm{s} $ with velocity fluctuation $ {u}^{\prime }=2.4\;\mathrm{m}/\mathrm{s} $ , surrounded by two laminar coflows where burnt gas flows at a bulk velocity $ {U}_c=3.6\;\mathrm{m}/\mathrm{s} $ . The injection of turbulence at the central inlet corresponds to homogeneous and isotropic turbulence using a Passot-Pouquet turbulence spectrum (Passot and Pouquet, Reference Passot and Pouquet1987) with an integral length scale $ {l}_t=2 $ mm. The domain is rectangular with periodic boundary conditions in the $ z $ -direction. Adiabatic walls are specified in the $ y $ -direction. Both inlets and outlets are specified in the $ x $ -direction. Exact dimensions are annotated on Figure 2. The configuration is very similar to that used in Berger et al. (Reference Berger, Attili and Pitsch2022c), Coulon et al. (Reference Coulon, Gaucherand, Xing, Laera, Lapeyre and Poinsot2023), and Gaucherand et al. (Reference Gaucherand, Laera, Schulze-Netzer and Poinsot2024a) to study $ {\mathrm{H}}_2 $ -air and $ {\mathrm{H}}_2 $ / $ {\mathrm{NH}}_3 $ -air premixed turbulent flames. This configuration is computed for five different equivalence ratios (Table 1) to assess the ability of the CNN to approximate the reaction rate over a range of equivalence ratios. All other parameters are kept constant. The Reynolds number of the central inlet is about $ \mathrm{10,000} $ for all cases.

Figure 2. Diagram of the slot burner configuration used to generate the DNS database. The flame is depicted by an iso-surface at a progress variable $ c=0.5 $ colored by $ \dot{\omega}=-{\dot{\omega}}_{{\mathrm{H}}_2} $ . The dimensions of the domain are annotated on the right. The length $ {L}_x $ is adapted to the length of the turbulent flame brush, which is a function of the global equivalence ratio $ {\phi}_g $ _.

Table 1. Global equivalence ratio $ {\phi}_g $ for the five different DNS cases

The laminar flame characteristic values are also indicated. $ {S}_L^0 $ is the laminar flame speed, $ {\delta}_{th,L}^0 $ the laminar thermal flame thickness, $ {T}^{ad} $ the adiabatic flame temperature.

The data generation and training strategy is based on processing multiple snapshots from the five DNS (Figure 3):

1. a. The five DNS cases (Table 1) are computed for $ 59 $ ms from a statistically steady state. Three-Dimensional (3D) instantaneous solutions are saved every ms. Solutions with an index multiple of 10 are set aside for testing after training. All other solutions are used for training and validation data. This generates 6 solutions for testing and 54 solutions for training and validation data for each DNS case.

2. b. A spatial filtering is performed on each solution to obtain $ \overline{\rho} $ , $ \overline{\rho {Y}_{{\mathrm{H}}_2}} $ , $ \overline{\rho \xi} $ and $ {\overline{\dot{\omega}}}_{{\mathrm{H}}_2} $ . The filtered variables are downsampled to a coarser grid, representing an LES grid. The filtering and downsampling operation is described in Section 2.3.2.

3. c. The filtered progress variable $ \tilde{c} $ is calculated using Eq. (2.3) with $ {\tilde{Y}}_{{\mathrm{H}}_2}=\overline{\rho {Y}_{{\mathrm{H}}_2}}/\overline{\rho} $ and $ \tilde{\xi}=\overline{\rho \xi}/\overline{\rho} $ . The filtered equivalence ratio $ \tilde{\phi} $ is calculated using Eq. (2.5) with $ \tilde{\xi}=\overline{\rho \xi}/\overline{\rho} $ . The filtered burning rate $ \overline{\dot{\omega}} $ is calculated using Eq. (2.1) with $ {\overline{\dot{\omega}}}_{{\mathrm{H}}_2} $ .

4. d. It is important that the CNN does not learn the specific configuration of the slot burner, but only learns generic flame elements. Depending on the effective receptive field of the CNN, training using the emulated LES solutions directly could encode in the network information on the distance to the inlets, outlet, or the wall for example. This would drastically reduce CNN’s ability to generalize and must be avoided. To address this, training is performed on randomly selected cubes in the domain as in Lapeyre et al. (Reference Lapeyre, Misdariis, Cazard and Poinsot2018, Reference Lapeyre, Misdariis, Cazard, Veynante and Poinsot2019). This also enables the CNN to be invariant to translation (Biscione and Bowers, Reference Biscione and Bowers2021).

5. e. The CNN is trained to approximate $ \overline{\dot{\omega}} $ (Eq. (2.6)) using data from the cubic samples. 10% of samples are randomly selected and set aside for validation, the rest are used for training. The complete description of the training is given in Section 2.3.3.

Figure 3. Diagram of the strategy used to generate data and train the CNN. $ {\delta}_L^0/{\delta}_L^1 $ is the inverse of the laminar flame thickening due to filtering (Section 2.3.2). $ \tilde{\phi} $ is calculated from $ \tilde{\xi} $ using Eq. (2.5).

2.3.1. Direct numerical simulation

DNS of the slot burner cases (Table 1) are performed using the AVBP (Schonfeld and Rudgyard, Reference Schonfeld and Rudgyard1999; Gicquel et al., Reference Gicquel, Gourdain, Boussuge, Deniau, Staffelbach, Wolf and Poinsot2011) massively parallel code solving the compressible multi-species Navier–Stokes equations. A third-order accurate Taylor–Galerkin scheme is adopted for the discretization of the convective terms (Colin and Rudgyard, Reference Colin and Rudgyard2000). Navier–Stokes Characteristic Boundary Conditions (NSCBCs) (Poinsot and Lelef, Reference Poinsot and Lelef1992) are imposed at the inlets (relaxation factor of $ 1000 $ $ {\mathrm{s}}^{-1} $ ) and at the outlet (relaxation factor of $ 200 $ $ {\mathrm{s}}^{-1} $ ). Dynamic viscosity $ \mu $ follows a power law function of temperature $ T $ _.

(2.7) $$ \mu ={\mu}_0{\left(\frac{T}{T_0}\right)}^{\gamma }, $$

with $ {\mu}_0=8.062\times {10}^{-5} $ $ \mathrm{kg}/\mathrm{m}/\mathrm{s} $ , $ {T}_0=2.645\times {10}^3 $ K and $ \gamma =6.481\times {10}^{-1} $ . Thermal diffusivity is computed from the viscosity using a constant Prandtl number: $ \Pr =0.66 $ . Species diffusivities are computed using a constant Schmidt number specific for each species (Table 2). This approach takes into account non-unity Lewis numbers and preferential diffusion between the different species. It is verified that the errors made by the simplified transport description are negligible by comparing the results with simulations using a mixture-averaged transport model, see Supplementary Material. Soret and Dufour transport processes are ignored in the simulations of the present work. Including these transport, processes could enhance thermodiffusive instabilities and increase source terms and flame surface curvature (Schlup and Blanquart, Reference Schlup and Blanquart2018). Although this would require another training on the new database, it does not affect the conclusions of this work. Indeed, it has been shown that neglecting these effects still leads to simulations that capture much of the peculiar behavior of lean-premixed $ {\mathrm{H}}_2 $ -air flames including thermodiffusive instabilities and enhanced burning rates (Aspden et al., Reference Aspden, Day and Bell2015; Aspden, Reference Aspden2017; Schlup and Blanquart, Reference Schlup and Blanquart2018; Pitsch, Reference Pitsch2024). Hydrogen chemical kinetics relies on the San Diego mechanism (Saxena and Williams, Reference Saxena and Williams2006), already successfully used for $ {\mathrm{H}}_2 $ -air premixed combustion in Coulon et al. (Reference Coulon, Gaucherand, Xing, Laera, Lapeyre and Poinsot2023). This mechanism comprises $ 9 $ species and $ 21 $ reactions.

Table 2. Schmidt $ {\mathrm{Sc}}_k $ and Lewis $ {\mathrm{Le}}_k $ numbers for the species used for the DNS of $ {\mathrm{H}}_2 $ -air flames

Table 3. Flame characteristics for the three sets of LES parameters ( $ \sigma $ , DSF) used in this work

$ {\Delta}_x $ is the grid cell size, equal to DSF times the original cell size of $ 0.1 $ mm. $ {\delta}_{c,L}^0 $ and $ {\delta}_{c,L}^1 $ are the laminar flame thickness based on the progress variable of a non-filtered and filtered 1D flame (Eq. (2.9)). $ \mathrm{RI}={\delta}_{c,L}^1/{\Delta}_x+1 $ is the resolution index of the filtered flame on the LES grid.

The mesh is a homogeneous cartesian grid with constant element size $ {\Delta}_x=0.1 $ mm. This gives a resolution index $ \mathrm{RI}={\delta}_{th,L}^0/{\Delta}_x+1 $ that varies between $ 12.5 $ and $ 5.2 $ for $ {\phi}_g $ between $ 0.35 $ and $ 0.7 $ . $ {\delta}_{th,L}^0 $ is the laminar thermal flame thickness. A mesh sensitivity study was carried out where we verified that turbulent combustion statistics are not affected by mesh refinements. This study is included as Supplementary Material. The length of the domain in the $ x $ -direction $ {L}_x $ is adapted to the length of the turbulent flame brush. It varies from $ 76 $ mm for $ {\phi}_g=0.35 $ to $ 36 $ mm for $ {\phi}_g=0.7 $ (Figure 2). The number of grid points accordingly varies from $ 39 $ to $ 19 $ million. Calculations were carried out on the new Alps Research Infrastructure at the Swiss National Supercomputing Centre (CSCS) using AMD EPYC 7742 CPUs. A total of $ 1.6 $ million core hours were used to generate the database. Depending on the grid size, computations are performed using between $ \mathrm{5,120} $ and $ \mathrm{2,560} $ computing cores with good efficiency thanks to the massively parallel computing ability of the AVBP code.

2.3.2. Filtering and downsampling

A spatial filtering is performed to obtain $ \overline{\rho} $ , $ \overline{\rho {Y}_{{\mathrm{H}}_2}} $ , $ \overline{\rho \xi} $ and $ \overline{{\dot{\omega}}_{{\mathrm{H}}_2}} $ , in order to calculate $ \tilde{c} $ , $ \tilde{\phi} $ and $ \overline{\dot{\omega}} $ (Eqs (2.3), (2.5), and ((2.1)). The filtered quantity of a variable $ \varphi $ from the DNS field is defined by Eq. (1.1). The filtered 3D fields are obtained by successively filtering in each direction $ x $ , $ y $ and $ z $ using a 1D Gaussian filter. The 1D filter function is written in discrete form as

(2.8) $$ F(u)=\left\{\begin{array}{ll}{e}^{-\frac{1}{2}{\left(\frac{u}{\sigma}\right)}^2}& \mathrm{if}\;u\in \left[-z..+z\right]\\ {}0& \mathrm{otherwise}\end{array}\right.\hskip0.22em \mathrm{with}\hskip0.22em z=4\sigma, $$

and then normalized by its sum $ {\sum}_{u=-z}^zF(u) $ . The filtered fields are then downsampled to a coarser grid to emulate LES fields. Three different standard deviations $ \sigma $ are used in this work to emulate three different LES resolutions (Table 3). For each $ \sigma $ , a downsampling factor DSF is associated. It is chosen so as to have a LES grid with a resolution index $ \mathrm{RI}={\delta}_{c,L}^1/{\Delta}_x+1 $ of about 6. $ {\delta}_{c,L}^0 $ and $ {\delta}_{c,L}^1 $ are the laminar flame thickness based on the progress variable of a non-filtered and filtered 1D flame

(2.9) $$ {\delta}_{c,L}^0={\left(\max \left|\nabla c\right|\right)}^{-1}\hskip0.1em \mathrm{and}\hskip0.22em {\delta}_{c,L}^1={\left(\max \left|\nabla \tilde{c}\right|\right)}^{-1}. $$

3.1. Training and testing

The training and validation database comprises $ 40\times 54\times 5\times 3=\mathrm{32,400} $ cubic samples ( $ 40 $ random cubes per solution, $ 54 $ solutions per case, $ 5 $ different equivalence ratios, $ 3 $ different sets of LES parameters (filtering and downsampling)). $ \mathrm{3,240} $ ( $ 10 $ %) of these cubes are randomly set aside for validation. Training the CNN model for $ \mathrm{1,500} $ epochs on this database takes $ 47 $ minutes. $ \mathrm{1,500} $ epochs are sufficient because the loss function evaluated on the validation dataset no longer decreases significantly (Figure 4). The set of model parameters giving rise to the lowest loss is selected as the best model (black circle in Figure 4). This best model is used to assess the ability of the CNN to approximate burning rates on the test solutions (full domain, never seen during training). The Normalized Mean Absolute Error (NMAE) is used to estimate the accuracy of the CNN model (Figure 4). It is defined as follows:

(3.1) $$ \mathrm{NMAE}=\frac{1}{\mathrm{mean}\left({\overline{\dot{\omega}}}^{\ast}\right)}\frac{1}{N}\sum \limits_{i=1}^N\left|{\overline{\dot{\omega}}}_i^{\ast }-{\overline{\dot{\omega}}}_i^{\mathrm{NN}}\right|, $$

where $ {\overline{\dot{\omega}}}^{\ast } $ is the filtered burning rate from the DNS dataset and $ {\overline{\dot{\omega}}}^{\mathrm{NN}} $ is the filtered burning rate modeled by the CNN. The mean absolute error is only about $ 5 $ % of the average burning rate for all cases. The very high accuracy of the CNN model is confirmed by the distribution of CNN-modeled burning rate $ {\overline{\dot{\omega}}}^{\mathrm{NN}} $ versus ground-truth filtered burning rate $ {\overline{\dot{\omega}}}^{\ast } $ (Figure 5). The log-normalized Two-Dimensional (2D) histograms show a large majority of points on the $ x=y $ axis (that is perfect agreement). Figure 6 shows a visualization of the result using a planar cut through a test solution for the case $ {\phi}_g=0.4 $ and for the three sets of LES parameters $ \left(\sigma, \mathrm{DSF}\right) $ . One can see the ability of the CNN to retrieve the extremely complex turbulent flame topology for all filter sizes. Even the burning pockets detached from the main flame brush are correctly modeled by the CNN. This is true for all other equivalence ratios. The visualization of the cases $ {\phi}_g=0.35 $ , $ 0.5 $ , $ 0.6 $ and $ 0.7 $ are given as Supplementary Material to avoid overloading the manuscript.

Figure 4. Left: Evolution of the RMSE during training, evaluated over the training dataset (red circles) and the validation dataset (blue squares). Black solid lines are moving averages of the RMSE. The black circle shows the lowest RMSE over the validation dataset. The model parameters at this specific epoch are selected. Right: Normalized mean absolute error over the testing solutions (Eq. (3.1)) for the different equivalence ratios and LES parameters used for building the training dataset. Error bars show the first and third quartiles of the data points.

Figure 5. Scatter plots with 2D histograms: CNN-modeled burning rate $ {\overline{\dot{\omega}}}^{\mathrm{NN}} $ versus ground-truth filtered burning rate $ {\overline{\dot{\omega}}}^{\ast } $ . Individual values are normalized by the maximum burning rate in the datasets. The points used for the histograms have a progress variable $ c $ : $ 0.05\le c\le 0.95 $ . Histogram values below the color scale are transparent. Gray dashed line indicates $ x=y $ (that is zero error). Each column corresponds to a global equivalence ratio. Each row corresponds to a set of LES parameters (filtering and downsampling). Data are collected from the testing solutions.

Figure 6. Planar cut normal to the $ z $ -axis, in the middle of the domain, colored by the ground-truth filtered burning rate $ {\overline{\dot{\omega}}}^{\ast } $ and the CNN-modeled burning rate $ {\overline{\dot{\omega}}}^{\mathrm{NN}} $ for three sets of LES parameters. The global equivalence ratio is $ {\phi}_g=0.4 $ . The complex turbulent flame topology of the burning rates is remarkably well reproduced by the CNN.

3.2. Generalization ability

The ability of the CNN model to generalize outside the cases for which it has been trained is now assessed. For this, two new sets of LES parameters $ \left(\sigma =6,\mathrm{DSF}=3\right) $ and $ \left(\sigma =12,\mathrm{DSF}=6\right) $ , and two new equivalence ratios $ {\phi}_g=0.45 $ and $ 0.55 $ are used.

Filtering and downsampling with the two new sets of LES parameters are performed on test solutions of the existing DNS with global equivalence ratios $ {\phi}_g=0.35 $ , $ 0.4 $ , $ 0.5 $ , $ 0.6 $ and $ 0.7 $ . Although the distribution of points $ {\overline{\dot{\omega}}}^{\mathrm{NN}} $ versus $ {\overline{\dot{\omega}}}^{\ast } $ is a little more scattered (Figure 7), the majority of the CNN-modeled burning rates are in very good agreement with their ground-truth counterparts. This demonstrates the ability of the model to generalize between two LES filter sizes for which it has been trained. It is important to note, however, that training with all three filter sizes (Table 3) is necessary to achieve this result. We have tried to train the CNN with only two filter sizes: $ \sigma =4 $ and $ 16 $ , and found that the resulting model has a significant bias for $ \sigma =8 $ , see Supplementary Material.

Figure 7. Scatter plots with 2D histograms: CNN-modeled burning rate $ {\overline{\dot{\omega}}}^{\mathrm{NN}} $ versus ground-truth filtered burning rate $ {\overline{\dot{\omega}}}^{\ast } $ . Individual values are normalized by the maximum burning rate in the datasets. The points used for the histograms have a progress variable $ c $ : $ 0.05\le c\le 0.95 $ . Histogram values below the color scale are transparent. Gray dashed line indicates $ x=y $ (that is zero error). Each column corresponds to a global equivalence ratio. Each row corresponds to a set of LES parameters (filtering and downsampling). Data are collected from solutions with two sets of LES parameters that were not included in the training dataset. The CNN approximates the burning rates with good accuracy, demonstrating the ability to generalize to other LES parameters from which it has been trained.

The ability of the model to approximate burning rates for other global equivalence ratios is now assessed. For this, two new DNSs are performed with a global equivalence ratio $ {\phi}_g=0.45 $ and $ 0.55 $ and six solutions per DNS (snapshot every $ 1 $ ms) are used to evaluate the model on the filter sizes $ \sigma =4 $ , $ 8 $ and $ 16 $ (Figure 8). The distribution of points $ {\overline{\dot{\omega}}}^{\mathrm{NN}} $ versus $ {\overline{\dot{\omega}}}^{\ast } $ is a little more scattered compared to inference on solutions with a global equivalence ratio used for training. The majority of the points are nevertheless clustered along the $ x=y $ axis: the approximation of the burning rates is still very good. The planar cuts show a flame topology very well retrieved by the CNN model. These excellent results are due to the large number of equivalence ratios included in the training dataset. We have trained the CNN with only $ {\phi}_g=0.35 $ , $ 0.4 $ , $ 0.6 $ and $ 0.7 $ (that is without $ {\phi}_g=0.5 $ ) and that model is flawed for the global equivalence ratio $ {\phi}_g=0.5 $ , see Supplementary Material.

Figure 8. Scatter plots and 2D histograms together with planar cuts comparing CNN-modeled burning rate $ {\overline{\dot{\omega}}}^{\mathrm{NN}} $ to ground-truth filtered burning rate $ {\overline{\dot{\omega}}}^{\ast } $ . See caption of Figure 5 for a full description of how the histograms are constructed. The two global equivalence ratios $ {\phi}_g=0.45 $ and $ 0.55 $ were not included in the training dataset. The CNN approximates the burning rates with high accuracy, demonstrating the ability to generalize to other equivalence ratios from which it has been trained.

These generalization results are very promising and show that a CNN model can generalize with a reasonable amount of training cases. This is particularly important with regard to the global equivalence ratio $ {\phi}_g $ . Computing a DNS for each global equivalence ratio has a significant computing cost. In contrast, the computational cost of a filtering and downsampling operation to produce a new set of LES parameters $ \left(\sigma, \mathrm{DSF}\right) $ is negligible.

3.3. Comparison with a filtered tabulated chemistry approach

A filtered tabulated chemistry modeling is implemented to compare CNN results with an established approach for modeling filtered burning rates. In the filtered tabulated chemistry framework, the flame structure in the direction normal to the flame front is assumed identical to the structure of a planar 1D freely propagating premixed laminar flame (Fiorina et al., Reference Fiorina, Vicquelin, Auzillon, Darabiha, Gicquel and Veynante2010). For a given filter size, it is then possible to apply the filter operation introduced in Section 2.3.2 (Eqs (1.1) and (2.8)) to a precomputed 1D flame at a given global equivalence ratio $ {\phi}_g $ and extract the burning rates directly as

(3.2) $$ \overline{\dot{\omega}}\approx {\overline{\dot{\omega}}}^{+}\left(\tilde{c},{\phi}_g\right), $$

where $ {\overline{\dot{\omega}}}^{+} $ is the filtered burning rate of the 1D flame, mapped as a function of the filtered progress variable $ \tilde{c} $ for a given filter size $ \sigma $ . To account for the local fluctuation of the equivalence ratio in the domain, it is possible to retrieve the reaction rate in a 1D flame computed with the local equivalence ratio $ \phi $ used to define the fresh gas composition. In this case, we write

(3.3) $$ \overline{\dot{\omega}}\approx {\overline{\dot{\omega}}}^{+}\left(\tilde{c},\tilde{\phi}\right). $$

Flame wrinkling at the subfilter-scale level is lost in LES. A subfilter-scale wrinkling factor $ \Xi $ needs to be applied to compensate for the loss in flame surface (Poinsot and Veynante, Reference Poinsot and Veynante2011). The filtered burning rate for filtered tabulated chemistry modeling therefore reads

(3.4) $$ {\overline{\dot{\omega}}}^F=\Xi \hskip0.1em {\overline{\dot{\omega}}}^{+}\left(\tilde{c},{\phi}_g\right)\;\mathrm{or}\;{\overline{\dot{\omega}}}^{FC}=\Xi \hskip0.1em {\overline{\dot{\omega}}}^{+}\left(\tilde{c},\tilde{\phi}\right), $$

depending on whether the global or local equivalence ratio is used to define the fresh gas composition in the corresponding 1D flame. The wrinkling factor $ \Xi $ is usually modeled by algebraic expressions. Models based on a fractal description of the flame front (Gouldin, Reference Gouldin1987; Gouldin et al., Reference Gouldin, Bray and Chen1989) are very common. They assume that in a range of physical scales bounded by an inner cutoff $ \eta $ and an outer cutoff $ L $ , the flame front is a fractal surface of dimension $ 2\le {D}_f\le 3 $ . The wrinkling factor is given by

(3.5) $$ \Xi ={\left(\frac{L}{\eta}\right)}^{D_f-2}. $$

$ L $ corresponds to the size of the largest unresolved wrinkles, which is of the order of the thickness of the LES flame estimated as $ {\delta}_{c,L}^1 $ (Eq. (2.9)). $ \eta $ is the size of the smallest wrinkles and varies from $ \eta =L $ (that is no unresolved wrinkling, $ \Xi =1 $ ) to $ \eta ={\delta}_{c,L}^0 $ (that is filling of the space by a flame surface of fractal dimension $ {D}_f $ between cut-off scales $ {\delta}_{c,L}^0 $ and $ L $ , $ \Xi ={\left(L/{\delta}_{c,L}^0\right)}^{D_f-2} $ ). The latter situation is reached when turbulence intensity is sufficiently high, which is often the case in practice (Veynante et al., Reference Veynante, Schmitt, Boileau and Moureau2012; Veynante and Moureau, Reference Veynante and Moureau2015). We then assume $ \eta ={\delta}_{c,L}^0 $ in Eq. (3.5), maximizing the flame surface density at the subfilter-scale level. The fractal dimension is estimated as $ {D}_f=2.5 $ following the work of Charlette et al. (Reference Charlette, Meneveau and Veynante2002). The set of 1D flames is computed using the San Diego mechanism (Saxena and Williams, Reference Saxena and Williams2006) for the chemical kinetics (same as the training database for the CNN) and a mixture-averaged transport model for the diffusivities.

CNN-based modeling significantly outperforms the filtered tabulated chemistry one, especially for the lean cases (e.g. $ {\phi}_g=0.35 $ in Figure 9). Due to the extremely high thermodiffusive effects in very lean cases, the burning rates of the turbulent 3D field can be much higher than those in a freely propagating planar flame at the same equivalence ratio. For example, for the case $ {\phi}_g=0.35 $ , $ \sigma =4 $ , $ \mathrm{DSF}=2 $ , the maximum burning rate in the 3D field $ \max \left({\overline{\dot{\omega}}}^{\ast}\right) $ is about $ 17 $ times greater than that in a laminar planar flame $ \max \left({\overline{\dot{\omega}}}^{+}\right) $ . This is due to the fact that flame strain is not considered in the tabulation method applied and that the subfilter wrinkling factor $ \Xi $ only takes into account turbulence-induced wrinkling, ignoring that induced by thermodiffusive instabilities and their interactions with turbulence (Berger et al., Reference Berger, Attili and Pitsch2022c). This $ {\phi}_g=0.35 $ case is a truly extreme test case to illustrate the strength of the CNN model. Tabulated chemistry modeling does not fail so spectacularly for equivalence ratios close to stoichiometry, where preferential and differential diffusion effects are less pronounced (Berger et al., Reference Berger, Attili and Pitsch2022b) (e.g. $ {\phi}_g=0.7 $ in Figure 9). The CNN results are still much more accurate thanks to the consideration of the 3D flame structure via successive convolutions. To further demonstrate the CNN’s ability to learn from flame topology, training has been performed without information on the equivalence ratio for the case $ {\phi}_g=0.35 $ . This study is given as Supplementary Material. It is found that the burning rate approximation is still very accurate, proving our initial claim that a CNN is capable of recognizing complex flame topology and inferring associated burning rates by learning spatial hierarchies of features.

Figure 9. Scatter plots of the CNN-modeled burning rate $ {\overline{\dot{\omega}}}^{\mathrm{NN}} $ and the filtered tabulated chemistry burning rate $ {\overline{\dot{\omega}}}^{\mathrm{F}} $ or $ {\overline{\dot{\omega}}}^{\mathrm{FC}} $ (Eq. (3.4)) versus ground-truth filtered burning rate $ {\overline{\dot{\omega}}}^{\ast } $ . Individual values are normalized by the maximum burning rate in the datasets. Gray dashed line indicates $ x=y $ (that is zero error). Two cases of LES parameter sets are presented for two global equivalence ratios: the lowest equivalence ratio $ {\phi}_g=0.35 $ for which the filtered tabulated chemistry is significantly flawed due to very strong thermodiffusive effects; the highest equivalence ratio $ {\phi}_g=0.7 $ for which thermodiffusive effects are less important.