A Priori Direct Numerical Simulation Assessment of MILD Combustion Modelling in the Context of Reynolds Averaged Navier–Stokes Simulations

A priori Direct Numerical Simulation (DNS) assessment of mean reaction rate closures for reaction progress variable in the context of Reynolds Averaged Navier–Stokes (RANS) simulations has been conducted for MILD combustion of homogeneous (i.e., constant equivalence ratio), methane-air mixtures. The reaction rate predictions according to statistical (e.g., presumed probability density function), phenomenological (e.g., eddy-break up (EBU), eddy dissipation concept (EDC) and the scalar dissipation rate (SDR) based approaches), and flame surface description (e.g., Flame Surface Density) based closures are compared. The performance of the various reaction rate closures has been assessed by comparing the models’ predictions to the corresponding quantities extracted from DNS data. It has been found that the usual presumed probability density function (PDF) approach using the beta-function predicts the PDF of the reaction progress variable in homogenous mixture MILD combustion throughout the flame brush for all cases considered here provided that the scalar variance is accurately predicted. The accurate estimation of scalar variance requires the solution of a modelled transport equation, which depends on the closure of Favre-averaged SDR. A linear relaxation based algebraic closure for the Favre-averaged SDR has been found to capture the behaviour of the Favre-averaged SDR in the current homogenous mixture MILD combustion setup. It has been found that the EBU, SDR and FSD-based mean reaction rate closures do not adequately predict the mean reaction rate of the reaction progress variable for the parameter range considered here. However, a variant of the EDC closure, with model coefficients expressed as functions of micro-scale Damköhler and turbulent Reynolds numbers, has been found to be more successful in predicting the mean reaction rate of reaction progress variable compared to other modelling methodologies for the range of turbulence intensities and dilution levels considered here.


