# A Data-Driven Approach to Refine the Partially Stirred Reactor Closure for Turbulent Premixed Flames

Lorenzo Piu, Arthur Péquin, Rodolfo S M Freitas, Salvatore Iavarone, Heinz Pitsch, Alessandro Parente

► To cite this version:

Lorenzo Piu, Arthur Péquin, Rodolfo S M Freitas, Salvatore Iavarone, Heinz Pitsch, et al., A Data-Driven Approach to Refine the Partially Stirred Reactor Closure for Turbulent Premixed Flames. Flow, Turbulence and Combustion, In press, 10.1007/s10494-024-00626-3, hal-05225694

HAL Id: hal-05225694

https://hal.science/hal-05225694v1

Submitted on 27 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Data-Driven Approach to Refine the Partially Stirred Reactor Closure for Turbulent Premixed Flames

Lorenzo Piu¹,² · Arthur Péquin¹,² · Rodolfo S. M. Freitas³ · Salvatore Iavarone¹,²,⁴ · Heinz Pitsch⁵ · Alessandro Parente¹,²,⁶

Received: 11 September 2024 / Accepted: 4 December 2024
© The Author(s) 2025

## Abstract

Accurately predicting turbulent combustion processes is fundamental for optimizing efficiency, reducing pollutant emissions, and ensuring operational safety in combustion systems. To this purpose, computational fluid dynamics (CFD) simulations are widely employed. In particular, large eddy simulations (LES) balance prediction accuracy with computational efficiency by resolving only the most energy-containing scales of turbulence and rely on modeling the turbulence-chemistry interactions (TCI) occurring at the smallest scales. Among the existing closures, the partially stirred reactor (PaSR) model incorporates finite-rate chemistry and estimates a cell reacting fraction based on the local Damköhler number to account for the subfilter-scale TCI. Although widely validated in CFD computations, the PaSR model was found limited by the way it computes the cell reacting fraction. To tackle this point, our study proposes a machine learning (ML) enhanced partially stirred reactor model for LES. A fully connected neural network is trained on direct numerical simulation (DNS) data of turbulent premixed jet flames to compute a correction coefficient for the cell reacting fraction. Maintaining the original model shape, this ML-enhanced closure aims at bridging the gap between physics-based models and advanced data-driven techniques. The proposed formulation not only improves the prediction accuracy of quantities of interest such as the heat release rate but also features computational feasibility and generalisation capabilities over a large range of LES grid refinement. This demonstrates the significant potential of ML-aided TCI closures in future applications of combustion engineering.

**Keywords** Machine learning · Neural networks · Turbulent combustion modeling · Direct numerical simulation (DNS) · Large eddy simulations (LES) · Partially stirred reactor (PaSR)

## 1 Introduction

A correct description of turbulent combustion processes is essential for optimizing their performance in terms of efficiency, emissions, and operational safety. Due to its ability to combine turbulence and chemistry modeling, computational fluid dynamics (CFD) is a

<footnote>Extended author information available on the last page of the article</footnote>
key instrument for simulating turbulent reacting flows. Reynolds Averaged Navier-Stokes (RANS) equations have been traditionally used for turbulent combustion modeling Brankovic et al. (2000); Poinsot and Veynante (2005); Amaduzzi et al. (2021); Huang et al. (2022) due to the relatively low computational cost and ability to provide time-averaged solutions. In the RANS approach, the governing equations of fluid motion are ensemble averaged, leading to a set of equations for the mean quantities. Turbulent fluctuations are modeled using various approaches, such as the $κ-ε, κ-ω$, and Reynolds stress models (RSM) Wilcox (2006); Pope (2004). These models are generally robust and computationally efficient, making RANS suitable for engineering applications where quick turnaround times are necessary Veynante and Vervisch (2002). However, RANS models struggle to accurately capture the complex and unstable nature of turbulent combustion, especially in flows with strong interactions between turbulence and chemistry. Moin and Apte (2006). In contrast, large eddy simulations (LES) provide a more detailed representation of turbulent flows by directly resolving the large, energy-containing eddies while modeling only the smaller, subfilter-scale (SFS) motions Sagaut (2001); Garnier et al. (2009). This approach significantly improves the accuracy of simulations for unsteady and complex flow phenomena, which are common in turbulent combustion. Peters (2000). LES can capture transient features of the flow Pitsch (2006), such as vortex shedding, flame instabilities, and local extinction/ignition events, which are averaged out in RANS simulations. Meneveau and Katz (2000). The primary challenge with LES has historically been its high computational cost, as resolving the large eddies requires a finer spatial and temporal resolution compared to RANS. Sagaut (2001). The exponential growth in processing capabilities, coupled with the development of high-performance computing (HPC) architectures, has made LES more accessible for practical engineering applications. Pérez Arroyo et al. (2020); Fureby et al. (2022). Moreover, the advent of parallel computing techniques and more efficient solvers have further reduced the computational burden associated with LES. Choi and Moin (2012). Consequently, LES is gaining momentum as a method for turbulent combustion modeling, offering an acceptable trade-off between accuracy and computational feasibility.

Both simulation approaches require closure for the unresolved terms. Reactor-based models Péquin et al. (2023) play a crucial role in turbulent combustion modeling, allowing the treatment of finite-rate chemistry at an affordable computational cost. Two primary reactor-based models used in turbulent combustion are the eddy dissipation concept (EDC) and the partially stirred reactor (PaSR) model. The EDC, originally proposed by Magnussen and Hjertager (1977); Magnussen (1989), is an extension of the eddy break-up (EBU) model, which describes the interaction between turbulence and combustion based on the assumption that the reaction rate is controlled by the rate at which turbulent eddies mix reactants at the molecular level. In the EDC model, small turbulent structures, or fine scales, are treated as a series of perfectly stirred reactors where chemical reactions occur. Gran and Magnussen (1996). This approach allows the EDC model to account for detailed chemical kinetics within the turbulent flow, providing a more accurate representation of the combustion process compared to simpler models like EBU Magnussen and Hjertager (1977). The EDC model is particularly useful in flows at high Reynolds numbers, where the interaction between turbulence and chemistry is complex and highly nonlinear. The partially stirred reactor model, first introduced in the works of Chomiak and Karlsson (1996); Golovitchev et al. (1998, 2000); Golovitchev (2001), and Karlsson (1995), naturally generalizes the EDC, providing robustness and independence on the combustion regime. Péquin et al. (2023). The PaSR model also relies on the assumption that chemical reactions take place inside the fine structures, with characteristic length scales smaller or equal than the computational cell size. Each cell volume is divided into two zones: a well-mixed reactive
zone modeled as an ideal reactor, and a non-reactive zone where mixing with the surround-
ing fluid takes place. In LES, the PaSR model is particularly advantageous because it can
explicitly account for subfilter scale mixing, modeling the cell reacting fraction from the
resolved quantities.

Turbulent combustion closures can present intrinsic deficiencies due to modeling approximations. A method to evaluate existing and novel SFS closures is through a priori validation Trisjono and Pitsch (2015). Such an approach compares the modeled quantities with those filtered from direct numerical simulations (DNS), which resolve the whole turbulence scale spectrum and do not require any turbulence-chemistry interactions modeling. Many studies used this technique to assess the PaSR closure. Iavarone et al. (2021) performed an analysis on a three-dimensional DNS to evaluate the PaSR model in moderate or intense low-oxygen dilution (MILD) combustion regime. Péquin et al. (2022) estimated a priori the performance of PaSR in LES, highlighting model shortcomings due to the definition of the cell reacting fraction. Quadarella et al. (2023) presented and validated a priori a generalized PaSR model, which involves directly manipulating the species production rates by employing a modal decomposition of the Jacobian matrix to account for the entire range of chemical timescales present in a reactive system.

Machine learning (ML) is a useful tool for improving the predictive capabilities of SFS closures in fluid dynamics. The ML algorithms' ability to handle large datasets is well suited to extract information regarding the correlations between turbulence and chemistry from DNS simulations. Swaminathan and Parente (2023) have extensively discussed the role of ML in the development of sustainable combustion technologies. A review of the applications of ML in reactive flows was recently written by Ihme et al. (2022). Bode et al. (2019, 2021) implemented physics-informed enhanced super-resolution generative adversarial networks (PIESRGANs) to model subgrid quantities in combustion LES. They used PIESRGANs to accurately reconstruct fully resolved turbulence fields from filtered data, combining adversarial and physics-informed losses to improve generalization capabilities. The network was tested in complex reactive and inert LES scenarios, demonstrating its potential for accurate subfilter modeling in high Reynolds number flows. Nista et al. (2023) further investigated the generalizability of generative adversarial networks (GANs) for predicting the scalar-flux and unresolved subfilter-stress tensors. The study utilized a super-resolution GAN to model these quantities in LES of turbulent premixed reacting flows, demonstrating the network's ability to generalize across different Reynolds and Karlovitz numbers. The results indicated that the GAN could accurately reconstruct high-resolution fields from low-resolution inputs, outperforming traditional models under various conditions and even in scenarios not represented in the training dataset. Henry de Frahan et al. (2019) evaluated various machine learning algorithms, i.e., a random forest, a deep neural network, and a variational autoencoder, for representing joint PDFs of mixture fraction and progress variable, considering a single snapshot of a DNS from a lean premixed low-swirl burner. In a study by Chen et al. (2021), deep neural networks were shown to outperform algebraic models at large filter sizes when modeling the joint filtered density functions (FDFs) of mixture fraction and progress variable from DNS of MILD combustion. The generalization capabilities of the networks trained on MILD combustion data were subsequently tested by Yang et al. (2023) on five DNS of premixed flames: three H₂-air (planar, V-shaped, swirled) flames and two CH₄-air (twin V-, stagnation) flames. The progress variable source terms were predicted fairly well, although some discrepancies were observed for the H₂-air cases as the networks were trained on reacting CH₄-air-exhaust gases mixtures Doan et al. (2018). Lapeyre et al. (2019) used a convolutional neural network (CNN) to regress the subgrid-scale flame surface density in premixed turbulent
flames, with training data generated from filtered DNS results sampled onto a coarse mesh to train a U-Net, a popular CNN architecture. Additionally, CNNs were explored in a priori studies for deconvolution tasks and approximating the progress variable variance of filtered DNS data of a freely propagating turbulent premixed flame Nikolaou et al. (2019). Péquin et al. (2024) used supervised clustering to select the optimal formulation for the chemical timescale. The methodology's robustness was demonstrated across different combustion regimes, showcasing its potential to enhance the accuracy of PaSR combustion closures in LES. Such work underlines the importance that the submodel selection for computing the cell reacting fraction has in the PaSR model. Recently, Freitas et al. (2023) used ML to demonstrate potential improvements to the PaSR closure. By optimizing the mathematical definition of the cell reacting fraction using sparse regression techniques, they enhanced the model's accuracy and applicability in turbulent combustion simulations. Finally, Jigjid et al. (2023) proposed a PaSR approach where the cell reacting fraction was predicted by two neural networks to handle both local premixed and non-premixed conditions in MILD combustion. The model was validated a priori using two MILD combustion DNS and the assessment indicated the model's effectiveness in capturing subfilter-scale contributions accurately.

This research seeks to leverage traditional turbulent combustion closures by introduc-
ing a correction term that can compensate for the deficiencies associated with the baseline
PaSR model. The correction term is built using a fully connected neural network (FCNN)
to refine the computation of the cell reacting fraction. The aim is to improve predictive
capabilities and to demonstrate the potential of machine learning to enhance existing model
closures in LES within the limits of a physics-based model. The formulation we propose
still relies on prior knowledge of the analytical form of the cell reacting fraction, while
only employing the FCNN to fine-tune inaccuracies due to model assumptions.

The present paper is organized as follows. Section 2 describes the DNS datasets used to train and test the neural networks. In Sect. 3, we describe the overall methodology employed in this study. In Sect. 4, we show and analyze the results obtained. A discussion is then conducted regarding the main deficiencies of the functional form of the cell reacting fraction in the original PaSR model, the generalizability of the FCNN across filter sizes and DNS datasets, and the importance of training parsimonious ML algorithms. Finally, the main conclusions of this research are summarized in Sect. 5.

## 2 DNS Datasets Description

The availability of high-fidelity DNS data is of fundamental importance to obtain infor-
mation regarding the subfilter-scale interactions between turbulence and chemistry. In the
present work, two different DNS cases are considered for the training and the validation
of the ML model. ‘Case 1’ refers to a 3D DNS simulation of a slot turbulent premixed
methane/air jet flame Luca et al. (2019); Attili et al. (2021) of moderate Reynolds number
Re = 5600. Further details on this dataset are presented in Sect. 2.1. To assess the mod-
el’s performance on out-of-sample data, the trained neural network is tested on a differ-
ent methane/air premixed jet flame Brouzet et al. (2021) employing a skeletal COFFEE
mechanism for methane Coffee (1984). This second DNS case, referred as ‘Case 2’ in the
following, is detailed in Sect. 2.2. Table 1 summarizes the characteristics of the two DNS
cases.
## 2.1 Case 1: Slot Jet Flame

The central jet includes a mixture of methane and air, with an equivalence ratio $\phi = 0.7$ at a temperature of 800 K, and is enclosed by a coflow of burnt gases Luca et al. (2019); Attili et al. (2021). The background pressure in the computational domain is set to 4 atm. The laminar flame speed and thermal thickness were estimated to be $S_L = 1 \ m \cdot s^{-1}$ and $\delta_L = 110 \ \mu m$ based on one-dimensional simulations of a freely propagating flame. From these quantities, the flame characteristic time scale can be defined as $\tau_L = \delta_L/S_L \approx 110 \ \mu s$. The domain is discretized with a uniform mesh comprising 1440 x 962 x 256 points along the x-, y-, and z-axis, which are the streamwise, crosswise, and spanwise directions, respectively. The lateral boundaries perpendicular to the y-axis use no-slip conditions while along the z-axis a periodic boundary condition is employed. The outlet uses an open boundary condition. In the simulation, a skeletal mechanism for methane Luca et al. (2018) of 16 species and 72 reactions was used. In the present work, a subset was used to perform the analysis, comprising $\approx 20$ million points out of the $\approx 354$ million in the original dataset. Both the original and reduced dataset are shown in Fig. 1.

## 2.2 Case 2: Round Jet Flame

The flame consists of a turbulent premixed methane/air round jet flame in an open environment Brouzet et al. (2021). The central jet has a stoichiometric equivalence ratio, i.e., $\phi = 1$, and a temperature of 700 K, at Reynolds number $Re = 5300$. The central jet is surrounded by a coflow to stabilize the flame. The laminar flame speed was estimated $S_L = 1 \ m \cdot s^{-1}$, while the thermal thickness is $\delta_L = 292 \ \mu m$. It follows that the flame characteristic timescale is

<figure><img src="image_2.png" /><figcaption>**Fig. 1** DNS dataset visualization of Case 1 (adapted from Freitas et al. (2023)). **a** Full domain of the premixed methane/air jet flame Luca et al. (2019); Attili et al. (2021). **b** Subset used in this study, showing a 2D temperature colormap and a 3D isosurface of the progress variable $C_{O_2} = 0.73$</figcaption></figure>
$\tau_L$ = 158 µs. The Karlovitz number is $Ka$ = 15.43. The turbulent flame can be classified in the 'thin reaction zone' regime.

The computational domain comprises 1546 grid points in the x direction (streamwise), and 676 points in the transverse directions (y and z). This accounts for a total of approximately 706 million grid points. The chemical kinetics is solved using a COFFEE skeletal mechanism, comprising 14 species and 38 reactions Coffee (1984). A subset comprising 300×300×300 grid points was extracted to perform the analysis. The full dataset and the subset are represented in Fig. 2.

# 3 Methodology

## 3.1 Extraction of ‘LES-like’ Field from DNS

The data processing aims to obtain ‘LES-like’ fields from the DNS snapshot. To accomplish this, a filtering procedure is applied to the original DNS instantaneous field. A filtered quantity $\bar{\phi}$ is obtained using a Gaussian filtering operation. The filtering function for a general variable $\phi$ is defined as:

$$ \overline{\phi(\mathbf{x}, t_0)} = \int_{\Omega} G(\mathbf{r})\phi(\mathbf{x} - \mathbf{r})d\mathbf{r}, \qquad (1) $$

where $\Omega$ is the 3D domain of interest and $G(\mathbf{r})$ is a homogeneous Gaussian filter function defined as

$$ G(\mathbf{r}) = \frac{1}{\sigma\sqrt{2\pi}}\exp\left(\frac{\mathbf{r}^2}{2\sigma^2}\right). \qquad (2) $$

According to Pope (2004), the second moment of the Gaussian $\sigma^2$ is defined in such a way that the Gaussian function maintains the same second moment as a box filter function with amplitude $\delta$:

<figure><img src="image_3.png" /><figcaption>Fig. 2 Graphical depiction of DNS data from Brouzet et al. (2021): (a) 3D contour plot showing the mass fraction of OH, (b) Temperature distribution of the flame in the z midplane, and (c) Subset utilized for testing the FCNN</figcaption></figure>
$$ \sigma^2 = \frac{1}{12} \delta^2, \qquad (3) $$

where $\delta$ is the filter amplitude and

$$ \Delta = \frac{\delta}{\delta^{\text{DNS}}} \qquad (4) $$

is the ratio between the filter amplitude and the grid spacing in the original DNS simulation. Once the filtered density has been obtained using this procedure, the thermochemical variables are determined using Favre filtering, denoted by $\tilde{\cdot}$:

$$ \tilde{\phi} = \frac{\overline{\rho \phi}}{\overline{\rho}}. \qquad (5) $$

From the filtered 'LES-like' fields, the source terms and heat release rate can be computed using a closure model.

## 3.2 The Partially Stirred Reactor (PaSR) Model

The formulation of the PaSR proposed by Chomiak and Karlsson (1996) divides the cell into two parts: fine structures and surroundings. Chemical reactions are supposed to happen only in fine structures. In the original model, the source terms are computed by integrating a perfectly stirred reactor over a residence time $\tau^*$. This leads to the following closure equation for the source terms:

$$ \overline{\dot{\omega}_k^{\text{PaSR}}} = \gamma_{\text{PaSR}} \cdot \frac{\overline{\rho(Y_k^* - Y_k^0)}}{\tau^*}, \qquad (6) $$

where $\gamma_{\text{PaSR}}$ is the cell reacting fraction, computed as

$$ \gamma_{\text{PaSR}} = \frac{\tau_c}{\tau_c + \tau_m}, \qquad (7) $$

with $\tau_c$ and $\tau_m$ representing the chemical and mixing timescales, respectively. $Y_k^0$ and $Y_k^*$ are the mass fractions of the k-th species in the fine structures before and after solving the reactor step, respectively. The initial mass fraction in the fine structures $Y_k^0$ is approximated with the average mass fraction inside the computational cell $\tilde{Y}_k$. $\tau^*$ is the residence time inside the reactor, considered equal to the mixing timescale $\tau_m$ in the original model Sabelnikov and Fureby (2013); Nordin (2001). In this work, we consider another formulation for the PaSR model:

$$ \overline{\dot{\omega}_k^{\text{PaSR}}} = \gamma_{\text{PaSR}} \cdot \overline{\dot{\omega}_k^{\text{LFR}}}, \qquad (8) $$

where $\overline{\dot{\omega}_k^{\text{LFR}}}$ are the source terms from the Laminar finite rate (LFR) model Péquin et al. (2022), directly computed from the filtered fields via the Arrhenius equation in every cell of the domain:

$$ \overline{\dot{\omega}_k^{\text{LFR}}} = \dot{\omega}_k(\bar{\rho}, \tilde{T}, \tilde{Y}_1, \tilde{Y}_2, \dots, \tilde{Y}_{N_{sp}}). \qquad (9) $$
Equation (8) will serve as a bridge between the LFR approach and the data-driven closure that is the main focus of this research. The time integration of ideal reactors is here replaced by the instantaneous Arrhenius rates computed from the filtered thermochemical space, significantly reducing the associated computational effort. The considered assumption neglects crucial aspects of the interaction between mixing processes and chemical reactions, which are accounted for in Eq. (6) by using a perfectly stirred reactor (PSR) or a plug flow reactor (PFR). Such a model approximation leaves an additional responsibility to the parameter $\gamma_{PaSR}$, which must also compensate for the fact that the mixing of reactants in the fine structures is not considered by $\bar{\omega}_k^{LFR}$. A strategy for the closure of Eq. (8) that leverages machine learning to refine the computation of $\gamma_{PaSR}$ is proposed. In Appendix 6, we show that selecting the LFR-based formulation of Eq. (8) as the baseline PaSR model enhances the predictive capabilities of the proposed machine learning strategy, which is described in the following section.

## 3.3 ML-Enhanced PaSR

The model we propose aims to overcome the limitations of the existing physics-based formulations presented in the previous section. Equation (8) can be enhanced with an additive term that encompasses the model deficiencies inherent in the cell reacting fraction formulation:

$$ \bar{\omega}_k^{\text{FCNN}} = (\gamma_{\text{PaSR}} + \gamma_{\text{FCNN}}) \cdot \bar{\omega}_k^{\text{LFR}}, \quad (10) $$

where $\gamma_{FCNN}$ is computed using a fully connected neural network (FCNN) trained on the DNS dataset presented in Sect. 2.1. The choice of using the laminar rates appears now more straightforward in Eq. (10): the model deficiencies are addressed by the correction term, which justifies the choice of not integrating ideal reactors in time. The term $\gamma_{FCNN}$ was formulated as an additive correction rather than a multiplicative one, to ensure compatibility with potential sparsity-promoting techniques applied to this term. The significance of defining sparse additive machine learning correction terms for CFD applications is further discussed in Sect. 4.3.

The selection of the neural network’s architecture type was based on the following considerations:

1. Integrating neural networks into LES solvers entails processing millions of cells, causing even slight computational times to become substantial when a forward pass in the NN is applied to each cell. Therefore, the computational cost of the neural network is crucial, favouring simpler architectures without convolutional or recursive layers.

2. Implementing a posteriori more complex architectures, such as Convolutional Neural Networks (CNNs), can be challenging, particularly with unstructured grids.

The choice of an FCNN overall appears more straightforward to support future applications, allowing this approach to be implemented in CFD. The architecture consists of 6 hidden layers, each comprising 64 neurons. A sigmoid activation function is selected for the first layer after the input one, while the other layers implement a rectified linear unit (ReLU). This architecture was optimized via grid search Hastie et al. (2009) with the number of layers ranging between 2 and 20 and the number of neurons per layer varying
between 16 and 256. A schematic representation of the methodology is represented in Fig. 3.

Determining the optimal set of input variables for the NN remains an open research question, which is highly dependent on the ultimate goal of the method. The variables selected in this study were chosen to minimize correlations among the inputs, to ensure alignment with the baseline PaSR model, and in such a way that they can be computed consistently over a wide range of combustion regimes. The chemical and mixing time-scales were chosen because they describe the cell reacting fraction in the original PaSR formulation, see Eq. (7). The progress variable was chosen since it can provide knowledge regarding the chemical reactions' advancement state in premixed flames. The variables used as input are defined as follows:

* the filtered progress variable $\bar{c}$ based on molecular oxygen ($O_2$) mass fraction. This formulation is consistent with the one proposed in Luca et al. (2018), and is computed from the filtered field:

$$ \bar{c} = 1 - \frac{\tilde{Y}_{O_2} - \tilde{Y}_{O_2,b}}{\tilde{Y}_{O_2,u} - \tilde{Y}_{O_2,b}}, \qquad (11) $$

with the subscripts b and u referring to burned and unburned conditions, respectively. Additionally, two alternative formulations of the progress variable, one based on temperature and the other on the mass fractions of the main combustion products were considered. Details on the model's sensitivity to the choice of the progress variable can be found in Appendix 7.

* the average chemical timescale computed from the filtered DNS field with the Chomiak formulation Chomiak and Karlsson (1996):

$$ \frac{1}{\tau_c} = \max\left(\frac{-\bar{\omega}_F}{\tilde{Y}_F}, \frac{-\bar{\omega}_O}{\tilde{Y}_O}\right) / \bar{\rho}, \qquad (12) $$

where $\bar{\omega}_F$ and $\tilde{Y}_F$ are the fuel reaction rate and mass fraction, while $\bar{\omega}_O$ and $\tilde{Y}_O$ refer to the oxidizer;

* the mixing timescale computed with the scalar-based model proposed by Raman and Pitsch (2007),

$$ \tau_m = \frac{\bar{c}''}{\bar{\chi}_c} = \frac{C_z |\nabla \bar{c}|^2}{2(D + D_t) |\nabla \bar{c}|^2}. \qquad (13) $$

In Eq. (13), $\bar{c}''$ and $\bar{\chi}_c$ are the variance and the dissipation rate of the filtered progress variable, respectively. $C_z$ is a model constant, $D$ is the mixture diffusivity, and $D_t$ is the mixture turbulent diffusivity.

Some of the FCNNs that will be presented in the following section include as additional inputs the module of the gradient of the filtered progress variable $|\nabla \bar{c}|$, and the Laplacian $\nabla^2 \bar{c}$. These operators are available in CFD solvers, and well defined on all kind of mesh, making them suitable as spatial information inputs.

The quantity used to define the error is the heat release rate, computed as
$$ \dot{Q} = - \sum_{i=1}^{N_{sp}} \Delta h_k^f \dot{\omega}_k, \qquad (14) $$

where $\Delta h_k^f$ is the enthalpy of formation of the k-th species, and $\dot{\omega}_k$ is the rate of production/consumption of the k-th species, which depends on the model used. Equation (14) highlights how the heat release rate is calculated as a weighted average of all species reaction rates. This justifies the choice of $\dot{Q}$ as the reference quantity for optimizing the FCNN. The loss function is defined as:

$$ \mathcal{L} = \| \bar{Q}_{DNS} - \bar{Q}_{FCNN} \|_2, \qquad (15) $$

where

$$ \bar{Q}_{DNS} = - \sum_{i=1}^{N_{sp}} \Delta h_k^f \dot{\omega}_k^{DNS}, \qquad (16) $$

and

$$ \bar{Q}_{FCNN} = - \sum_{i=1}^{N_{sp}} \Delta h_k^f \bar{\omega}_k^{FCNN}. \qquad (17) $$

The model was trained using an Apple M2 chip. The training process employed the Adam optimizer, with learning rate equal to $10^{-3}$, for loss minimization via gradient descent. Learning rates ranging from $10^{-1}$ to $10^{-5}$ were tested in the grid search mentioned above. The dataset was partitioned into training and testing subsets using a 70-30% split ratio. Due to the small size of the training dataset, batching was not necessary. The training was conducted in serial mode on the CPU, without GPU acceleration, using PyTorch version 2.1.2. Convergence was achieved after 1000 epochs.

# 4 Results and Discussions