Introduction
Moderate or Intense Low oxygen Dilution (MILD) combustion (Cavaliere and de Joannon 2004;Wünning and Wünning 1997;Katsuki and Hasegawa 1998) is a methodology that can achieve low emissions without compromising the thermal efficiency.There are two definitions for MILD combustion based on the perfectly stirred reactor (PSR) theory and thus these definitions are applicable for homogenous mixture combustion (in this analysis, the term "homogenous mixture" refers to mixtures with no equivalence ratio variation).The first definition suggests that there are no critical ignition or extinction points in MILD combustion (Cavaliere and de Joannon 2004).The second, and most commonly used definition, states that MILD combustion occurs when the initial temperature of the reactants is higher than the mixture self-ignition temperature T 0 > T ign and the maximum temperature rise due to combustion is less than the mixture self-ignition temperature ΔT| max < T ign (Cavaliere and de Joannon 2004).Experimental techniques and direct numerical simulations (DNS) have provided significant physical insights into the fundamental understanding of MILD combustion.Experimental investigations using Planar Laser-Induced Fluorescence images of OH radicals (OH-PLIF) revealed the existence of flame fronts in MILD combustion, while the distributed nature of combustion was observed from temperature measurements using Rayleigh thermometry (Plessing et al. 1998).DNS investigations (Minamoto et al. 2013;Minamoto et al. 2014a, b;Minamoto andSwaminathan 2014, 2015;Awad et al. 2021) attributed the appearance of distributed combustion to the significant amount of reaction zone interactions.Moreover, it has been found that the extent of reaction zones interaction is significantly affected by O 2 concentration (dilution level), while the effect of turbulence intensity remains insignificant (Minamoto et al. 2014a, b).There have been several attempts at the numerical modelling of MILD combustion and the closure of the mean reaction rate was achieved using either flamelet (Göktolga et al. 2017;Coelho and Peters 2001;Huang et al. 2022) or Eddy dissipation (Aminian et al. 2011;De et al. 2011;Parente et al. 2016;Awad and Eldrainy 2021) approaches.These studies provided mostly satisfactory agreement with experimental data in terms of the predictions of mean values of temperature and major species mass fractions.However, when it comes to the predictions of peak temperature and minor species mass fractions (e.g., CO and OH concentrations), significant discrepancies were reported.Thus, despite its advantages, MILD combustion is still not well-understood, and its modelling remains a challenging task.To date, relatively limited effort has been focused on the modelling of MILD combustion in contrast to the vast body of literature on conventional combustion approaches.
Conventional combustion models can be classified into three broad categories, which are: models based on flame surface description (e.g., Flame Surface Density (FSD)), statistical approaches (e.g., presumed PDF approach), and phenomenological models (e.g., eddy-break up (EBU) model, eddy dissipation model (EDC) and the scalar dissipation rate (SDR) based mean reaction rate closure).This paper aims to compare the performance of these combustion models in terms of the prediction of the mean reaction rate of reaction progress variable in the context of Reynolds Averaged Navier-Stokes (RANS) simulations using a priori DNS analysis for methane-air, homogenous mixture (i.e., constant equivalence ratio) MILD combustion.
Several previous studies (Huang et al. 2022;Chen et al. 2017;Minamoto and Swaminathan 2015) considered mean reaction rate closure in MILD combustion using the presumed PDF approach.Minamoto and Swaminathan (2015) reported that the reaction progress variable ( c ) distribution in MILD combustion can be adequately approximated by the beta-PDF and several numerical investigations for MILD combustion used the beta-PDF along with different canonical solutions for tabulating the reaction rate (Huang et al. 2022;Chen et al. 2017).The main modelling challenge in the presumed PDF approach based on beta-PDF is the accurate prediction of the second moment of the scalar in question.One of the unclosed terms in the second moment transport equation involves the scalar dissipation rate (SDR), which is either modelled by an algebraic expression (Kolla et al. 2009;Chakraborty and Swaminathan 2011;Dunstan et al. 2013;Gao et al. 2014;Gao et al. 2015a, b) or a modelled transport equation is solved for SDR (Mantel and Borghi 1994;Chakraborty et al. 2008;Gao et al. 2015a, b).Closures of the SDR have been investigated and assessed in conventional homogeneous mixture combustion and interested readers are referred to Chakraborty et al. (2011) for further information in this regard.However, this aspect has not been analysed for homogeneous mixture MILD combustion.Therefore, it is worthwhile to assess the use of the beta-PDF and different algebraic closures for the SDR in MILD combustion.It can be appreciated from the foregoing that the assessments of the validity of beta-PDF parameterisation of the reaction progress variable and SDR closures in homogenous mixture MILD combustion are necessary for evaluating the capability of the presumed PDF approach in predicting the mean reaction rate of reaction progress variable.
A phenomenological model that has been extensively used in the literature for the modelling of mean reaction rate in MILD combustion is the eddy dissipation concept (EDC) (Magnussen 1981), which is an extended version of the modified eddybreakup (EBU) model (Magnussen and Hjertager 1977).The EDC approach with its standard model constants shows a good agreement with experimental data for conventional combustion (Gran and Magnussen 1996).However, when it is compared with experimental data for MILD combustion, discrepancies have been observed, and caseby-case adjustment of model constants was needed for accurate predictions (Aminian et al. 2011;De et al. 2011;Parente et al. 2016;Awad and Eldrainy 2021).The standard EDC model assumes that all chemical reactions occur in "reacting fine structures" and can be treated as Perfectly Stirred Reactor (PSR).Thus, the model constants (for the reactor residence time/timescale and the size/mass fractions of the fine structures) were proposed based on the classical cascade model, which assumes that all the dissipation processes occur at the smallest scales and are developed for high turbulent Reynolds number Re t situations with a clear separation between large and small scales (Ertesvåg 2022).However, in MILD combustion chemical reaction occurs over a wide range of scales, and the turbulent Reynolds numbers involved are generally low (Minamoto et al. 2014a, b;Minamoto et al. 2013;Minamoto andSwaminathan 2014, 2015).Thus, there is a need to revise the EDC model constants to account for the short cascade at low turbulent Reynolds numbers (Parente et al. 2016;Evans et al. 2019;Lewandowski et al. 2020).To account for the effects of low Reynolds and Damköhler numbers, Parente et al. (2016) proposed expressions for the EDC model parameters (related to the reactor timescale and the mass fraction of fine structures) as functions of the local micro-scale Damköhler and turbulent Reynolds numbers.These expressions were then revised further by Evans et al. (2019), Lewandowski et al. (2020), Romero-Anton et al. (2020) and Fordoei et al. (2021).The proposed functional expressions for reactor timescale and mass fraction of fine structures resulted in better predictions for the peak temperature and minor species compared to the EDC with constant (standard or modified) coefficients, particularly at low and medium values of turbulent Reynolds number Re t (Aminian et al. 2011;De et al. 2011).A recent theoretical investigation by Ertesvåg (2022) analysed the behaviours of the functional expressions for model parameters for the reactor timescale and mass fraction of fine structures proposed by Parente et al. (2016), Evans et al. (2019), Lewandowski et al. (2020) and Romero-Anton et al. (2020) in the context of RANS and provided important insights into the effects of micro-scale Damköhler and turbulent Reynolds numbers, as well as the limiting behaviours of the coefficients related to the reactor timescale and mass fraction of fine structures.It is worth noting that the EDC model is usually used to estimate the mean reaction rate for individual species in a multispecies system, but the current study only assesses the predictive abilities of the EDC model for the mean reaction rate of a suitably chosen reaction progress variable.
Another phenomenological model is the SDR based mean reaction rate closure.Modelling the mean reaction rate based on the SDR is not commonly used in MILD combustion because the model is based on an infinitely fast chemistry assumption.However, the closure of the mean reaction rate with the help of SDR has been found to hold reasonably well in the thin reactions zone regime combustion (Peters 2000) of homogeneous mixtures where the assumption of infinitely fast chemistry does not hold (Chakraborty and Cant 2011).The FSD based modelling, which belong to the flame surface description approaches, is also not commonly used for the modelling of the mean reaction rate in MILD combustion but is one of the most used methods of mean reaction rate closure for premixed combustion (Chakraborty and Cant 2011;Cant and Bray 1988;Cant et al. 1990;Candel et al. 1990;Trouvé and Poinsot 1994).A recent DNS investigation (Awad et al. 2021) compared the behaviour of the reaction progress variable gradient |∇c| for homogeneous mixture MILD combustion with turbulent pre- mixed flames in the thin reaction zones regime and indicated that FSD models for turbulent premixed flames could be extended for homogenous mixture MILD combustion which warrant assessment of the FSD based closures in the case of homogeneous mixture MILD combustion.
From the previous presentation, it can be appreciated that it is worthwhile to assess the use of the beta-PDF and different algebraic closures for the SDR in MILD combustion as well as the performance of the EBU, EDC, FSD and SDR-based closures in predicting the mean reaction rate of reaction progress variable for homogenous mixture MILD combustion in the context of RANS.To achieve these objectives, a Direct Numerical Simulation (DNS) database for homogenous mixture MILD combustion of CH 4 has been considered which includes cases with different turbulence intensities and different O 2 concentrations (4.8% and 3.5% by volume) at an equivalence ratio of ϕ = 0.8.The main objectives for the current analysis are: (a) To carry out a priori DNS assessment of different closures for mean reaction rate of reaction progress variable in the context of RANS for homogeneous mixture (i.e., constant equivalence ratio) MILD combustion, and (b) To investigate the role of turbulence intensities and dilution levels on the performance of these models.
The paper will be structured as follows.Different turbulent combustion models and SDR closures will be discussed in the next section.This is followed by a discussion of the DNS dataset.Then, the results will be presented and discussed.Finally, the conclusions are summarised.