We conduct an a priori investigation to determine whether the additional data-driven correction term proposed in Eq. (10) can improve the PaSR formulation of Eq. (8), indicated as ‘PaSR’ hereafter. Our assessment focuses on evaluating various strategies based on predicting the heat release rate and species reaction rates. To assess the generalizability of the FCNN model to different filter sizes, we apply various filter widths to the DNS data. All machine learning models discussed in this paper are trained on data from the slot premixed jet flame Case 1. The thermal thickness of Case 1 is $\delta_L = 110 \ \mu m$, which is 5.5 times the DNS filter width ($\delta^{DNS} = 20 \ \mu m$). Based on this thermal thickness, we tested different filter widths: $\Delta = 2, 3$ for fine LES grids, $\Delta = 6, 9$ for moderate grids, and $\Delta = 12, 18$ for coarse grids. Results from Case 1, used for both training and testing, are detailed from Sect. 4.1 to 4.4. Additionally, DNS Case 2, with a thermal thickness of $\delta_L = 292 \ \mu m = 12 \cdot \delta^{DNS}$ was analyzed using filter widths of $\Delta = 8, 12, 16, 24$. In Sect. 4.5, we evaluate one of the ML models trained on Case 1 against the data from Case 2 to assess its performance on unseen data.
## 4.1 A Priori Validation Results

This section presents a comparative analysis of five distinct models. The first two, LFR and PaSR, are fully physics-based models. The last three, referred to as M1 (Model 1), M2 (Model 2), and M3 (Model 3), adopt a hybrid approach, combining physics-based principles with data-driven techniques, see Eq. (10). These three models utilize an FCNN to predict a correction factor $\gamma_{\text{FCNN}}$ for the cell reacting fraction, see Sect. 3.3. The key distinctions among them lie in the input parameters and the filter sizes employed during the training. Table 2 provides a detailed overview of each model's characteristics and nomenclature.

In addition, a relative error metric, defined as:

$$ \Delta \text{NMSE}_{\text{log}} = \frac{\log_{10}(\text{NMSE}_{\text{ref}}) - \log_{10}(\text{NMSE}_{\text{pred}})}{\log_{10}(\text{NMSE}_{\text{ref}})} \cdot 100, \quad (18) $$

where $\text{NMSE}_{\text{ref}}$ represents the normalized mean squared error (NMSE) of the PaSR model, is used in the table. The metric is computed on the predictions of the filtered heat release rates $\bar{Q}$ at filter $\Delta = 6$ and quantifies the error reduction obtained from the use of the ML models with respect to the PaSR model baseline, see Eq. (8). The predictions are illustrated on Fig. 4a on the values of the filtered progress variable, comparing the LFR and PaSR methods and the machine learning-enhanced model M1. The LFR approach tends to overestimate the heat release rate calculation. The PaSR model addresses this issue by introducing the concept of cell reacting fraction to mitigate the overestimation inherent in the original model. The M1 model, referred as 'FCNN' in the figure, further refines this calculation using a neural network, showing an additional reduction of the estimation error.

Figure 4b displays the NMSE for some selected species. While the M1 model improves the prediction of $\bar{Q}$, using the same cell reacting fraction for every species presents limitations: the machine-learning-enhanced closure improves the prediction of some species, but it worsens others (e.g., CO and $CO_2$). Specifically, referring to the relative error metrics $\Delta \text{NMSE}_{\text{log}}$ defined in Eq. (18), the M1 model degraded the predictions of CO by 4.68%, and $CO_2$ by 7.35%. A more in-depth discussion on the limitations of using a single value of $\gamma$ can be found in several previous studies Iavarone et al. (2019); Péquin et al. (2022); Quadarella et al. (2023), where PaSR versions considering multiple cell reacting fractions for different species were introduced and shown to improve the original model. Nevertheless, the present study aims to enhance model accuracy without resorting to multiple reacting fractions, which may not assure mass balance Iavarone et al. (2019). The machine-learning-enhanced PaSR ensures instead mass conservation and, due to its simplicity, is well-suited for integration into widely utilized CFD solvers such as OpenFOAM Weller et al. (1998) or ANSYS Fluent Shabani et al. (2022).

Figure 5a showcases the different model performances at filter size $\Delta = 6$ using parity plots, where the models' predictions of $\bar{Q}$ are plotted against the reference DNS values. The points are coloured by the value of the progress variable $\bar{c}$, indicative of the reaction progress. The differences among the five models in the figure lie solely in the calculation of the cell reacting fraction, highlighting how different approaches to computing $\gamma$ impact model capabilities.
## 4.2 Generalizability Across Different Filter Sizes

Figure 5b provides a quantitative representation of the overall NMSE of heat release rate at different filter sizes, with models arranged from highest to lowest overall error. When tested on the field at $\Delta = 6$, M2 improves upon M1 by incorporating spatial information into the FCNN. This enhancement is achieved by including the modulus of the gradient of the progress variable $|\nabla \bar{c}|$ and its Laplacian $\nabla^2 \bar{c}$ as additional input parameters. Although this strategy demonstrates a large improvement on the dataset at $\Delta = 6$, the improvement is less pronounced on the finer grids. This behaviour might be attributed to the fact that both gradient and Laplacian operators inherently contain information regarding the filter size, as larger filters tend to smooth gradients during the filtering operation. The model utilizing these inputs thus becomes more specialized on the data used in the training phase but loses generalization capabilities over various filter kernels. This can become a drawback when applying the same strategy *a posteriori* since the grid size can change. By training Model 3 on data from multiple filter sizes, we managed to leverage the robustness of the model across multiple filters.

An additional comparison between M2 and M3 is shown in Fig. 6, where the parity plots indicate that training M3 on different filter sizes improved the generability of the closure over a broader range of filter sizes, including $\Delta = 12$ and $\Delta = 18$. These filter sizes were not included in the training, thus serving as an extrapolation test for M3.

A comment is necessary at this point. The 'LES-Like' fields used in the training were downsampled from the original DNS. Hence, the number of available points for training the machine learning models decreases when the filter size increases. The number of points from each field must be equal to avoid introducing biased responses in the neural network. In the present work, the highest filter size determines the number of points that can be used. This consideration led to training M2 with ≈27000 points, while M3 was trained with ≈32000 points. It is worth mentioning that the increment in the training set size is not substantial to lead to the improvements obtained with M3, but rather they may be attributed to the diversity of the training dataset. The datasets used to train the models are graphically represented in Fig. 7.

## 4.3 Model Sparsity

One of the main limitations of neural networks is their poor extrapolation capability. This strongly limits the application of such algorithms in CFD, since in practical applications the flow conditions can be different from the training data. To build trustworthy and interpretable predictive models for science and engineering, it is crucial to integrate known physical principles into machine learning models. Once the two components are integrated, the ML/data-driven part should be activated sparingly. This approach ensures that the ML component, which can be prone to extrapolation errors, only intervenes when necessary.

Key chemical reactions primarily occur across the flame front. Therefore, we expect our FCNN to adjust the cell reacting fraction most significantly in highly reactive regions. Figure 8 depicts the output, $\gamma_{FCNN}$, of the neural network on the x-y midplane for both M2 and M3. In M2, deviations from zero are frequently observed within the non-reacting region of the flame. On the other hand, M3 demonstrates a sparser pattern, featuring values near zero in the non-reacting region and higher magnitude correction values concentrated around the
flame front. In simpler terms, the model M2 fails in identifying the crucial regions that require correction.

Figure 8 also highlights how M3 identifies the major errors in the baseline PaSR model. An over-prediction of the cell reacting fraction happens near the inner flame side, while an under-prediction occurs closer to the oxidation zone. The increased sparsity of M3 was achieved without incorporating a regularization term into the loss function. Instead, it resulted from the variability introduced by merging multiple filter sizes in the input dataset. A natural question that arises is how incorporating sparsity-promoting terms or regularization would impact the final results. Further investigating this topic could offer a compelling direction for future studies.

## 4.4 Model form Limitations

The structure of Eq. (7) describing the cell reacting fraction poses some limitations. The parameter $\gamma_{PaSR}$ is uniquely a function of the Damköhler number. This limits the degrees of freedom that can be leveraged to account for the lack of subgrid information in the LES.

Figure 9 shows the functional graphs of the relationship between cell reacting fraction and Damköhler number in the baseline model (red line) and the ML-enhanced PaSR (scattered points). We also show two alternative formulations of Eq. (7): $\gamma_{PaSR^*} = 1/(1 + \text{Da}^{1/3})$ (orange line), and $\gamma_{PaSR^{**}} = 1/(1 + \text{Da})^{1/3}$ (cyan line). Following the work of Freitas et al. (2023), different exponents are applied to the mathematical formulation of $\gamma_{PaSR}$. The baseline model, which describes $\gamma_{PaSR}$ solely as a function of the Damköhler number, is limited in its accuracy compared to the ML-PaSR model that incorporates an additional degree of freedom along $\bar{c}$. The shape of the baseline model's curve $\gamma_{PaSR} = 1/(1 + \text{Da})$ aligns well with the ML-PaSR model for high values (closer to 1) of the progress variable. However, for the lower values of the progress variable, the scattered points deviate from the red line and are better represented by the alternative form $\gamma_{PaSR^*} = 1/(1 + \text{Da}^{1/3})$. A significant detail is that the monotonic decreasing behaviour of the mathematical formulations appears to be sub-optimal to describe points in the high Damköhler region. Some scattered points, with values of $\bar{c}$ approximately in the range 0.5–0.9, present an increasing trend that is not accounted for in the analytical expressions of the cell reacting fraction. An additional advantage of the ML-PaSR model over the physics-based models, is the possibility to consider values of $\gamma$ that exceed unity. The baseline model structure approximates the source terms computation, due to a lack of knowledge of the fine structures species concentration. The values of $\gamma$ that exceed unity can compensate for underestimations of the heat release rate in the LFR approach, while this is not accounted for when confining the cell reacting fraction between 0 and 1.