Combustion Modelling Approaches
The focus of the current analysis is the prediction of the reaction rate of the reaction progress variable in the context of RANS for homogeneous mixture MILD combustion.In this analysis, the reaction progress variable ( c ) is defined as (Minamoto et al. 2014a, b;Mina- moto et al. 2013;Minamoto andSwaminathan 2014, 2015;Awad et al. 2021): where Y F is the fuel (i.e., CH 4 ) mass fraction and subscripts R and P refer to values in unburned reactants and fully burned products, respectively.Thus, the reaction rate of the reaction progress variable ωc is given by: where ωF is the fuel's reaction rate.The following subsections include descriptions of the various closures for the mean reaction rate of reaction progress variable ωc analysed in this paper.Here q , q = q∕ and q �� = q − q are the Reynolds averaged, Favre-averaged and Favre fluctuation of a general quantity q , respectively.

Presumed PDF approach
In the context of tabulated chemistry, the presumed PDF approach provides a closure for the mean reaction rate of the reaction progress variable ωc .This is done by integrating the tabu- lated values of ωc with a presumed PDF for the reaction progress variable c .The closure of ωc in the presumed PDF approach is given as (Huang et al. 2022; Chen et al. 2017): The beta function has been used in different MILD combustion RANS investigations due to its flexibility (Huang et al. 2022;Chen et al. 2017).The PDF according to the beta distribution is defined as (Poinsot and Veynante 2005): where Γ(x) is the gamma function and it is defined as (Poinsot and Veynante 2005): In Eq. 4, a and b are defined as (Poinsot and Veynante 2005): Thus, the presumed PDF approach requires the knowledge of the Favre-averaged reaction progress variable variance c′′2 .The transport equation for c′′2 is given by (Chakraborty and Swaminathan 2011;Poinsot and Veynante 2005): where D is the reaction progress variable diffusivity.The first, third and fifth terms on the right-hand side of Eq. 7 are unclosed and require modelling but the present analysis focuses on the last term on the right-hand side, which includes the SDR involving scalar fluctuations (Kolla et al. 2009 ; Chakraborty and Swaminathan 2011): In the context of passive scalar mixing, ̃ c is often modelled using the linear relaxation model, given as (Kolla et al. 2009;Chakraborty et al. 2011): where ∕ ρ are the Favre-averaged turbulent kinetic energy and its dissipation rate respectively, and C is a proportionality constant of the order of unity (Kolla et al. 2009;Chakraborty et al. 2011).Kolla et al. (2009) proposed an SDR closure for high Damköhler numbers ( Da ≫ 1 ) premixed flames, which includes the effects of chemical and turbulent timescales.Their proposed model also includes the physics which affects the scalar gradient alignment with strain rate eigenvalues and captures the competition between the strain rate induced by heat release and turbulent straining.The model by Kolla et al. (2009) is given as: where th is the thermal flame thickness defined as th = T ad − T 0 ∕max |∇T| L (with T, T ad and T 0 being the instantaneous, adiabatic flame and reactant temperatures, respectively and the subscript 'L' refers to the values in the corresponding 1D unstretched premixed flame) and S L is the unstretched laminar burning velocity.In Eq. 10, ) dc is a thermochem- ical parameter and the other model parameters are defined as, 4 and = 6.7 with Ka L being the local Karlovitz number defined as Kolla et al. 2009;Chakraborty and Swaminathan 2011).

The Eddy Break-Up (EBU) and Eddy Dissipation Concept (EDC) closures
According to the EBU model the mean reaction rate of reaction progress variable is expressed in the following manner (Magnussen and Hjertager 1977): (7) where ỸL is the mean concentration of the limiting reactant ỸL = min( ỸF , Ỹo ∕s), s is the sto- chiometric oxygen to fuel mass ratio and A is a model constant, which takes a value of 4.0 (Magnussen and Hjertager 1977).
The EDC model includes the effect of finite rate chemistry by separating the computational cell into two regions: reacting fine structures, and non-reacting surrounding fluids.The model assumes that all chemical reactions occur at the fine structure which are assumed to be of the same order as the Kolmogorov scale.The fine structure residence time * and its mass fraction γ λ are given as (Gran and Magnussen 1996; Parente et al.  2016).
where C and C are model constants and ̃ is the Favre-averaged kinematic viscosity.At these scales, the turbulent mixing is assumed to be considerably fast in comparison to the chemical reaction rate.Thus, in the present study, the fine structures are modelled as a perfectly stirred reactor (PSR) according to the following set of equations (Gran and Magnussen 1996): where Y * i and ω * i are the species fine structure mass fraction and reaction rate, respectively, h is the specific enthalpy and p is the pressure.The initial conditions for the above govern- ing equations are the mean mass fractions in the reactant mixture ⟨Y i ⟩ , initial temperature T o and the residence time res .The mean mass fractions ⟨Y i ⟩ are estimated based on the vol- ume average of the species in the DNS initial field and res is taken as the mean convective time from the inflow to the outflow DNS boundaries.This set of equations is integrated in time using COSILAB (2007) until attaining a steady state solution following the same process previously used by Minamoto and Swaminathan (2015).
The mean reaction rate closure for species i is modelled in the EDC methodology in the following manner (Parente et al. 2016;Chakraborty et al. 2008): The standard values for C and C are 0.4083 and 2.1377 respectively (Gran and Mag- nussen 1996).It is worth noting that Eq. 17 is valid for every species in the chemical mechanism.However, the present analysis focuses only on the closure of mean reaction rate of reaction progress variable, as the FSD and SDR-based mean reaction rate closures are usually expressed for reaction progress variable c .Thus, to compare the EDC model ( 12) performance with the other combustion models considered here, the EDC mean reaction rate is expressed in terms of the mean reaction rate of reaction progress variable.Using i = F (i.e., fuel mass fraction) in Eq. 17 and using the resulting expression in Eq. 2 yields the expression for the mean reaction rate of reaction progress variable ωc .
It is also worth noting that the EDC model suffers from inherent deficiencies at low turbulent Reynolds numbers.This can be illustrated by expressing the fine structure mass fraction as a function of turbulent Reynolds number Re t = k2 ∕ ν ε (i.e., = C Re t −1∕4 ) (Ertesvåg 2022).Under low Re t , approaches unity (or assumes values even greater than unity), resulting in singularities (or unrealistic values) in the EDC mean reaction rate closure (i.e., Eq. 17).To avoid these undesired consequences (i.e., singularities and unrealistic mean reaction rates), the fine structure mass fraction is bounded by an arbitrary value within the range of 0.7-0.95(i.e.,   < 0.7 − 0.95 ) in practical applications (Ertesvåg 2022).To address these deficiencies, Parente et al. (2016) revised the standard model and proposed functional expressions for C and C based on the turbulent Reynolds and the local Kolmogorov Damköhler numbers which can capture some of the specific features of MILD combustion (e.g., distributed burning, smooth scalar gradients, etc.).The proposed functional expressions by Parente et al. (2016) and Evans et al. (2019) are given as follows: where Re t = � k 2 ∕( ε ν) is the local turbulent Reynolds number and Da is the local Kolmogo- rov Damköhler number defined as the ratio between the Kolmogorov mixing time scale t to the chemical time scale t c (Parente et al. 2016):   In the present analysis, the Kolmogorov mixing time is defined as t  = ( ν∕ ε) 1∕2 and the chemical time scale is defined as et al. 2016) where i is the species used to define the reaction progress variable.Lewandowski et al. (2020) applied more stringent assumptions when revising the classical energy cascade model and proposed modifications to the functional expressions of C and C in the following manner: The functional expressions by Parente et al. (2016), and their modification by Lewandowski et al. (2020), extend the validity of the EDC at low Re t conditions.How- ever, constraints on C and C cannot be ignored to avoid singularities or unrealistic predictions (Ertesvåg 2022).

The Flame Surface Density (FSD) and SDR-based closures
The closure of ωc is obtained in the Flame Surface Density (FSD) modelling approach (Chakraborty and Cant 2011; Cant and Bray 1988;Cant et al. 1990;Candel et al. (1990); Trouvé and Poinsot 1994) through the following expression: where Σ gen = |∇c| is the generalised FSD (Boger et al. 1998) and S d s = S d |∇c| ∕ |∇c| is the surface weighted value of the density-weighted displacement speed S d = (Dc∕Dt)∕|∇c| .In the context of RANS, the assumption ωc ≫ ∇.(D∇c) is often invoked and S d s is often modelled as o S L (i.e., S d s = 0 S L ) (Chakraborty and Cant 2011;Cant and Bray 1988;Cant et al. 1990;Candel et al. 1990;Trouvé and Poinsot 1994;Boger et al. 1998) where o is the unburned gas density.Under this assumption, the stretch factor, defined as I o = ωc /( o S L Σ gen ) , is assumed to be unity (i.e., I o = 1.0 ) and the mean reaction rate is closed according to the following expression (Chakraborty and Cant 2011;Cant and Bray 1988;Cant et al. 1990;Candel et al. (1990); Trouvé and Poinsot 1994;Boger et al. 1998): Bray (1980) proposed the following closure for the mean reaction rate ωc using SDR for high Damköhler number (i.e., Da ≫ 1 ) combustion, which was subsequently shown to hold even in the thin reaction zones regime combustion based on a scaling analysis (Chakraborty and Cant, 2011): (Bray 1980).

Numerical Implementation
The simulations have been conducted using the compressible DNS code SENGA2 (Cant 2012) which solves the standard governing equations of energy, momentum, mass, and species for turbulent reacting flows.In SENGA2, the spatial derivatives at internal grid points are approximated using a tenth-order central difference scheme and the order of accuracy gradually decreases to a one-sided fourth-order scheme at non-periodic boundaries (Cant 2012).The time advancement has been achieved using a fourth order Runge-Kutta scheme.A skeletal methane-air chemical mechanism consisting of 16 species and 25 reactions (Smooke and Giovangigli 1991) has been used for the present analysis.The Navier-Stokes Characteristic Boundary Condition (NSCBC) methodology has been used for specifying the boundary conditions (Poinsot and Lele 1992).The left x-direction boundary is a tur- bulent inflow with specified density, species, and velocity components, whereas the right x-direction boundary is taken to be partially non-reflecting outflow.Periodic boundary conditions have been assigned for all other boundaries.A cubic domain of size 10 mm × 10 mm × 10 mm is used for the simulations considered in this paper with a cartesian (23) grid of 252 × 252 × 252 used to discretise the domain.This grid spacing ensures that th is resolved by using at least 12 grid points and the Kolmogorov length scale is resolved by at least 1.5 grid points for all cases investigated.Two thermochemical conditions are investigated in the present analysis corresponding to high O 2 concentration (low dilution case, 4.8% O 2 by volume) and low O 2 concentration (high dilution case, 3.5% O 2 by vol- ume).The thermochemical conditions, which include the reactants mole fractions and initial temperature, are given in Table 1.The low and high dilution cases are simulated for two turbulence intensities u � ∕S L at the same length scale ratio l∕ th where l is the integral length scale and u ′ is the root-mean-square turbulent velocity fluctuations.The turbulence conditions, which include the initial velocity and length scale ratios along with the values of Damköhler number Da = lS L ∕u � th and Karlovitz number Ka = u � ∕S L 3∕2 l∕ th −1∕2 , are given in Table 2. Simulating a real-life MILD combustor using DNS is still not possible with the available computational power.Thus, the current DNS of MILD combustion is conducted according to the methodology proposed by Minamoto et al. (2013).In this methodology, the MILD combustion process is split into two phases: 1) the initial and inflow field preparation phase and 2) the combustion phase.The first phase, which mimics the exhaust gas recirculation, provides the initial and inflow fields for the second phase.The procedure for preparing the first phase is summarised in the following steps (Minamoto et al. 2014a, b;Minamoto et al. 2013;Minamoto andSwaminathan 2014, 2015;Awad et al. 2021): 1.The pseudo-spectral method of Rogallo (1981) has been used to create a freely decaying homogenous isotropic turbulence following the Batchelor-Townsend spectrum (Batchelor and Townsend 1948).2. 1D laminar premixed flames for the low and high dilution cases under the conditions shown in Table 1 are simulated, and the species mass fractions are then specified as functions of the methane mass fraction-based reaction progress variable.3. The pseudo-spectral methodology proposed by Eswaran and Pope (1988) is used to generate a bimodal distribution for c with prescribed integral length scale l c = 1.50 mm .A bimodal distribution for the various species fields is then developed using the laminar flame parameterisation in step 2. 4. The generated freely decaying isotropic turbulence and bimodal scalar fields are then allowed to evolve without reaction for about 1 mixing time scale (l∕u � ) at atmospheric pressure and an unburned gas temperature of 1500 K (i.e., T 0 = 1500K ).The mean and variance of the c fields at the end of this step are ⟨c⟩ ≈ 0.50 and ⟨c �2 ⟩ ≈ 0.09 respectively.The temperature in this pre-processed mixture exhibits a variation of about ±0.3% of its mean value.
The scalar fields from phase one are then used as an initial field for phase two of the simulation as well as being fed into the domain at the inflow boundary condition on the left x-boundary with a mean velocity U in = 20 m∕s .The simulations have been continued for 2.5 through pass times (i.e., t sim = 2.5L x ∕U mean ), which amounts to 11.0l∕u � and 22.0 l∕u � ( 6.0 l∕u � and 12.0 l∕u � ) for the LD (HD) case with initial u � ∕S L = 4.0 and 8.0, respectively.To extract Reynolds/Favre averaged and Favre fluctuation quantities from DNS data, the simulation data has been time-averaged based on time sequences between 1.0 and 2.5 through pass times following the procedure previously adopted by Minamoto et al. (2013).