In Freitas et al. (2023), the best fit model was shown to be $\gamma_{PaSR^{**}} = 1/(1 + \text{Da})^{1/3}$. However, in that study, the authors computed the reaction rates integrating in time a PSR, while in the present study the instantaneous Arrhenius rates are used. This structural difference in the baseline model leads to an inability to directly compare the two studies. As already discussed in Sect. 3.2, in the present research the cell reacting fraction must account for the fact that calculating the instantaneous reaction rates from the filtered quantities completely ignores the mixing processes in the fine structures. In the original formulation of the PaSR, instead, the mixing process is modeled leveraging a PSR.

Although it is impossible to compare the two, both studies underscore the potential of proposing a novel cell reacting fraction formulation, that provides additional flexibility by adjusting the exponents in the equation. This research also highlights the potential benefits
of locally selecting the appropriate formulation between $\gamma_{PaSR}$ and $\gamma_{PaSR^*}$. Péquin et al. (2024) first used clustering methods to locally select optimal submodels to increase the accuracy of the PaSR approach. Figure 9 shows that a further advancement could be represented by utilizing clustering techniques to locally adjust the exponents of the cell reacting fraction in the PaSR model.

Figure 10 shows the cell reacting fraction $\gamma_{PaSR}$ as a function of the progress variable $\bar{c}$ and compares it with the term $\gamma_{PaSR} + \gamma_{FCNN}$. In the former case, the behaviour results from the formulation in Eq. (7). The data-driven approach overcomes the traditional mathematical limitations that restrict the cell reacting fraction to predefined boundaries, thereby unlocking the full potential of the PaSR model. Confining $\gamma$ between 0 and 1 only allows for the accounting of overpredictions, limiting the model's ability to correct underpredictions effectively. In the original PaSR model Chomiak and Karlsson (1996), a bounded value of $\gamma$ reflects the assumption that chemical reactions occur only within the fine structures, which are a fraction of the computational cell. In contrast, the PaSR formulation considered in this study provides greater flexibility for $\gamma$. The term $(\gamma_{PaSR} + \gamma_{FCNN})$ no longer describes the amount of mixture available for combustion per se but is the sum of a physical contribution and a data-driven correction, which embeds the effect of turbulent fluctuations on the Arrhenius rate. The species concentrations are still bounded by physical values, while only the species source terms are not upper bounded anymore by their 'laminar' values. Furthermore, mass conservation is ensured since $(\gamma_{PaSR} + \gamma_{FCNN})$ has the same value across all species.

## 4.5 FCNN Generalization on Unseen Data

Neural networks are typically known for their limited extrapolation capabilities, making it essential to perform tests on unseen data to evaluate the model's robustness and applicability to new cases. Model 3, which provided the smallest error among the proposed models, is tested here on unseen data from Case 2. To compare the filter widths used to process the two DNS cases, we introduce the normalised filter width $\Delta^+$ defined as

$$ \Delta^+ = \Delta \frac{\delta_L}{\delta^{\text{DNS}}} \qquad (19) $$

where $\Delta^+ = 1$ condition specifies that the LES grid cell is of a similar dimension to the flame thickness. The considered LES grids under-resolve the flame for $\Delta^+$ values above 1. On the other hand, values below 1 indicate that grid cells are thinner than the characteristic length of the flame. Table 3 presents the normalised filter width values for both cases for comparison of the results. Case 2 was filtered with values of $\Delta^+$ consistent with the ones used for Case 1. In particular, $\Delta^+ = 0.5$ represents a fine LES grid, $\Delta^+ = 1, 1.5$ correspond to a moderately fine LES filter, while $\Delta^+ = 24$ is representative of a coarse grid.

In a priori tests, the filtering operation is explicit: high-resolution data is filtered directly to separate resolved and sub-filter scales. However, during LES runs, the filtering happens implicitly, as it is determined by the numerical discretization and the grid resolution. The robustness of the model should hence be tested across different filter sizes and types. To assess the model's independency on the explicit filter type parameter, the data from Case 2 were filtered with both a Gaussian and a Box filter.

Figure 11 shows the out-of-sample testing results of model M3 on unseen data from Case 2, which differ from the training data mainly in terms of pressure and equivalence ratio, see Table 1. Both the parity and bar plots demonstrate that the hybrid M3 model consistently
<table><caption>**Table 1** Characteristics of the DNS datasets used in the present study</caption><thead><tr><th>Parameter</th><th>Case 1 Luca et al. (2019); Attili et al. (2021)</th><th>Case 2 Brouzet et al. (2021)</th></tr></thead><tbody><tr><td>Reynolds number (Re)</td><td>5600</td><td>5300</td></tr><tr><td>Jet width (H)</td><td>2 mm</td><td>1.85 mm</td></tr><tr><td>Jet mixture composition</td><td>Methane and air</td><td>Methane and air</td></tr><tr><td>Equivalence ratio (φ)</td><td>0.7</td><td>1.0</td></tr><tr><td>Karlovitz number (Ka)</td><td>23.00</td><td>15.43</td></tr><tr><td>Jet temperature</td><td>800 K</td><td>700 K</td></tr><tr><td>Background pressure</td><td>4 atm</td><td>1 atm</td></tr><tr><td>Laminar flame speed (S<sub>L</sub>)</td><td>1.0 m/s</td><td>1.85 m/s</td></tr><tr><td>Thermal thickness (δ<sub>L</sub>)</td><td>110 μm</td><td>292 μm</td></tr><tr><td>Flame characteristic time scale (τ<sub>L</sub>)</td><td>110 μs</td><td>158 μs</td></tr><tr><td>DNS minimum grid spacing (δ<sup>DNS</sup>)</td><td>20 μm</td><td>24 μm</td></tr><tr><td>Mesh size (xyz)</td><td>1440 x 962 x 256</td><td>1546 x 676 x 676</td></tr><tr><td>Full dataset grid points</td><td>≈354,000,000</td><td>≈706,000,000</td></tr><tr><td>Subset used for analysis</td><td>≈20,000,000</td><td>≈27,000,000</td></tr><tr><td>Kinetic mechanism</td><td>16 species, 72 reactions<br/>Luca et al. (2018)</td><td>14 species, 38 reactions<br/>Coffee (1984)</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Methodology to obtain the machine-learning enhanced closure for the chemical source term</figcaption></figure>

outperforms the physics-based PaSR approach across a wide range of filter sizes. Notably, M3 exceeds PaSR's performance on unseen data filtered with both box and Gaussian filter kernels, with a slight edge in accuracy on the box-filtered data. These tests show the independence of M3 on the artificial filter parameter, enhancing its potential to deliver accurate and reliable predictions in diverse LES scenarios.
<table><caption>**Table 2** Comparison of training strategies across the three ML-enhanced models introduced. M1 uses 3 input parameters and a single filter size ($\Delta = 6$), representative of a moderately fine LES grid. M2 incorporates additional spatial information via $|\nabla \bar{c}|$ and $\nabla^2 \bar{c}$. M3 extends this approach by training on multiple filter sizes. The error metric $\Delta \text{NMSE}_{\text{log}}(\bar{Q})$ represents the NMSE reduction percentage in the computation of $\bar{Q}$, in logarithmic scale, with respect to the PaSR model, and was computed on the dataset filtered with $\Delta = 6$. The training time represents the total time spent until convergence, while the testing time is associated with the time required to perform inference on a single DNS sample</caption><thead><tr><th></th><th>M1</th><th>M2</th><th>M3</th></tr></thead><tbody><tr><td>Input parameters</td><td>$\tau_c, \tau_m, \bar{c}$</td><td>$\tau_c, \tau_m, \bar{c}, |\nabla \bar{c}|, \nabla^2 \bar{c}$</td><td>$\tau_c, \tau_m, \bar{c}, |\nabla \bar{c}|, \nabla^2 \bar{c}$</td></tr><tr><td>Filter size $\Delta$</td><td>6</td><td>6</td><td>2, 3, 6, 9 (merged)</td></tr><tr><td>Training points</td><td>27,518</td><td>27,518</td><td>32,614</td></tr><tr><td>$\Delta \text{NMSE}_{\text{log}}(\bar{Q})$ [%]</td><td>-37.10</td><td>-78.31</td><td>-73.98</td></tr><tr><td>Training time [s]</td><td>2.403</td><td>4.307</td><td>4.909</td></tr><tr><td>Testing time [s]</td><td>5e-7</td><td>6e-7</td><td>6e-7</td></tr></tbody></table>