Results and Discussion
Figure 1 shows the instantaneous views of the c = 0.8 isosurface and slices of the tempera- ture field for both LD and HD cases at u � ∕S L = 4.0.The instantaneous views of the c = 0.8 iso-surface can be considered as representative of the flame surface since the maximum heat release rate for the present thermochemistry takes place at around c = 0.8 .It can be seen from Fig. 1 that MILD combustion cases exhibit significant interactions of reaction zones for both LD and HD cases.This is consistent with previous MILD combustion DNS studies (Minamoto et al. 2014b).Figure 1 also shows that the temperature rise in the HD case is lower than that in the LD case.This is expected since the higher concentrations of H 2 O and CO 2 in HD cases compared to LD cases (see Table 1) result in a higher mixture specific heat capacity and the lower O 2 concentration, which in turn give rise to a lower heat release rate.

Parameterisation of the PDF and Closure of SDR of Reaction Progress Variable
The PDFs of c for different values of c (i.e., exemplarily for c = 0.35, 0.6 and 0.80 to illus- trate the behaviour across the flame brush) obtained from the DNS and the predictions of beta-function parameterisation are shown in Fig. 2 for both LD and HD cases at initial turbulence intensities of u � ∕S L = 4.0 and 8.0.In the present analysis, the time aver- aging between 1.0 and 2.5 through pass times is used to obtain the Reynolds and Favreaveraged quantities and the samples at every location within the domain at different time instants during the sampling period are used to calculate the PDF of c for a given value of c .It can be seen from Fig. 2 that the beta-function captures the qualitative and quantita- tive behaviour of the PDF of c throughout the flame brush for all cases investigated.Thus, the beta-function can be considered as an adequate representative for the reaction progress variable PDF in homogenous mixture MILD combustion processes.This observation is consistent with previous DNS analysis by Minamoto and Swaminathan (2015).However, where, O c is the burning mode contribution.In the limit of infinitely fast chemistry, where the PDF of c can be assumed to behave as a double-delta function with impulses at c = 0 and c = 1.0 , the contribution O c can be neglected and therefore, c′′2 can be mod- elled as c 1 − c .
Figure 3 shows c′′2 and c 1 − c conditioned upon c for all cases considered.It can be seen from Fig. 3 that c 1 − c assumes higher values compared to c′′2 in all investigated cases suggesting a non-negligible contribution for the burning mode contribution O c .This is expected since the PDFs of c do not exhibit bimodal distributions for the cases con- sidered here (see Fig. 2).Thus, one can conclude that c′′2 cannot be modelled as c 1 − c in MILD combustion and a transport equation for c′′2 is needed.However, the closure of c′′2 based on its transport equation depends on the modelling of the Favre-averaged scalar dissipation rate (SDR) Ñc , which will be discussed next in this paper.
Figure 4 shows the normalised Favre averaged SDR Ñc × th ∕S L conditioned upon c calculated using Eqs. 9 and 10 after adding the resolved SDR component (i.e., Ñc = Eq.9/10 + D∇c.∇c ) for all cases investigated.It can be seen from Fig. 4 that the SDR decreases with increasing c for all cases investigated because c′′2 decreases with increasing c (see Fig. 3).Moreover, it can be observed that the SDR assumes higher magnitudes for the higher turbulence intensity cases, which is due to the increase in the scalar fluctuation gradients as a result of increasing the turbulence intensity.It can also be seen from Fig. 4 that the model proposed by Kolla et al. (2009) (i.e., Eq. 10) underpredicts the SDR at low values of c for all cases investigated and the level of disagreement is more significant for the HD cases (3.5% O 2 ).However, the model shows a reasonably good prediction for the SDR at � c > 0.7 for the lower dilution cases (4.8% O 2 ).To explain the previous observations, the local Dam- köhler number Da L = t f ∕t c (with t f being the integral flow time scale t f = 0.25 k∕̃ = l t ∕u rms based on k = 3u 2 rms ∕2 and ̃ = 0.37u 3 rms ∕l t with u rms and l t being the local values of root- mean-square velocity and integral length scale, respectively, and t c is defined in the same manner as shown in the context of Eq. 20) conditioned upon c for the LD case with inlet turbulence intensities of u � ∕S L =4.0 and 8.0 are shown in Fig. 5.It can be seen from Fig. 5 that the Da L assumes small values (i.e., Da L < 1.0 ) for � c < 0.7 , but Da L > 1 is obtained for � c > 0.7 .As the model proposed by Kolla et al. (2009) (i.e., Eq. 10) is strictly valid for Da L ≫ 1 , the model predicts the SDR behaviour for � c > 0.7 and underpredicts the SDR in regions where Da L assumes values smaller than unity (i.e., � c < 0.7 ).The prediction of the SDR by a linear relaxation model (i.e., Eq. 9), which is used for passive scalar mixing situations, is also shown in Fig. 4 for optimised values of the model parameter C .It has been found that the linear relaxation model (i.e., Eq. 9) captures the qualitative nature of the variation of Ñc and a satisfactory level of quantitative agreement can be obtained when C is parameterised in the following manner based on a regression analysis: It can be seen from Fig. 4 that the linear relaxation model (i.e., Eq. 9) with optimised C is able to capture the trending behaviour of the SDR for all cases investigated suggesting that the behaviour of the SDR in MILD combustion and passive scalar mixing share some similarities.The small change in density and temperature in MILD combustion leads to weakening of the effects of strain rate induced by the flame normal acceleration and dilatation rate and therefore the effects of chemical timescale in the SDR modelling can be ignored and the micro-mixing is principally determined by turbulent mixing, which also justifies the satisfactory performance of the linear relaxation model especially for low Damköhler number (i.e., Da < 1 ) condition.
In practical combustors both MILD and conventional combustion modes can coexist.For example, Li et al. (2014) stated that Da ranges between 0.01 and 5.35 under MILD combustion conditions, which is consistent with the current findings (see Fig. 5).Thus, it is worth proposing an algebraic SDR closure that is valid for both high Da premixed combus- tion and low Damköhler number MILD combustion.Therefore, Eqs. 9 and 10 have been combined based on the suggestion by Mura et al. (2007) in the following manner so that the model expression remains valid for both high Da premixed combustion and low Dam- köhler number MILD combustion: Here s * = c��2 ∕ c 1 − c is the segregation factor.For Da ≫ 1 where c��2 ≈ c 1 − c the segregation factor will be close to 1.0 and the first term on the right hand side of Eq. 28 will be the leading order term, whereas for low Damköhler number (i.e., Da < 1 ) combus- tion one obtains s * ≪ 1 , and thus the second term on the right-hand side of Eq. 28 will be ( 27) the leading order term.The prediction of Eq. 28 is also shown in Fig. 4, which shows that Eq. 28 prediction follows the same trend as that of the linear relaxation model prediction for all cases investigated, which is expected since c′′2 assumes smaller values than c 1 − c in MILD combustion (see Fig. 3).The coefficient of determination (i.e., R-squared) between the discussed SDR models (i.e., Eq. 9, Eq. 10 and Eq.28) and Ñc extracted from the DNS data for all cases investigated are shown in Fig. 6.It can be seen from Fig. 6 that the coefficient of determination increases slightly with increasing the turbulence intensity for all cases investigated apart from the predictions of Eq. 10 in HD cases.However, it significantly decreases with increasing the dilution level.This, in turn, suggests that all the discussed models are significantly affected by the level of dilution, and they are expected to show better performance at low dilution conditions.The coefficient of determination for Eq. 9 and Eq.28 also confirms the above discussion (i.e., passive scalar SDR models can captures the qualitative nature of Ñc in MILD combustion) with a satisfactory R-squared values ( R 2 > 0.5) for all cases investigated.

Assessment of Mean Reaction Rate Closures by EBU and EDC Methodologies
Figure 7 shows the prediction of the mean reaction rate ωc using the models listed in Table 3 for both LD and HD cases with turbulence intensities of u � ∕S L = 4.0 and 8.0.It can be seen from Fig. 7 that the modified EBU and the EDC with standard model constants (i.e., C = 0.4083 and C = 2.1377) overpredict the mean reaction rate ωc by an order of magnitude for all cases investigated.In MILD combustion the low O 2 concentration slows down the chemical reaction, thus models that are based on infinitely fast chemistry assumptions (i.e., the modified EBU) are no longer valid, particularly when using coefficients calibrated for conventional combustion processes.This explains the overprediction of the mean reaction rate by the modified EBU model.The overprediction of ωc also occurs for the conventional EDC model with standard model constants.This is expected since the values of the standard model constants are based on the classical energy cascade model, but a clear separation of scales, as assumed in the model derivation, is not realised for moderate values of the turbulent Reynolds number (such as those in the current DNS cases).Moreover, the effects of chemical reactions in MILD combustion are not confined to the fine structures (Swaminathan and Chakraborty 2021).Different RANS studies (De et al. 2011;Evans et al. 2015) suggested increasing C from 0.4083 to 3.00 and decreasing C from 2.1377 to 1.0 to account for the dilution effect in MILD combustion.The predic- tions of the EDC model with these modified constants are also shown in Fig. 7, which shows that the EDC model with modified model constants underpredicts the mean reaction Fig. 6 Coefficient of determination between SDR models and DNS for all cases investigated rate for all cases investigated.Increasing C from 0.4083 to 3.0 overpredicts the fine structure residence time * (see Eq. 17), thus underpredicts the mean reaction rate for all cases investigated.Figure 7  It is worth noting that Eq. 30 allows for a maximum value of C of 1.0 instead of 2.1377 suggested in the analysis by Ertesvåg (2022).The upper limit of C is taken to be 1.0 in this analysis since this value was shown to be more appropriate for MILD combustion by previous studies by De et al. (2011) and Evans et al. (2015).It can be seen from Fig. 7 that the EDC model with C and C , as given by Eqs. 29 and 30, shows the best qualitative agreement with the DNS data among all the variants.The effect of dilution changes locally in MILD combustion, thus the chemical time scale and consequently Da (see Fig. 5) change across the flame brush.The C and C expressions include the effect of Damköhler number variations and low turbulent Reynolds number effects and thus, this model performs well in comparison to the EDC variants with fixed model constants.In order to gain further insights into the performance of the M4 model, the PDFs of Da and Re t in the region given by 0.1 ≤ c ≤ 0.9 are shown in Fig. 8a and b.It can be seen from Fig. 8a and b that the most probable values of Da and Re t are about 0.05-0.1 and 50 -100, respectively, for all cases considered here.Equations 28 and 29 and the analysis by Ertesvåg (2022) indicate that C is determined by 0.5 Da 2 Re T + 1 −0.5 while C assumes a value of unity (due to the limiting condition in Eq. 30) for the range of Da and Re t obtained in the cases considered here.This can be further substantiated from the PDFs of C and C (evaluated according to Eqs. 18 and 19) in the region given by 0.1 ≤ c ≤ 0.9 , which are also shown in Fig. 8c and d, respectively.This suggests that for the parameter range considered here, C can be taken to be unity alongside the C expression given by Eq. 29, but a different outcome might be obtained for a different parameter range.It has indeed been found that using 2.1377 as the upper limit of C worsens the performance of M4 model in comparison to the predictions obtained from Eq. 30 for the cases considered here and thus are not shown explicitly in this paper.Similar to Eqs. 28 and 30, the functional expressions of C and C in the EDC model by Lewandowski et al. (2020) (i.e., M5 in Table 3) can be subjected to the following combination: (31) C = 0.5 Da 2 0.09Re T + 1 −0.5  De et al. (2011;Evans et al. 2015) Parente et al. (2016); Evans et al. (2019) (EDC-C and C depend on local properties) C in this paper given by Eq. 29 C in this paper given by Eq. 30 M5 (originally proposed by Lewandowski et al. (2020) (EDC-C and C depend on local properties) C in this paper given by Eq. 31 C in this paper given by Eq. 32 The predictions of M5 with C and C given by Eqs. 31 and 32 are also shown in Fig. 7 and the performance of this model remains comparable to the M4 model predictions with limiting conditions for the model parameters.Equation 32 also yields C ≈ 1.0 for the range of Da and Re t obtained in the cases considered here because of the limiting condi- tions (i.e., Eq. 22 leads to C  > 1.0 as shown in Ertesvåg ( 2022)).Thus, the differences between M4 and M5 model predictions originate due to the differences in C values.The coefficient of determination (i.e., R-squared) between the models shown in Table 3 and the mean reaction rate estimated from the DNS data for all cases investigated are shown in Fig. 9.It can be seen from Fig. 9 that the EDC model variants M4 and M5 exhibit good performance ( R 2 > 0.85) regardless of the change in turbulent intensity or dilution level.Moreover, Fig. 9 shows that the other models considered (i.e., M1, M2 and M3) do not show a satisfactory performance ( R 2 < 0.3) for all cases investigated.where K Σ is a constant.Thus, both models are expected to show similar behaviours for the prediction of ωc (i.e., ωc decreases as c increases).However, the discrepancies between the predicted ωc and that calculated from the DNS are relatively lower for the SDR-based mean reaction rate closure compared to the FSD based mean reaction rate modelling.Moreover, it can be observed that the prediction of the FSD and SDR-based modelling approaches remains insensitive to the change in turbulence intensity and O 2 concentration level.

Conclusions
The prediction of the mean reaction rate of reaction progress variable for homogeneous mixture (i.e., constant equivalence ratio) MILD combustion of methane in the context of Reynolds Averaged Navier-Stokes (RANS) has been investigated based on a priori DNS analysis.Three different combustion modelling approaches (i.e., statistical, phenomenological, and flame surface description approaches) have been analysed.The assessment has been performed using DNS data for different initial turbulence intensities (i.e., u � ∕S L = 4.0 and 8.0) and different dilution levels (4.8% and 3.5% O 2 by volume) at equivalence ratio = 0.8 .The main findings of this analysis are: • In the context of the statistical approach of mean reaction rate closure of reaction progress variable, the presumed PDF method requires parameterisation of c .It has been found that (33) εc = K Σ S L Σ gen a beta-function based parameterisation captures the PDF of c both qualitatively and quan- titatively throughout the flame brush for all cases investigated.However, the closure of scalar variance is needed to use the beta-function parameterisation of the reaction progress variable PDF.The algebraic closure of scalar variance using the BML expression, which is derived based on a presumed bimodal PDF with impulses at unburned reactants and fully products under high Damköhler number conditions, has been found to be inadequate for MILD combustion.This indicates that a modelled transport equation for scalar variance is needed, which requires the closure of Favre-averaged scalar dissipation rate (SDR).It has been found that the SDR in a homogenous mixture MILD combustion can be satisfactorily modelled using a linear relaxation-based model usually used for passive scalar mixing SDR closure suggesting a similarity between both regimes.Moreover, a new algebraic SDR closure expression, utilising a bridging function based on a segregation factor, that changes the leading order contribution in the model according to the local combustion mode (i.e., MILD or premixed combustion) is proposed.
• The flame surface density (FSD), which belongs to the flame surface description modelling approaches, overpredicts (underpredicts) the mean reaction rate at low (high) values of Favre-averaged reaction progress variable c for all cases investigated.This suggests that S d s cannot be approximated as o S L and accurate modelling of the stretch factor I 0 is needed for extending the FSD-based closure to MILD combustion.
• The SDR-based mean reaction rate closure of the reaction progress variable, which is one of the phenomenological modelling approaches considered here, showed the same qualitative behaviour as the FSD modelling approach.This closure method has been found not to provide an adequate prediction of the mean reaction rate of the reaction progress variable for the parameter range considered here.• The other phenomenological modelling concepts, such as the eddy dissipation concept (EDC) and Eddy-Break Up (EBU) closures have also been investigated in the present analysis.The modified EBU and the EDC models with standard values of model parameters overpredict the mean reaction rate by an order of magnitude.However, modifying the EDC model parameters related to the reactor timescale and fine-structure mass fraction based on the previous suggestion by De et al. (2011) and Evans et al. (2015) leads to underprediction of the mean reaction rate of the reaction progress variable.The EDC model variants with functional expressions for the model parameters in terms of micro-scale Damköhler number and turbulent Reynolds number (including appropriate limiting condition) have been found to exhibit reasonably good agreement with the mean reaction rate of reaction progress variable extracted from DNS data for the range of parameters considered here.
The EDC model variants that provide the satisfactory prediction of the mean reaction rate for reaction progress variable needs to be validated further for higher values of turbulent Reynolds number and these model variants need to be implemented in RANS simulations for the purpose of a posteriori assessment.It is worth noting that the performance of the EDC model could be different for different species but the assessment of EDC model performance for all species in a multispecies system is beyond the scope of the current analysis.Furthermore, the validity of the modelling methodologies discussed in this paper for LES of MILD combustion is yet to be ascertained based on a priori DNS analysis.However, the model parameters in the EDC model in the context of RANS were proposed based on modelling of the full energy cascade, which will not be valid in the context of LES due to the resolution of a part of the energy cascading process.This suggests that the model parameters of the EDC model may need to be altered from the values suggested for RANS in the context of LES.It is worth noting that both FSD and SDR closures were originally proposed for conventional premixed combustion under the assumption of fast chemistry (i.e., high values of Damköhler number), despite these models showing some possibility to provide reasonable predictions within the thin reaction zones (Peters 2000) combustion where the assumption of infinitely fast chemistry does not hold (Chakraborty and Cant 2011).However, low Damköhler number effects in MILD combustion are likely be stronger than the conventional thin reaction zones regime combustion, which means that additional considerations need to be accounted for extending FSD and SDR-based mean reaction rate closures for MILD combustion.This includes accurate modelling of surface-averaged values of the density-weighted displacement speed for the FSD-based closure and accounting for the additional terms in the SDR-based reaction rate due to low Damköhler number effects, as suggested by Bray et al. (2011).Some of these aspects will form the basis of future investigations.

Fig. 1
Fig. 1 Instantaneous views of the c = 0.8 isosurface (top row) and slices showing the non-dimensional temperature = (T − T o )∕(T p − T o ) at the mid X-Y plane normalised by its maximum value in the shown slice (bottom row) for LD (left) and HD (right) cases for the inlet turbulence intensity of u � ∕S L = 4.0

Fig. 2
Fig. 2 PDFs of c for c = 0.35, 0.60 and 0.80 from DNS (solid line) and − function parameterisation (dashed line) for LD (up) and HD (bottom) cases for inlet turbulence intensities of u � ∕S L =4.0 (left) and 8.0 (right)

Fig. 3
Fig. 3 Variation of c′′2 (black) and c 1 − c (red) conditioned upon c for LD (up) and HD (bottom) cases for inlet turbulence intensities of u � ∕S L = 4.0 (left) and 8.0 (right)

Fig. 4
Fig. 4 Variation of the normalised Favre averaged SDR Ñc × th ∕S L conditioned upon c obtained from DNS along with the predictions of Eqs. 9, 10 and 28 for LD (up) and HD (bottom) cases in the case of inlet turbulence intensities of u � ∕S L = 4.0 (left) and 8.0 (right)

Fig. 5
Fig. 5 Variation of Da L conditioned upon c for the LD cases with inlet turbulence intensities of u � ∕S L = 4.0 (left) and 8.0 (right) also shows the prediction of the EDC model with variables C and C (i.e., M4 model in Table 3) based on the functional expressions proposed by Parente et al. (2016) and Evans et al. (2019) (see Eqs. 18 and 19) but with the following combination informed by a recent analysis by Ertesvåg (2022):

Fig. 7
Fig. 7 Variation of the normalised mean reaction rate ωc ×  th ∕ 0 S L conditioned upon c obtained from DNS along with the predictions of models summarised in Table 3 for LD (top row) and HD (bottom row) cases in the case of inlet turbulence intensities of u � ∕S L = 4.0 (left) and 8.0 (right).Note: M1 and M2 are multiplied by a factor of 0.1

Fig. 8
Fig. 8 PDFs of a Da , b Re t , c C (according to Eq. 18) and d C (according to Eq. 19) in the region given by 0.1 ≤ c ≤ 0.9 for both LD (black lines) and HD (red lines) cases in the case of inlet turbulence intensities of u � ∕S L = 4.0 (solid line) and 8.0 (broken line)

Fig. 9 Fig. 10
Fig.9Coefficient of determination between the modified EBU and EDC models (see Table3) and the mean reaction rate ωc extracted from DNS for all cases investigated.Note: unrecognized R 2 values are close to zero R 2 ≈ 0)

Fig. 11
Fig. 11 Variation of the stretch factor I o = ωc ∕ o S L Σ gen conditioned upon c for all cases investigated

Table 1
Thermochemical conditions for the 1D unstretched laminar premixed flames

Table 3
List of considered EBU and EDC models