# 5 Conclusions and Future Perspectives

This study investigates the potential for enhancing the predictive capability of the partially stirred reactor (PaSR) model by integrating a data-driven correction term in a key model parameter, namely, the cell reacting fraction, which provides the *partially stirred* condition. We compared the modeled heat release rates and chemical source terms against the filtered quantities extracted from two DNS datasets of three-dimensional premixed methane/ air flames with detailed kinetics. Our evaluation included five distinct models: two fully physics-based, namely, the laminar finite rate (LFR) and PaSR, and three hybrid models, indicated as M1, M2, and M3, that combine physics-based principles with a fully connected neural network (FCNN). The results obtained with the ML-enhanced approach were consequently used for further discussion on the model's functional relationship between cell reacting fraction and Damköhler number. Alternative physics-based formulations were proposed by referring to the optimized results of the hybrid ML-physics based model.

Key findings from our *a priori* comparative analysis are the following:

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Comparative analysis of LFR, PaSR, and FCNN (M1) results, at filter size $\Delta = 6$. **a** Conditional mean of the heat release rate $\bar{Q}$, plotted as a function of progress variable $\bar{c}$. The scattered grey points in the background represent DNS values. **b** Normalized mean squared error (NMSE) in species reaction rate predictions for the three models</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 5** a Parity plot comparing the modeled heat release rate, obtained using a Gaussian filter kernel with $\Delta = 6$, to the reference DNS results. The black dashed line represents the region of perfect agreement ($\bar{Q}_{\text{model}} = \bar{Q}_{\text{DNS}}$). The grey region is included to aid visual comprehension of the plot. It encompasses points with a relative error of 20%, that is $(\bar{Q}_{\text{model}} - \bar{Q}_{\text{DNS}})/\bar{Q}_{\text{DNS}} \le 0.2$. b Normalized mean square error (NMSE) of the 5 different models, tested on a wide range of filter sizes</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Parity plots of the heat release rate. The horizontal axis represents the DNS reference values, while the vertical axis the modeled values. The performance of M2 and M3 is compared across different filter sizes</figcaption></figure>

1. The hybrid models demonstrated a significant reduction in estimation errors for the heat release rate compared to the fully physics-based models. M1, M2, and M3 effectively reduce the overpredictions observed in the LFR model and further refine the PaSR model's predictions. Notably, the relative error reduction achieved are 37%, 78%, and 74% respectively for the three models.

2. While the improvement in heat release rate predictions is significant, the enhancement in modeling species reaction rates is restricted to certain compounds. Some reactants, such as CO and CO₂, showed worsened predictions. The loss of accuracy from the PaSR model is of 4.68% in the prediction of CO, and 7.35% in the prediction of CO₂.
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Scatter plot showing the training dataset distribution for the machine-learning-enhanced models M1, M2, and M3. The x-axis represents the base 10 logarithm of $\tau_c$ and the y-axis the base 10 logarithm of $\tau_m$. Green dots correspond to the data used in models M1 and M2, while purple dots correspond to the data used to train model M3 (see Table 2). Kernel Density Estimate (KDE) plots on the top and right sides illustrate the marginal distributions of $\log_{10}(\tau_c)$ and $\log_{10}(\tau_m)$ for the respective datasets. The mixing timescale can assume values lower than the Kolmogorov scale, because they are associated with the dissipation of the scalar $\bar{c}$, describing the chemical reaction progress, as outlined in Eq. (13)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8** 2D spatial representation of the FCNN output on the x – y midplane, tested on filter size $\Delta = 2$</figcaption></figure>

This result highlights the limitations of using a unique cell reacting fraction for all the species.

3. The generalization capabilities of the FCNN across different filter sizes are limited when training on a single value of $\Delta$. This shortcoming is more pronounced with M2 although the model is more accurate on the filter size used for the training ($\Delta = 6$). The addition of spatial information in the input parameters seems to be beneficial for the prediction of seen data but shows a negative effect on extrapolation capabilities.

4. M3, trained on data from multiple filter sizes, demonstrated better generalizability and robustness across a broader range of grid sizes, addressing the model specialization
<figure><img src="image_10.png" /><figcaption>**Fig. 9** Plot of the cell reacting fraction as a function of the Damköhler number. The scattered points represent $\gamma_{FCNN} + \gamma_{PaSR}$, obtained with the Machine Learning correction computed with M3 on filter size $\Delta = 6$. The continuous red line shows the definition of $\gamma$ from the baseline PaSR formulation, $\gamma_{PaSR} = \frac{1}{1+Da}$. The continuous orange and cyan lines represent two alternative formulations of the cell reacting fraction proposed by Freitas et al. (2023).</figcaption></figure>

issue observed for M2. The improvements were obtained without explicitly informing the network of the filtering kernel. Additionally, the M3 model demonstrated strong generalization capabilities on out-of-sample data, showing adaptability across different flame conditions.

5. The study highlights the importance of modeling a sparse ML correction for practical applications in CFD. M3 displayed a more localized adjustment of the cell reacting

<figure><img src="image_11.png" /><figcaption>**Fig. 10** Behavior of the cell reacting fraction as a function of the progress variable. The left plot illustrates $\gamma_{PaSR}$, calculated using the mathematical formulation. The right plot displays the corrected cell reacting fraction, which incorporates the additional term $\gamma_{FCNN}$, computed with M3, into the baseline model</figcaption></figure>

<table><caption>**Table 3** Comparison of filter widths between case 1 and case 2</caption><thead><tr><th>Case</th><th>δ<sub>L</sub><br>[μm]</th><th>δ<sup>DNS</sup><br>[μm]</th><th>Δ|<sub>Δ+</sub>≈0.3</th><th>Δ|<sub>Δ+</sub>≈0.5</th><th>Δ|<sub>Δ+</sub>≈1</th><th>Δ|<sub>Δ+</sub>≈1.5</th><th>Δ|<sub>Δ+</sub>≈2</th><th>Δ|<sub>Δ+</sub>≈3</th></tr></thead><tbody><tr><td>1</td><td>110</td><td>20</td><td>2</td><td>3</td><td>6</td><td>9</td><td>12</td><td>18</td></tr><tr><td>2</td><td>292</td><td>24</td><td>-</td><td>6</td><td>12</td><td>18</td><td>24</td><td>-</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Fig. 11** Results of the DNS Case 2 testing: the top section displays parity plots comparing the modeled heat release rate to the reference DNS filtered values at a filter size of $\Delta = 12$. The grey region includes points with a relative error of 20%. The bottom section shows the normalized mean square error (NMSE) for three models (LFR, PaSR, and M3) across varying filter sizes. Plots in (a) are obtained by filtering the DNS with a Gaussian filter kernel, while those in (b) are based on a box filter kernel</figcaption></figure>
fraction, concentrating higher absolute values of the correction near the flame front and maintaining near-zero values in non-reacting regions, thus improving the model's physical interpretability and reliability. This improvement in M3's sparsity was obtained without additional regularization in the loss function. Instead, it was shown to be the consequence of combining multiple filter sizes in the training dataset.

6. The data-driven approach successfully overcame the mathematical constraints of the original PaSR model, which limited the cell reacting fraction to specific boundaries defined by the hypotheses on the term $γ_{PaSR}$. The hybrid models unlocked additional degrees of freedom for this parameter, enhancing the model's capability to account for subgrid information.

Overall, this work demonstrates that integrating data-driven corrections into physics-based models can substantially enhance the accuracy of predictions in reactive flow simulations of premixed flames. Future research will aim to further refine these hybrid models by incorporating multiple progress variables as input features and multi-objective loss functions that account for errors in species predictions. Additionally, investigating the transferability of the neural network correction in a *posteriori* applications would be valuable. Another promising direction is to incorporate multiple cell reacting fractions into the baseline equations, as in the generalized PaSR formulation proposed in Quadarella et al. (2023), allowing for the extraction of a greater amount of information from DNS data. Our research not only validates the potential of ML-enhanced PaSR models but also opens up new possibilities for improving model accuracy through innovative formulations and clustering techniques. These advancements ultimately contribute to better predictive capabilities of turbulent combustion closures, paving the way for more precise simulations of turbulent reacting flows and the consequential development of sustainable combustion technologies.

## Appendix A PaSR Model Selection

The PaSR model approximates TCI by integrating in time a PSR to compute the reaction rates. In the present study, however, we relied on an alternative formulation of PaSR as a baseline model, outlined in Eq. (8). In this section, we demonstrate that not only the chosen formulation reduces the computational time, but it also increases the accuracy of the ML-enhanced model. Considering Eq. (6) as the initial model, the ML-enhanced model takes the following form

$$ \bar{\omega}_k^{\text{PaSR}} = (\gamma_{\text{PaSR}} + \gamma_{\text{FCNN}}) \cdot \bar{\omega}_k^{\text{QLFR}}, \quad (A1) $$

where $\bar{\omega}_k^{\text{QLFR}}$ represents the Quasi Laminar Finite Rate (QLFR) source terms

$$ \bar{\omega}_k^{\text{QLFR}} = \frac{\bar{\rho}(Y_k^* - Y_k^0)}{\tau^*}. \quad (A2) $$

We refer to the model of Eq. (A1) as QLFR-based, while the model of Eq. (10) is here termed LFR-based. To determine the most suitable closure, we tested both formulations by training the Model 1 (M1) architecture on each approach. Figure 12 illustrates the conditional mean of the heat release rate $\bar{Q}$ given $\bar{c}$ obtained by the LFR, QLFR, LFR-based Model 1, and QLFR-based Model 1 formulations in comparison with the DNS results. The LFR-based Model 1 provides a NMSE on the heat release rate of 2.58 · 10⁻³ while the
<figure><img src="image_13.png" /><figcaption>Fig. 12 Conditional mean of the heat release rate $\bar{Q}$ as a function of the progress variable $\bar{c}$. Results from the DNS, M1 trained on the LFR-based formulation (FCNN$_{LFR}$), M1 trained on the QLFR-based formulation (FCNN$_{QLFR}$), the LFR and QLFR models are reported</figcaption></figure>

QLFR-based Model 1 gives a NMSE error of $4.83 \cdot 10^{-3}$. The higher accuracy, together with the lower complexity and computational time that characterize the LFR-based PaSR, contributed to the decision of selecting the LFR-based formulation.

## Appendix B Alternative Definitions of the Progress Variable

Two alternative definitions of the progress variable are here tested and compared to the oxygen-based progress variable to assess the sensitivity of the model to this input parameter. We use a progress variable based on temperature

$$ \bar{c}_T = \frac{\bar{T} - \bar{T}_u}{\bar{T}_b - \bar{T}_u}, \qquad (B3) $$

and one based on the concentration of the products H$_2$O and CO$_2$

$$ \bar{c}_{H_2O+CO_2} = \frac{\bar{Y}_{H_2O} + \bar{Y}_{CO_2} - (\bar{Y}_{H_2O} + \bar{Y}_{CO_2})_u}{(\bar{Y}_{H_2O} + \bar{Y}_{CO_2})_b - (\bar{Y}_{H_2O} + \bar{Y}_{CO_2})_u}, \qquad (B4) $$

with the subscripts b and u referring to the burned and unburned conditions, respectively.

Figure 13 shows the correlation between the three formulations of the progress variables proposed. It can be deduced that the three variables are almost perfectly correlated with each other.

The ML-enhanced model that achieved the best performance, i.e., Model 3 (M3), was trained on both temperature-based and products-based progress variables. The testing results obtained on filter size $\Delta = 6$ are compared in Fig. 14, reporting the modeled heat release rate values and the reference DNS filtered values. Model 3 does not show a significant sensitivity to the progress variable definition, with the progress variable based on the main combustion products performing slightly better than the others.
<figure><img src="image_14.png" /><figcaption>**Fig. 13** Two-dimensional bin plot illustrating the relationships among the three progress variables analyzed in the sensitivity study. Pearson and Spearman correlation coefficients are provided for each variable pair. The colormap indicates the density of points within each bin. The subfigures display the correlations between: **a** $\bar{c}_{H_2O+CO_2}$ and $\bar{c}_T$, **b** $\bar{c}_{O_2}$ and $\bar{c}_{H_2O+CO_2}$, and **c** $\bar{c}_{O_2}$ and $\bar{c}_T$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14** Predictions of $\bar{Q}$ by the ML-enhanced PaSR against the reference DNS values for M3 trained on (a) the oxygen-based progress variable $\bar{c}_{O_2}$, (b) the products-based progress variable $\bar{c}_{H_2O+CO_2}$, (c) the temperature-based progress variable $\bar{c}_T$</figcaption></figure>

**Acknowledgements** The ENCODING project has received funding from the European Union's Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101072779. The results of this publication/presentation reflect only the author(s) view and do not necessarily reflect those of the European Union. The European Union cannot be held responsible for them. A. Péquin acknowledges the financial support of the Fonds National de la Recherche Scientifique (FRS-FNRS) through an ASPIRANT fellowship.

**Author Contributions** All authors conceptualised the study. L.P., A.P., and R.F. conducted the numerical study. L.P. drafted the manuscript text. A.P., R.F., S.I., H.P. and A.P. supervised and reviewed the manuscript.

**Funding** Not applicable

**Data Availability** The data that support the findings of this study are accessible from the corresponding authors upon reasonable request.

**Code Availability** Not applicable
**Declarations**

**Conflict of interest** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Ethics Approval and Consent to Participate** Not applicable

**Consent for publication** Consent for publication has been obtained from all authors. This manuscript does not include any identifiable personal data that requires separate consent.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**References**

Brankovic, A., Porter, L., McKinney, R., Ouyang, H., Kennedy, J., Madabhushi, R., Colket, M.: Comparison of measurements and predictions of flow in a gas turbine engine fuel nozzle. In: 38th Aerospace Sciences Meeting and Exhibit. American Institute of Aeronautics and Astronautics, Reston, VA, USA (2000)

Poinsot, T., Veynante, D.: Theoretical and Numerical Combustion, vol. 28. R.T. Edwards Inc, Philadelphia, PA, USA (2005)

Amaduzzi, Ruggero, Ferrarotti, Marco, Parente, Alessandro: Strategies for hydrogen-Enriched methane flameless combustion in a Quasi-industrial furnace. Front. Energy Res. (2021). https://doi.org/10.3389/fengr.2020.590300

Huang, X., Tummers, M.J., Veen, E.H., Roekaerts, D.J.E.M.: Modelling of MILD combustion in a lab-scale furnace with an extended FGM model including turbulence-radiation interaction. Combust. Flame **237**, 111634 (2022)

Wilcox, D.: Turbulence Modeling for CFD, 3rd edn. DCW Industries Inc, La Canada, CA, USA (2006)

Pope, S.B.: Ten questions concerning the large-eddy simulation of turbulent flows. New J. Phys. **6**, 35–35 (2004)

Veynante, D., Vervisch, L.: Turbulent combustion modeling. Prog. Energy Combust. Sci. **28**(3), 193–266 (2002)

Moin, P., Apte, S.V.: Large-eddy simulation of realistic gas turbine combustors. AIAA J. **44**(4), 698–708 (2006)

Sabaut, P.: Large Eddy Simulation for Incompressible Flows. Scientific computation. Springer, Berlin, Germany (2001)

Garnier, E., Adams, N., Sabaut, P.: Large Eddy Simulation for Compressible Flows, 2009th edn. Scientific computation. Springer, Dordrecht, Netherlands (2009)

Peters, N.: Cambridge Monographs on Mechanics: Turbulent Combustion. Cambridge University Press, Cambridge, England (2000)

Pitsch, H.: Large-eddy simulation of turbulent combustion. Annu. Rev. Fluid Mech. **38**(1), 453–482 (2006)

Meneveau, C., Katz, J.: Scale-invariance and turbulence models for large-eddy simulation. Annu. Rev. Fluid Mech. **32**(1), 1–32 (2000)

Pérez Arroyo, C., Dombard, J., Duchaine, F., Gicquel, L., Odier, N., Exilard, G., Richard, S., Buffaz, N., Démolis, J.: Large-Eddy Simulation of an Integrated High-Pressure Compressor and Combustion Chamber of a Typical Turbine Engine Architecture. In: Proceedings of the ASME Turbo Expo 2020: Turbomachinery Technical Conference and Exposition. Volume 2C: Turbomachinery, pp. 02–35058. American Society of Mechanical Engineers, New York, NY, USA (2020). doi: https://doi.org/10.1115/ GT2020-16288
Fureby, C., Sahut, G., Ercole, A., Nilsson, T.: Large eddy simulation of combustion for high-speed air-breathing engines. Aerospace **9**(12), 785 (2022)

Choi, H., Moin, P.: Grid-point requirements for large eddy simulation: Chapman's estimates revisited. Phys. Fluids **24**(1), 011702 (2012)

Péquin, A., Evans, M.J., Chinnici, A., Medwell, P.R., Parente, A.: The reactor-based perspective on finite-rate chemistry in turbulent reacting flows: a review from traditional to low-emission combustion. Appl Energy Combust. Sci. **16**, 100201 (2023)

Magnussen, B.F., Hjertager, B.H.: On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion. Symp. Int. Combust. **16**(1), 719–729 (1977)

Magnussen, B.: Modeling of nox and soot formation by the eddy dissipation concept. International Flame Research Foundation First Topic Oriented Meeting Amsterdam (1989)

Gran, I.R., Magnussen, B.F.: A numerical study of a bluff-body stabilized diffusion flame. part 2. influence of combustion modeling and finite-rate chemistry. Combust. Sci. Technol. **119**(1–6), 191–217 (1996)

Chomiak, J., Karlsson, A.: Flame liftoff in diesel sprays. Symp. Int. Combust. **26**(2), 2557–2564 (1996)

Golovitchev, V.I., Nordin, N., Chomiak, J.: Neat dimethyl ether: Is it really diesel fuel of promise? In: SAE Technical Paper Series. SAE International, 400 Commonwealth Drive, Warrendale, PA, United States (1998)

Golovitchev, V.I., Nordin, N., Jarnicki, R., Chomiak, J.: 3-D diesel spray simulations using a new detailed chemistry turbulent combustion model. In: SAE Technical Paper Series, vol. 109, pp. 1391–1405. SAE International, 400 Commonwealth Drive, Warrendale, PA, United States (2000)

Golovitchev, V.I.: Revising “old” good models: Detailed chemistry spray combustion modeling based on eddy dissipation concept. In: SAE Technical Paper Series. SAE International, 400 Commonwealth Drive, Warrendale, PA, United States (2001)

Karlsson, A.: Modeling Auto-ignition, Flame Propagation and Combustion in Non-stationary Turbulent Sprays. Chalmers university of technology, ??? (1995)

Trisjono, P., Pitsch, H.: Systematic analysis strategies for the development of combustion models from DNS: a review. Appl. Sci. Res. **95**(2–3), 231–259 (2015)

Iavarone, S., Péquin, A., Chen, Z.X., Doan, N.A.K., Swaminathan, N., Parente, A.: An a priori assessment of the partially stirred reactor (PaSR) model for MILD combustion. Proc. Combust. Inst. **38**(4), 5403–5414 (2021)

Péquin, A., Iavarone, S., Malpica Galassi, R., Parente, A.: The partially stirred reactor model for combustion closure in large eddy simulations: physical principles, sub-models for the cell reacting fraction, and open challenges. Phys. Fluids (1994) **34**(5), 055122 (2022)

Quadarella, E., Péquin, A., Stagni, A., Parente, A., Faravelli, T., Im, H.G.: A generalized partially stirred reactor model for turbulent closure. Proc. Combust. Inst. **39**(4), 5329–5338 (2023)

Swaminathan, N., Parente, A.: Machine Learning and Its Application to Reacting Flows: ML and Combustion. Springer, ??? (2023)

Ihme, M., Chung, W.T., Mishra, A.A.: Combustion machine learning: principles, progress and prospects. Prog. Energy Combust. Sci. **91**(101010), 101010 (2022)

Bode, M., Gauding, M., Lian, Z., Denker, D., Davidovic, M., Kleinheinz, K., Jitsev, J., Pitsch, H.: Using physics-informed super-resolution generative adversarial networks for subgrid modeling in turbulent reactive flows (2019)

Bode, M., Gauding, M., Lian, Z., Denker, D., Davidovic, M., Kleinheinz, K., Jitsev, J., Pitsch, H.: Using physics-informed enhanced super-resolution generative adversarial networks for subfilter modeling in turbulent reactive flows. Proc. Combust. Inst. **38**(2), 2617–2625 (2021)

Nista, L., Schumann, C.D.K., Grenga, T., Attili, A., Pitsch, H.: Investigation of the generalization capability of a generative adversarial network for large eddy simulation of turbulent premixed reacting flows. Proc. Combust. Inst. **39**(4), 5279–5288 (2023)

Deep learning for presumed probability density function models: T. Henry de Frahan, M., Yellapantula, S., King, R., S. Day, M., W. Grout, R. Combust. Flame **208**, 436–450 (2019). https://doi.org/10.1016/j.combustflame.2019.07.015

Chen, Z.X., Iavarone, S., Ghiasi, G., Kannan, V., D'Alessio, G., Parente, A., Swaminathan, N.: Application of machine learning for filtered density function closure in MILD combustion. Combust. Flame **225**, 160–179 (2021). https://doi.org/10.1016/j.combustflame.2020.10.043

Yang, H., Kobayashi, T., Iavarone, S., Massey, J.C., Chen, Z.X., Minamoto, Y., Swaminathan, N.: Towards a generalised artificial neural network for sub-grid filtered density function closure in turbulent combustion. Appl Energy Combust. Sci. **14**, 100142 (2023). https://doi.org/10.1016/j.jaecs.2023.100142

Doan, N.A.K., Swaminathan, N., Minamoto, Y.: DNS of MILD combustion with mixture fraction variations. Combust. Flame **189**, 173–189 (2018). https://doi.org/10.1016/j.combustflame.2017.10.030
Lapeyre, C.J., Misdariis, A., Cazard, N., Veynante, D., Poinsot, T.: Training convolutional neural networks to estimate turbulent sub-grid scale reaction rates. Combust. Flame **203**, 255–264 (2019). https://doi.org/10.1016/j.combustflame.2019.02.019

Nikolaou, Z.M., Chrysostomou, C., Vervisch, L., Cant, S.: Progress variable variance and filtered rate modelling using convolutional neural networks and flamelet methods. Flow Turbulence Combust **103**, 485–501 (2019). https://doi.org/10.1007/s10494-019-00028-w

Péquin, A., Iavarone, S., Malpica Galassi, R., Parente, A.: Supervised clustering for optimal sub-model selection in reactor-based models. Flow Turbulence Combust **112**, 931–955 (2024)

Freitas, R.S.M., Péquin, A., Galassi, R.M., Attili, A., Parente, A.: Model identification in reactor-based combustion closures using sparse symbolic regression. Combust. Flame **255**, 112925 (2023)

Jigjid, K., Minamoto, Y., Doan, N.A.K., Tanahashi, M.: SGS reaction rate modelling for mild combustion based on machine-learning combustion mode classification: development and a priori study. Proc. Combust. Inst. **39**(4), 4489–4499 (2023)

Luca, S., Attili, A., Lo Schiavo, E., Creta, F., Bisetti, F.: On the statistics of flame stretch in turbulent premixed jet flames in the thin reaction zone regime at varying Reynolds number. Proc. Combust. Inst. **37**(2), 2451–2459 (2019)

Attili, A., Luca, S., Denker, D., Bisetti, F., Pitsch, H.: Turbulent flame speed and reaction layer thickening in premixed jet flames at constant Karlovitz and increasing Reynolds numbers. Proc. Combust. Inst. **38**(2), 2939–2947 (2021)

Brouzet, D., Talei, M., Brear, M.J., Cuenot, B.: The impact of chemical modelling on turbulent premixed flame acoustics. J. Fluid Mech. (2021). https://doi.org/10.1017/jfm.2020.1184

Coffee, T.P.: Kinetic mechanisms for premixed, laminar, steady state methane/air flames. Combust. Flame **55**(2), 161–170 (1984)

Luca, S., Al-Khateeb, A.N., Attili, A., Bisetti, F.: Comprehensive validation of skeletal mechanism for turbulent premixed methane-air flame simulations. J. Propuls. Power **34**(1), 153–160 (2018)

Sabelnikov, V., Fureby, C.: Extended LES-PaSR model for simulation of turbulent combustion. In: Progress in Propulsion Physics. EDP Sciences, Les Ulis, France (2013)

Nordin, P.A.N.: Complex chemistry modeling of diesel spray combustion, 1–55 (2001)

Hastie, T., Tibshirani, R., Friedman, J.: The Elements of Statistical Learning, 2nd edn. Springer series in statistics. Springer, New York, NY (2009)

Raman, V., Pitsch, H.: A consistent LES/filtered-density function formulation for the simulation of turbulent flames with detailed chemistry. Proc. Combust. Inst. **31**(2), 1711–1719 (2007)

Iavarone, S., Cafiero, M., Ferrarotti, M., Contino, F., Parente, A.: A multiscale combustion model formulation for NOx predictions in hydrogen enriched jet flames. Int. J. Hydrog. Energy **44**(41), 23436–23457 (2019). https://doi.org/10.1016/j.ijhydene.2019.07.019. (Accessed **2022-11-15**)

Weller, H.G., Tabor, G., Jasak, H., Fureby, C.: A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput. Phys. **12**(6), 620–631 (1998)

Shabani, S., Majkut, M., Dykas, S., Wiśniewski, P., Smolka, K., Cai, X., Zhang, G.: Numerical analysis of the condensing steam flow by means of ansys fluent and in-house academic codes with respect to the capacity for thermodynamic assessment. Front. Energy Res. **10**, 902629 (2022)

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Authors and Affiliations

Lorenzo Piu<fnref n="1" />,² · Arthur Péquin<fnref n="1" />,² · Rodolfo S. M. Freitas³ · Salvatore Iavarone<fnref n="1" />,²,⁴ · Heinz Pitsch⁵ · Alessandro Parente<fnref n="1" />,²,⁶

☒ Lorenzo Piu
lorenzo.piu@ulb.be

☒ Alessandro Parente
alessandro.parente@ulb.be

<footnote n="1">¹ Aero-Thermo-Mechanics Laboratory, École Polytechnique de Bruxelles, Université Libre de Bruxelles, 1050 Brussels, Belgium</footnote>
2 Brussels Institute for Thermal Energy (BRITE), 1050 Brussels, Belgium

3 School of Engineering and Materials Science, Queen Mary University of London, London E1 4NS, UK

4 EM2C laboratory, CentraleSupélec, CNRS, Université Paris-Saclay, 91190 Gif-sur-Yvette, France

5 Institut für Technische Verbrennung, RWTH Aachen University, 52062 Aachen, Germany

6 WEL Research Institute, 1300 Wavre, Belgium