Surface Density Function and Its Evolution in Homogeneous and Inhomogeneous Mixture n-Heptane MILD Combustion

ABSTRACT

Moderate or intense low-oxygen dilution (MILD) combustion is a combustion technique with potential to simultaneously improve thermal efficiency and reduce emissions. This paper focuses on the mean behavior of the reactive scalar gradient characterized by the surface density function (SDF = magnitude of the reaction progress variable gradient) and its evolution for exhaust gas recirculation (EGR) type, homogeneous and inhomogeneous mixture n-heptane combustion under MILD conditions using Direct Numerical Simulations (DNS) with reduced chemical mechanism. Two oxygen concentrations have been considered here for the homogeneous mixture case, namely, 3.0% and 4.5% by volume. The characteristics of the SDF and displacement speed were also analyzed in the reaction and propagating-flame-dominated regions of the domain to illustrate the effect of combustion mode on the mean variations of these quantities. The mean values of the SDF in turbulent MILD combustion cases are found to be smaller compared to the corresponding laminar flame cases. Moreover, the effect of mixture inhomogeneity on the SDF variation is found to be marginal for the parameters considered here. It is found that increasing the dilution factor reduces the percentage of flame thickening with respect to the corresponding laminar flame thickness under identical turbulence conditions. The effects of dilatation rate are found to be weak in the studied cases due to the expected low heat release under MILD conditions which leads to weak thermal expansion effects. The effective flame normal and tangential strain rates are found to be dominated by the additional flame normal and tangential strain rates due to flame propagation, respectively. The effective flame normal strain rate has positive values across the flame which promotes flame thickening and could explain the decrease in the SDF values, while the effective flame tangential strain rate is negative throughout the flame and indicates a fractional reduction in flame area that could be attributed to flame surface interactions. The reaction-dominated regions give rise to reduced level of the SDF in the homogeneous mixture cases and increased level of $S_d$ in all cases compared to those calculated in the propagating-flame region.

KEYWORDS:

• MILD combustion
• Surface density function
• Flame surface area
• Tangential strain rate
• Normal strain rate
• Displacement speed
• Curvature

Introduction

Combustion processes still maintain a major role in primary energy supply, particularly in applications requiring high energy density such as heavy vehicle transportation, industrial furnaces, and gas turbines. This trend is expected to continue for the foreseeable future (Masri 2021). As the environmental regulations are becoming increasingly stringent, it is necessary to develop a new generation of highly efficient and environmentally friendly combustors. In this context, moderate or intense low oxygen dilution (MILD) combustion has been demonstrated to have the potential for achieving both high-energy efficiency and ultra-low emissions (Cavaliere and de Joannon 2004; Li et al. 2011; Perpignan, Gangoli Rao, and Roekaerts 2018; Wünning and Wünning 1997).

In MILD combustion, the reacting mixture is preheated above its auto-ignition temperature ($T_r>T_{ign}$), while the maximum temperature rise remains smaller than the auto-ignition temperature ($\mathrm{\Delta }T{|}_{\mathrm{m}\mathrm{a}\mathrm{x}}<T_{ign}$) (Wünning and Wünning 1997). This can be achieved by recirculating the hot flue gas back into the combustion chamber in furnaces and boilers, which inherently increases their thermal efficiency.

Experimental investigations of MILD combustion using Planar Laser-Induced Fluorescence images of OH radicals (OH-PLIF) have shown the presence of flame fronts, and indications of distributed burning were revealed through temperature measurements using Rayleigh thermometry (Dally, Karpetis, and Barlow 2002; Plessing, Peters, and Wünning 1998; Özdemir and Peters 2001). Minamoto and coworkers conducted 3D DNS studies of Exhaust Gas Recirculation (EGR) type, homogeneous methane–air mixture, MILD combustion, where they investigated the distribution of species, temperature, and reaction rate fields, morphological and topological structures, and scalar gradient statistics in terms of temperature and equivalent OH-PLIF signal (Minamoto and Swaminathan 2014; Minamoto et al. 2013, 2014a, 2014b). They concluded that the distributed combustion seen in experimental visualization (Plessing, Peters, and Wünning 1998; Özdemir and Peters 2001) can be attributed to the interaction of thin-reaction zones (Minamoto et al. 2013). Doan and coworkers (Doan and Swaminathan 2019a, 2019b; Doan, Swaminathan, and Minamoto 2018) extended the previous analysis for stratified methane–air mixture (inhomogeneous) MILD combustion and provided important insights into the flame structure, reaction zone topology, and markers of combustion mode. These DNS analyses (Doan, Swaminathan, and Minamoto 2018; Minamoto et al. 2013) revealed that the convolution of the reaction zones in MILD combustion is strongly influenced by the dilution factor of the reacting mixture as well as the turbulence–chemistry interaction. This confirms the experimental finding reported by Dally, Riesmeier, and Peters (2004). A review of the physical insights gained from the previous DNS studies was compiled by Swaminathan (2019) and noted that homogeneous and inhomogeneous mixture MILD combustion for similar thermochemical and turbulence parameters behave similarly in terms of temperature rise and heat release rate.

Recent studies using DNS data indicated that flamelet-based models can be applied to gaseous fuel MILD combustion with some modifications (Minamoto and Swaminathan 2015; Minamoto et al. 2013). However, when applied to MILD combustion in Reynolds averaged Navier–Stokes (RANS) and large eddy simulations (LES) type investigations, the conventional flamelet, and eddy dissipation approaches provided good agreement with experimental data in terms of mean velocity and temperature fields, but discrepancies have been reported in terms of peak temperature and minor species concentrations (Aminian et al. 2011; Christo and Dally 2005). Thus, further investigations are needed to fully understand the behavior of flamelet-based approaches in MILD combustion. In this context, the natural first step is to investigate the statistical behavior of the magnitude of the gradient of the reaction progress variable (Surface Density Function (SDF) =$|\mathrm{\nabla }c|$ where $c$ is the reaction progress variable) since it is closely related to the Scalar Dissipation Rate (SDR) ($N_c=D_c|\mathrm{\nabla }c{|}^2$ with $D_c$ being the reaction progress variable diffusivity) and generalized Flame Surface Density (FSD) (${\mathrm{\Sigma }}_{gen}=\overline{|\mathrm{\nabla }c|}$ with the overbar suggesting a Reynolds averaging/LES filtering operation as appropriate). Particularly, it is important to understand the influence of the strain rates induced by fluid motion and flame propagation on the SDF evolution in homogeneous and inhomogeneous mixture MILD combustion. The present analysis aims to address this gap in the literature by analyzing data from three-dimensional DNS with reduced n-heptane/air chemical mechanism for EGR type homogeneous and inhomogeneous mixture MILD combustion and at two dilution levels (i.e. mole fraction of oxygen in the unburned gas $X_{O_2}=3.0\mathrm{\%},4.5\mathrm{\%}$). Previous work by the current authors and coworkers (Awad et al. 2021) focused on the effects of turbulence intensity on the statistical behavior of the SDF in homogeneous MILD combustion of diluted methane–air mixtures. Awad et al. (2021) compared the SDF statistics in homogeneous mixture MILD combustion with the corresponding statistics in conventional turbulent premixed flames and found that the mean value of SDF is significantly reduced in the turbulent MILD cases compared to the corresponding unstretched diluted laminar flames, while in the turbulent premixed cases, the mean value of SDF was comparable to the corresponding unstretched 1D laminar flame. Moreover, the low heat release in the MILD cases resulted in a negligible dilatation rate effect with $\mathrm{\nabla }c$ preferentially aligned with the most compressive principal strain rate eigndirection, whereas in the turbulent premixed cases, $\mathrm{\nabla }c$ showed greater propensity to align with the most extensive strain rate eigndirection due to the significant dilatation rate effect caused by the strong heat release. Furthermore, the mean displacement speed in the turbulent premixed cases was found to increase from the unburned gas side toward the burned gas side, while in the MILD cases, the opposite trend was observed. Finally, in the MILD cases, the contribution of the reaction to the displacement speed was found to take a non-negligible value across a wider range of $c$ values compared to the conventional premixed cases where the reaction contribution was significant only toward the burned gas side. Awad et al. (2021) also showed that the difference in the mean behavior of the displacement speed between the MILD and conventional cases has significant implications on the mean behavior of the normal and tangential strain rates induced by flame propagation. The present work extends the previous analysis by Awad et al. (2021) by focusing on the effects of dilution level and mixture inhomogeneity on the statistical behavior of the SDF. Furthermore, the current study looks at the combustion of n-heptane which has a radically different oxidization path compared to methane, but is more relevant to the fuel mixtures used in industry.

The main objectives of the current analysis are:

1. A comparison between the SDF statistics in EGR-type homogeneous and inhomogeneous mixture MILD combustion.

2. Highlighting the effect of oxygen dilution levels on SDF statistics.

The rest of the paper will be organized as follows. The mathematical background of this analysis will be presented in the next section. This will be followed by a brief discussion of the numerical implementation. Following that, results will be presented and discussed.

Mathematical background

The reaction progress variable ($c$) is used to quantify the extent of reaction completion. In this work, $c$ has been defined based on the oxidizer mass fraction ($Y_{O_2}$) as follows:

(1) $c=\frac{(1-\xi )Y_{O_2,2}-Y_{O_2}}{(1-\xi )Y_{O_2,2}-max\left[0,\frac{{\xi }_{st}-\xi }{{\xi }_{st}}{Y}_{O_2,2}\right]}$(1)

where $Y_{O_2,2}$ is the oxidizer mass fraction in the oxidizer stream, while $\xi$ and ${\xi }_{st}$ are the mixture fraction and the stoichiometric mixture fraction, respectively.

The choice of oxygen as the progress variable species is justified since larger hydrocarbons, like n-heptane, tend to break down to smaller ones at high temperatures. Thus, defining the progress variable based on fuel will be less representative of the whole combustion process.

Following Bilger (1989), both $\xi$ and ${\xi }_{st}$ were defined in terms of the elemental mass fractions and atomic masses as:

(2) $\xi =\frac{\beta -{\beta }_2}{{\beta }_1-{\beta }_2}$(2)

and

(3) ${\xi }_{st}=\frac{-{\beta }_2}{{\beta }_1-{\beta }_2}$(3)

where $\beta =2Z_c/W_c+0.5Z_H/W_H-Z_O/W_O$, $Z_j$ and $W_j$ are the elemental mass fractions and atomic masses for carbon, oxygen, and hydrogen atoms. The subscripts 1, 2 refer to the fuel and (diluted) oxidizer streams, respectively. The equivalence ratio $\varphi$ as a function of the mixture fraction is given by:

(4) $\varphi =\frac{\xi (1-{\xi }_{st})}{{\xi }_{st}(1-\xi )}$(4)

A transport equation for the reaction progress variable can be written as:

(5) $\rho \frac{\mathrm{\partial }c}{\mathrm{\partial }t}+\rho u_j\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\rho D_c\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}\right)+{\dot{w}}_c+A_{\xi }$(5)

where $\rho$ is the gas density, $D_c$ is the mass diffusivity of the species on which the progress variable is based, $u_j$ is the $j^{th}$ velocity component, ${\dot{w}}_c$ is the reaction rate of the progress variable and $A_{\xi }$ arises due to $\xi$ gradient. Both ${\dot{w}}_c$ and $A_{\xi }$ are expressed as functions of $\xi$ and $Y_{O_2,2}$ as follows:

(6) ${\dot{w}}_c=\left\{\begin{array}{cc}\frac{-{\xi }_{st}{\dot{w}}_{O_2}}{Y_{O_2,2}\xi (1-{\xi }_{st})}& \mathrm{f}\mathrm{o}\mathrm{r}\xi \le {\xi }_{st}\\ \frac{-{\dot{w}}_{O_2}}{Y_{O_2,2}(1-\xi )}& \mathrm{f}\mathrm{o}\mathrm{r}\xi >{\xi }_{st}\end{array}\right.$(6)

(7) $A_{\xi }=\left\{\begin{array}{cc}\frac{2\rho D_c\mathrm{\nabla }\xi .\mathrm{\nabla }c}{\xi }& \mathrm{f}\mathrm{o}\mathrm{r}\xi \le {\xi }_{st}\\ \frac{-2\rho D_c\mathrm{\nabla }\xi .\mathrm{\nabla }c}{(1-\xi )}& \mathrm{f}\mathrm{o}\mathrm{r}\xi >{\xi }_{st}\end{array}\right.$(7)

For a given $c$ iso-surface, eq. (5)(5) $\rho \frac{\mathrm{\partial }c}{\mathrm{\partial }t}+\rho u_j\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\rho D_c\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}\right)+{\dot{w}}_c+A_{\xi }$(5) can be rewritten as (Malkeson and Chakraborty 2010):

(8) $\frac{\mathrm{\partial }c}{\mathrm{\partial }t}+u_j\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}=S_d|\mathrm{\nabla }c|$(8)

where $S_d$ is the displacement speed of the $c$ iso-surface and is given by (Echekki and Chen 1996):

(9) $S_d=\frac{{\dot{w}}_c+\mathrm{\nabla }.(\rho D_c\mathrm{\nabla }c)+A_{\xi }}{\rho |\mathrm{\nabla }c|}$(9)

The displacement speed can be written in terms of its components as follows:

(10) $S_d=S_r+S_n+S_t+S_{\xi }$(10)

where $S_r,S_n,S_t\mathrm{a}\mathrm{n}\mathrm{d}{S}_{\xi }$ are the reaction, normal and tangential diffusion, and mixture inhomogeneity components of $S_d$, respectively. These components are defined as (Echekki and Chen 1999):

(11) $\begin{array}{cc}{S}_r& =\frac{{\dot{w}}_c}{\rho |\mathrm{\nabla }c|}\\ \\ S_n& =\frac{\vec{N}.\mathrm{\nabla }(\rho D_c\vec{N}.\mathrm{\nabla }c)}{\rho |\mathrm{\nabla }c|}\\ \\ S_t& =-2D_c{\kappa }_m\\ \\ S_{\xi }& =\frac{A_{\xi }}{\rho |\mathrm{\nabla }c|}\end{array}$(11)

where $\vec{N}=-\mathrm{\nabla }c/|\mathrm{\nabla }c|$ is the flame normal vector and ${\kappa }_m=0.5\mathrm{\nabla }.\vec{N}$ is the flame curvature. It is worth noting that while eqs. (5)(5) $\rho \frac{\mathrm{\partial }c}{\mathrm{\partial }t}+\rho u_j\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\rho D_c\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}\right)+{\dot{w}}_c+A_{\xi }$(5) to (11(11) $\begin{array}{cc}{S}_r& =\frac{{\dot{w}}_c}{\rho |\mathrm{\nabla }c|}\\ \\ S_n& =\frac{\vec{N}.\mathrm{\nabla }(\rho D_c\vec{N}.\mathrm{\nabla }c)}{\rho |\mathrm{\nabla }c|}\\ \\ S_t& =-2D_c{\kappa }_m\\ \\ S_{\xi }& =\frac{A_{\xi }}{\rho |\mathrm{\nabla }c|}\end{array}$(11) ) were initially developed for conventional flames propagating into homogeneous and stratified mixtures (i.e. arrangements dominated solely by flame propagation), their use remains valid in MILD combustion cases (which include both flame propagation and auto-ignition events) since it was demonstrated by Minamoto et al. (2014a) that MILD combustion exhibits significant amount of propagation-dominated behavior. Moreover, the concept of displacement speed was used previously in ignition-dominated problems (Chen et al. 2006; Desai, Sankaran, and Im 2020). It will be shown later in the paper that the cases considered here also exhibit a significant amount of propagation-dominated behavior which justifies the utilization of methodologies usually encountered in premixed combustion modeling for this analysis. The kinematic form of the reaction progress variable equation (eq. (8))(8) $\frac{\mathrm{\partial }c}{\mathrm{\partial }t}+u_j\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}=S_d|\mathrm{\nabla }c|$(8) can be used to obtain a transport equation for the Surface Density Function (SDF = $|\mathrm{\nabla }c|$) at a given $c$ iso-surface which can be written as (Chakraborty and Cant 2005; Dopazo et al. 2015; Sankaran et al. 2007):

(12) $\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }t}+\frac{\mathrm{\partial }(u_j|\mathrm{\nabla }c|)}{\mathrm{\partial }{x}_j}=a_T|\mathrm{\nabla }c|-\frac{\mathrm{\partial }\left(S_d{N}_j|\mathrm{\nabla }c|\right)}{\mathrm{\partial }{x}_j}+2S_d{\kappa }_m|\mathrm{\nabla }c|$(12)

where $a_T=({\delta }_{ij}-N_i{N}_j)\mathrm{\partial }{u}_i/\mathrm{\partial }{x}_j$ is the flame tangential strain rate due to fluid motion, $-\mathrm{\partial }\left(S_d{N}_j|\mathrm{\nabla }c|\right)/\mathrm{\partial }{x}_j$ is the SDF propagation term and $2S_d{\kappa }_m|\mathrm{\nabla }c|$ is the SDF curvature term. Equation (12)(12) $\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }t}+\frac{\mathrm{\partial }(u_j|\mathrm{\nabla }c|)}{\mathrm{\partial }{x}_j}=a_T|\mathrm{\nabla }c|-\frac{\mathrm{\partial }\left(S_d{N}_j|\mathrm{\nabla }c|\right)}{\mathrm{\partial }{x}_j}+2S_d{\kappa }_m|\mathrm{\nabla }c|$(12) can be rewritten in a reference frame attached to the flame as follows (Dopazo et al. 2015):

(13) $\frac{1}{|\mathrm{\nabla }c|}\frac{d|\mathrm{\nabla }c|}{dt}=-\frac{1}{\mathrm{\Delta }{x}_N}\frac{d\mathrm{\Delta }{x}_N}{dt}=-a_N-N_j\frac{\mathrm{\partial }{S}_d}{\mathrm{\partial }{x}_j}=-a_N^{eff}$(13)

where $d(\dots )/dt=\mathrm{\partial }(\dots )/\mathrm{\partial }t+V_j^c\mathrm{\partial }(\dots )/\mathrm{\partial }{x}_j$ is the total derivative in the reference frame attached to the flame, $V_j^c=u_j+S_d{N}_j$ is the $j^{th}$ component of the flame propagation velocity, $\mathrm{\Delta }{x}_N$ is the distance between two adjacent $c$ iso-surfaces. The effective flame normal strain rate $a_N^{eff}$ has contributions from the fluid-dynamic normal strain rate $a_N$ and the flame normal strain rate due to flame propagation $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$.

It is also worthwhile to investigate the evolution of the elemental reaction zone surface area ($\delta A$) since it is closely related to $|\mathrm{\nabla }c|$ (i.e. $A=\underset{V}{\int }\left|\mathrm{\nabla }c\right|dV$). The transport equation for $\delta A$ is given by (Candel and Poinsot 1990; Dopazo et al. 2015; Pope 1988):

(14) $\frac{1}{\delta A}\frac{d(\delta A)}{dt}=a_T+2S_d{\kappa }_m=a_T^{eff}$(14)

where $a_T^{eff}$ is the effective flame tangential strain rate, $a_T$ is the fluid-dynamic tangential strain rate, and $2S_d{\kappa }_m$ arises due to flame propagation.

From the previous presentation, it is evident that $a_N$, $a_T$, $2S_d{\kappa }_m$, $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$, $a_N^{eff}$ and $a_T^{eff}$ drive the statistical behaviors of the SDF and flame surface area. Thus, the statistical behaviors of these strain rates will be discussed in this work.

Finally, since MILD combustion includes regions of auto-ignition, interacting flames and propagating flames, it can be useful to characterize these regions based on the balance of convection, diffusion, and chemical reaction processes. This balance can be assessed using (Minamoto et al. 2014a):

(15) $\mathcal{B}\equiv |\mathcal{C}-\mathcal{D}|-|\mathcal{R}|$(15)

where $\mathcal{C},\mathcal{D}$ and $\mathcal{R}$ are taken from eq. (5)(5) $\rho \frac{\mathrm{\partial }c}{\mathrm{\partial }t}+\rho u_j\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\rho D_c\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}\right)+{\dot{w}}_c+A_{\xi }$(5) and written as (Minamoto et al. 2014a):

(16) $\begin{array}{cc}\mathcal{C}& =\frac{\mathrm{\partial }\rho u_{j}c}{\mathrm{\partial }{x}_j}\\ \mathcal{D}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\rho D_c\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}\right)\\ \mathcal{R}& ={\dot{w}}_c\end{array}$(16)

Following Minamoto et al. (2014a), regions dominated by autoignition and interacting flames (Reaction) will give rise to $\mathcal{B}<0$, while propagating-flame regions will yield $\mathcal{B}>0$. Conditioning SDF, $S_d$, $2S_d{\kappa }_m$ and $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$ on $\mathcal{B}$ value will give insight into the contribution of each combustion mode in driving the SDF behavior.

Numerical implementation

For the current study, the compressible DNS code SENGA+ (Cant 2012) has been used. In this code, the standard conservation equations of mass, momentum, energy, and species mass fractions for compressible reacting flows are solved in dimensional form with the thermo-physical properties such as viscosity, thermal conductivity, specific heat capacities, and mass diffusivities taken to be temperature dependent. The SENGA+ code employs a ${10}^{th}$ order central difference scheme for spatial differentiation at the internal grid points, but the order of accuracy gradually reduces to a one-sided fourth-order scheme at the non-periodic boundary (Cant 2012). A fourth-order explicit low storage Runge-Kutta scheme is used for time advancement. The boundary conditions are specified according to the Navier-Stokes Characteristic Boundary Condition (NSCBC) methodology. A reduced chemical mechanism comprising 22 species and 18-steps (Liu et al. 2004) has been taken to represent the chemical kinetics of combustion.

A cubic domain with an edge $L=20mm$ is discretized by a Cartesian grid with uniform spacing comprising 216 nodes in each direction. The resulting grid spacing ensures that both the thermal flame thickness (${\delta }_{th}=(T_p-T_r)/max|\mathrm{\nabla }T{|}_L$ where $T,T_p$ and $T_r$ are the instantaneous, products and reactants’ temperatures and the subscript $L$ refers to the 1D unstretched laminar flame), and the Kolmogorov length scale ($\eta$) are adequately resolved (e.g. a minimum of 15 and 1.1 grid points in ${\delta }_{th}$ and $\eta$, respectively). Turbulent inflow boundary with specified density, velocity, and species have been imposed at the left-hand $x_1$ boundary and a partially non-reflective outflow boundary condition had been specified at the right-hand $x_1$ boundary. Periodic boundary conditions have been imposed on all other boundaries. To illustrate the effects of dilution and mixture inhomogeneity, two homogeneous mixture cases at different dilution levels have been considered in this work (HM-A with ${X_O}_{{ }_2}=4.5\mathrm{\%}$ and HM-B with ${X_O}_{{ }_2}=4.5\mathrm{\%}$), as well as one inhomogeneous mixture case (IM-A with ${X_O}_{{ }_2}=4.5\mathrm{\%}$). The current simulations have been set up following the methodology detailed by Minamoto et al. (2013) for the homogeneous cases and Doan, Swaminathan, and Minamoto (2018) for the inhomogeneous case. This methodology consists of the following steps:

1. An incompressible homogeneous isotropic field is used to initialize the turbulent velocity fluctuations following the Batchelor-Townsend spectrum (Batchelor and Townsend 1948). This is done using the pseudo-spectral method of Rogallo (1981). For all cases, the initial turbulent conditions chosen for this study are $u{ }_0^{\prime }=2.0$ m/s and ${\mathrm{\ell }}_0=0.002$ m. These values are comparable to those reported by Oldenhof et al. (2011) for their experimental analysis.

2. For each homogeneous case, a corresponding one-dimensional unstretched freely propagating laminar premixed flame simulation was carried out using the thermochemical conditions shown in Table 1 and $\varphi =0.8$. The species mass fractions resulting from the laminar simulations were then written as functions of the progress variable based on oxygen mass fraction (i.e. $Y_{\alpha ,L}=f_{\alpha }^{HM}(c_{O2,L})$). For the inhomogeneous case, a number of laminar flame simulations were conducted covering the range of equivalence ratio given by $\varphi =0.3-1.3$ with a mean of $0.8$. The resulting species mass fractions were parametrized as functions of both oxygen mass fraction-based reaction progress variable and equivalence ratio (i.e. $Y_{\alpha ,L}=f_{\alpha }^{IM}(c_{O2,L},\varphi )$). The unstretched laminar burning velocity and Zeldovich flame thickness (${\delta }_f={\alpha }_T/S_L$ where ${\alpha }_T$ is the thermal diffusivity of the mixture) of the above laminar simulations are also shown in Table 1.

3. For each case (both homogeneous and inhomogeneous), an initial bimodal distribution of $c$ with length scale ${\mathrm{\ell }}_c$ (given in Table 1) was generated using the methodology by Eswaran and Pope (1988). For the inhomogeneous case, an additional bimodal distribution of the equivalence ratio with peaks at $\varphi =0.3$ and $1.3$ was also generated with the integral length scale ${\mathrm{\ell }}_{\varphi }$ given in Table 1. Figure 1 shows the probability density function of the equivalence ratio for the initial field of the inhomogeneous case.

   Figure 1. PDF of the equivalence ratio ϕ in the initial scalar field for case IM-A.

   

4. The functions $Y_{\alpha ,L}$ formulated in step 2 are then used to create the initial scalar field using the bimodal $c$ field for the homogeneous cases and the bimodal $c$ and $\varphi$ fields in the inhomogeneous case, generated in step 3, as an input parameter. Here, atmospheric pressure and an unburned gas temperature $T_r=1100$ K were also used. The chosen unburned gas temperature is comparable to that used in the experimental investigation by Ye et al. (2017).

5. The resulting scalar field was allowed to evolve with turbulence without reaction for about one turnover time in a periodic domain mimicking the mixing of fuel and oxidizer in an EGR-type combustor. The mixing time is well below the ignition delay time. For example, at the composition of case HM-A (Table 1) and $T_r=1100$ K, the ignition delay time, ${\tau }_{ign}$, estimated using a 0D ideal gas reactor simulation in Cantera (Goodwin et al. 2022) , is found to be ${\tau }_{ign}=37.83ms$ while the mixing time ${\tau }_{turb}={\mathrm{\ell }}_{{ }_0}/u{ }_0^{\prime }=1.0ms$.

6. Finally, the evolved scalar field was then used as an initial field for the reacting simulation as well as being fed into the domain through the inlet with a mean velocity $U_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}=6.0$ m/s.

Table 1. DNS initial conditions for MILD combustion. Values with * were calculated at the mean equivalence ratio of 0.8.

Case	$X_{O_2,2}$	$X_{C{O}_2,2}$	$X_{H_{2}O,2}$	ϕ	$S_L$ [m/s]	${\delta }_f$ [m]	${\mathrm{\ell }}_c/{\mathrm{\ell }}_0$	$\mathrm{\ell }$ϕ$/{\mathrm{\ell }}_c$
HM-A	0.045	0.097	0.112	0.8	0.420	4.55$\times {10}^{-4}$	1.225	-
HM-B	0.030	0.106	0.122	0.8	0.246	7.84$\times {10}^{-4}$	1.225	-
IM-A	0.045	0.097	0.112	0.3–1.3	0.420*	4.55$\times {10}^{-4\ast }$	1.225	1.34

The simulations have been run for 2.0 through-pass times (i.e. ${\tau }_{\mathrm{s}\mathrm{i}\mathrm{m}}=2.0L/U_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$), which amounts to 6.67 initial eddy turnover times.

Results and discussion

A view of the instantaneous field

The instantaneous views of $c$ iso-surfaces and the normalized instantaneous temperature field for all cases are shown in Figure 2. The iso-surfaces are taken at $c$ values corresponding to the maximum heat release for the unstretched laminar premixed flames used to initialize each case (shown in Table 1), and thus it can reasonably be considered to represent the flame surface. The $c$ iso-surfaces in Figure 2 show the distributed nature of the flame in MILD combustion as observed in OH-PLIF visualizations by Plessing, Peters, and Wünning (1998). Moreover, it can be seen from Figure 2 that a considerable amount of flame self-interactions is apparent from the shown $c$ iso-surface. The nondimensional temperature (i.e. $c_T=(T-T_r)/(T_p-T_r)$) fields shown in Figure 2 indicate that only modest variations in temperature are reported through the domain with no sharp gradients in temperature visible in all cases. This is consistent with expectations in MILD combustion as shown by Minamoto et al. (2013).

Figure 2. The iso-surfaces at $c$ values corresponding to the maximum heat release location in the unstretched laminar premixed flame used to initialise each case and the instantaneous normalised temperature $c_T$ field for all cases considered here.

Mean behavior of the SDF

The profiles of the mean values of the normalized SDF conditioned upon $c$ are shown in Figure 3. Here ${\delta }_{\mathrm{\ell }}=1/max(|\mathrm{\nabla }c{|}_L)$ is the inverse of the maximum value of SDF in the corresponding unstretched laminar flame. Since the flamelet thickness can be considered to be proportional to the inverse of the maximum of $|\mathrm{\nabla }c|\times {\delta }_{\mathrm{\ell }}$ conditioned on $c$, the lower values of $⟨|\mathrm{\nabla }c|{⟩}_c$ (where $⟨\dots {⟩}_c$ is the mean value conditioned upon $c$) are indicative of thickened reaction zones in the turbulent cases. A comparison between Figures 3a and 3c shows that increasing the dilution level leads to a reduced flame thickening (relative to the corresponding unstretched laminar flame) under turbulent conditions compared to the low-dilution case. It is worth noting that increasing dilution gives rise to thicker flames, in general, both under laminar and turbulent conditions, however, Figure 3 shows that under turbulent conditions the more diluted flame does not exhibit as much thickening in comparison to the corresponding ${\delta }_{\mathrm{\ell }}$, as in the low dilution case. In the inhomogeneous case, a slightly thicker flame toward the unburned gas side can be observed, but in general the level of flame thickening is comparable to that in the homogeneous case with the same dilution level (HM-A).

Figure 3. Profiles of the mean values of $|\mathrm{\nabla }c|\times {\delta }_{\mathrm{\ell }}$ conditioned upon $c$ for all cases considered here. Here ’LAM’ curves refers to the diluted laminar flame used to create the initial MILD field following the methodology by Minamoto et al. (2013).

Fluid-dynamic strain rates

Figure 4 shows the profiles of the mean values of dilatation rate, normal strain rate, and tangential strain rate due to fluid flow conditioned upon $c$ for all the cases considered here. It can be seen from Figure 4 that the mean value of dilatation rate $\mathrm{\nabla }\cdot \vec{u}$ conditioned upon $c$ is negligible in all cases. This is expected since the heat release in MILD combustion is characteristically low, and thus leads to weak thermal expansion effects. The low heat release also leads to a small change in temperature, as reported in MILD combustion experiments (Dally, Karpetis, and Barlow 2002; Dally, Riesmeier, and Peters 2004; Plessing, Peters, and Wünning 1998) and shown in Figure 2. Figure 4 shows that the mean value of $a_N$ conditioned upon $c$ remains negative across the flame front. In order to explain this behavior, $a_N$ can be expressed in the following manner (Awad et al. 2021; Dopazo et al. 2015):

(17) $a_N=(e_{\alpha }co{s}^2{\theta }_{\alpha }+e_{\beta }co{s}^2{\theta }_{\beta }+e_{\gamma }co{s}^2{\theta }_{\gamma })$(17)

Figure 4. Profiles of the normalised mean values of dilatation rate $\mathrm{\nabla }\cdot \vec{u}$, normal strain rate $a_N$ and tangential strain rate $a_T$ due to background fluid motion conditioned upon $c$ for all cases considered here.

where $e_{\alpha }$, $e_{\beta }$ and $e_{\gamma }$ are the most extensive, intermediate, and most compressive principal strain rates, respectively, and ${\theta }_{\alpha }$, ${\theta }_{\beta }$ , and ${\theta }_{\gamma }$ are angles between $\mathrm{\nabla }c$ and eigenvectors associated with $e_{\alpha }$, $e_{\beta }$ and $e_{\gamma }$, respectively. Awad et al. (2021) showed in a similar configuration to the present study, but using methane as the fuel, that $\mathrm{\nabla }c$ aligns preferentially with the most compressive principal strain rate eigendirection similar to passive scalar mixing and thus gives rise to a predominantly negative normal strain rate $a_N$ because of $|cos{\theta }_{\gamma }|\approx 1.0$ and $e_{\gamma }<0$. A similar qualitative behavior has been observed for all the cases considered here and thus are not shown explicitly for the sake of conciseness. The negative mean values of $a_N$ lead to positive mean values of tangential strain rate $a_T=\mathrm{\nabla }\cdot \vec{u}-a_N$ due to the negligible mean values of $\mathrm{\nabla }\cdot \vec{u}$. The predominance of positive values of $a_T$ is indicative of the flame area generation due to fluid-dynamic straining in all cases. Figure 4 shows that mixture inhomogeneity or changing the dilution level does not affect the mean values of fluid-dynamic strain rates.

Mean behavior of displacement speed and its components

Figure 5 shows the variations of the mean values of the normalized displacement speed $S_d$ and its components conditioned upon $c$. It can be seen from Figure 5 that, in all cases considered, the mean displacement speed has the same order of magnitude and remains positive toward the unburned gas side but with its value reducing with increasing $c$ to cross zero at about $c\approx 0.6$ which suggests that the leading edge of the flame front propagates into the reactants, whereas the trailing edge retreats into the products indicating local thickening in the reaction zone in the mean sense. For the high dilution case (HM-B), the mean displacement speed continues to decrease monotonically with $c$. However, for the low dilution cases (HM-A and IM-A), the mean value of $S_d$ starts to increase with $c$ after reaching a minimum. This is due to the significant increase in the mean value of the reaction component of displacement speed $S_r$ in these cases toward the burned gas side.

Figure 5. Profiles of the normalised mean values of $S_d$ and its components $S_r$, $S_n$, $S_t$ and $S_{\xi }$ conditioned upon $c$ for all cases considered here.

To fully understand the behavior of $S_d$ components, the mean contributions of ${\dot{w}}_c$, $\mathrm{\nabla }\cdot (\rho D_c\mathrm{\nabla }c)$ and $A_{\xi }$ (where present) conditioned upon $c$ are shown in Figure 6. Figure 6 shows that the mean reaction rate $⟨{\dot{w}}_c{⟩}_c$ remains small (but not zero) for small $c$ values, but its magnitude increases with $c$ and the peak value is obtained toward the burned gas side in all cases. However, it can be seen from Figure 6 that increasing the dilution level leads to reduced values of $⟨{\dot{w}}_c{⟩}_c$. The reduced reactivity with the increased dilution level is also evident when comparing ${\dot{w}}_c$ in the laminar flames corresponding to cases HM-A (low dilution, $X_{O_2}=4.5\mathrm{\%}$) and HM-B (high dilution, $X_{O_2}=3.0\mathrm{\%}$) as shown in Figure 7. Figure 6 also shows that, in the inhomogeneous mixture case (IM-A), an increased peak value for $⟨{\dot{w}}_c{⟩}_c$ can be observed in comparison to the homogeneous case (HM-A). This increase in $⟨{\dot{w}}_c{⟩}_c$ with $c$ drives the increase in $⟨S_r{⟩}_c$ values, which in turn leads to increased $⟨S_d{⟩}_c$ values toward the burned gas side. Figure 8 shows the probability density function of the equivalence ratio for the inhomogeneous case (IM-A) at different $c$ values. It can be seen from Figure 8 that, for case IM-A, the peak in $⟨{\dot{w}}_c{⟩}_c$ observed in Figure 6 occurs in mostly lean mixtures and that the equivalence ratio in the inhomogeneous case (IM-A) generally evolves toward lean mixtures as the reaction progresses.

Figure 6. Profiles of the normalised mean values of ${\dot{w}}_c$, $\mathrm{\nabla }\cdot (\rho D_c\mathrm{\nabla }c)$ and $A_{\xi }$ (where present) conditioned upon $c$ for all cases considered here.

Figure 7. Profiles of the normalised ${\dot{w}}_c$ in the laminar flames upon which the cases HM-A and HM-B are based.

Figure 8. Pdfs of the equivalence ratio ϕ in the inhomogeneous case (IM-A) at different values for the progress variable $c$.

Figure 6b shows that the mean contribution of the mixture inhomogeneity component $⟨A_{\xi }{⟩}_c$ remains negligible for all $c$ values which is consistent with previous findings for turbulent-stratified flames (Malkeson and Chakraborty 2010). This further suggests that the effects of mixture inhomogeneity are predominantly felt through the contributions of ${\dot{w}}_c$ and $\mathrm{\nabla }\cdot (\rho D_c\mathrm{\nabla }c)$. Both $⟨\vec{N}\cdot \mathrm{\nabla }(\rho D_c\vec{N}\cdot \mathrm{\nabla }c){⟩}_c$ and $⟨-2\rho D_c{\kappa }_m|\mathrm{\nabla }c|{⟩}_c$ have approximately the same magnitude and exhibit positive values toward the unburned gas side but assume negative values with magnitudes comparable to that of $⟨{\dot{w}}_c{⟩}_c$ toward the burned gas side in all cases. Consequently, $⟨S_t{⟩}_c$ and $⟨S_n{⟩}_c$ follow the qualitative trends of $⟨\vec{N}\cdot \mathrm{\nabla }(\rho D_c\vec{N}\cdot \mathrm{\nabla }c){⟩}_c$ and $⟨-2\rho D_c{\kappa }_m|\mathrm{\nabla }c|{⟩}_c$, respectively, and remain the dominant components of $⟨S_d{⟩}_c$ at $c⪅0.6$. This leads to the positive values of $⟨S_d{⟩}_c$ toward the unburned gas side. Figures 5 and 6 show that the dilution level and mixture inhomogeneity have limited impacts on $⟨\vec{N}\cdot \mathrm{\nabla }(\rho D_c\vec{N}\cdot \mathrm{\nabla }c){⟩}_c$ and $⟨-2\rho D_c{\kappa }_m|\mathrm{\nabla }c|{⟩}_c$ and consequently on the diffusion components of $⟨S_d{⟩}_c$ (i.e. $⟨S_n{⟩}_c$ and $⟨S_t{⟩}_c$).

Strain rates due to flame propagation

Normal strain rate

The profiles of the normalized mean values of the additional flame normal strain rate due to flame propagation $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$ and its components (i.e. $N_j\mathrm{\partial }{S}_t/\mathrm{\partial }{x}_j$, $N_j\mathrm{\partial }{S}_n/\mathrm{\partial }{x}_j$, $N_j\mathrm{\partial }{S}_r/\mathrm{\partial }{x}_j$ and $N_j\mathrm{\partial }{S}_{\xi }/\mathrm{\partial }{x}_j$) conditioned upon $c$ are shown in Figure 9. It can be seen from Figure 9 that $⟨N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j{⟩}_c$ exhibits positive values across the whole range of $c$ values with comparable magnitudes in all cases. Both $⟨N_j\mathrm{\partial }{S}_n/\mathrm{\partial }{x}_j{⟩}_c$ and $⟨N_j\mathrm{\partial }{S}_t/\mathrm{\partial }{x}_j{⟩}_c$ have positive values and the same order of magnitudes across the flame front with these terms being the dominant contributors to $⟨N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j{⟩}_c$ up to $c\approx 0.7$. The strengthening of $N_j\mathrm{\partial }{S}_t/\mathrm{\partial }{x}_j$ contributions in MILD combustion can be attributed to the importance of curvature effects for the cases with high $Ka$ (i.e. $Ka\gg 1$) (Peters 2000). Figure 9 shows that $⟨N_j\mathrm{\partial }{S}_r/\mathrm{\partial }{x}_j{⟩}_c$ has a negligible value up to $c\approx 0.7$ but then becomes negative and its contribution to $⟨N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j{⟩}_c$ becomes significant. The changes in $⟨N_j\mathrm{\partial }{S}_r/\mathrm{\partial }{x}_j{⟩}_c$ value at $c⪆0.7$ occur due to the increases in $⟨{\dot{w}}_c{⟩}_c$, shown in Figure 6. It can be seen from Figure 9b that the mixture inhomogeneity component $⟨N_j\mathrm{\partial }{S}_{\xi }/\mathrm{\partial }{x}_j{⟩}_c$ remains negligible across the flame front. Similarly, the effect of increased dilution remains mostly marginal apart from the decrease in $⟨N_j\mathrm{\partial }{S}_r/\mathrm{\partial }{x}_j{⟩}_c$ values at $c⪆0.7$ due to the decrease in $⟨{\dot{w}}_c{⟩}_c$ in the high dilution case, but this has a marginal effect on $⟨N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j{⟩}_c$.

Figure 9. Profiles of the normalised mean values of normal strain rate due to flame propagation $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$ and its components (i.e. $N_j\mathrm{\partial }{S}_r/\mathrm{\partial }{x}_j$, $N_j\mathrm{\partial }{S}_n/\mathrm{\partial }{x}_j$, $N_j\mathrm{\partial }{S}_t/\mathrm{\partial }{x}_j$ and $N_j\mathrm{\partial }{S}_{\xi }/\mathrm{\partial }{x}_j$ where appropriate) conditioned upon $c$ for all cases considered here.

Tangential strain rate

Figure 10 shows the normalized mean values of the tangential flame strain rate due to flame propagation $2S_d{\kappa }_m$ and its components (i.e. $2S_t{\kappa }_m$, $2S_n{\kappa }_m$, $2S_r{\kappa }_m$ and $2S_{\xi }{\kappa }_m$) conditioned on the progress variable $c$. It can be seen from Figure 10 that $⟨2S_d{\kappa }_m{⟩}_c$ has a negative value across the flame front with the negative definite value of $⟨2S_t{\kappa }_m{⟩}_c=-4⟨D_c{\kappa }_m^2{⟩}_c$ being the major contributor to $⟨2S_d{\kappa }_m{⟩}_c$. Moreover, $⟨2S_n{\kappa }_m{⟩}_c$ assumes negative values and its magnitude does not change appreciably across the range of $c$ values, while $⟨2S_r{\kappa }_m{⟩}_c$ has a negligible value up to $c\approx 0.7$ where it starts to increase and assumes positive values. Similarly to $⟨N_j\mathrm{\partial }{S}_r/\mathrm{\partial }{x}_j{⟩}_c$, the increase in $⟨2S_r{\kappa }_m{⟩}_c$ for $c>0.7$ corresponds to the increase in $⟨{\dot{w}}_c{⟩}_c$ value shown in Figure 6 and shows some dependence on the dilution level since $⟨{\dot{w}}_c{⟩}_c$ is dependant on the dilution level as well. Figure 10b shows that the effect of the mean value of the curvature stretch due to mixture inhomogeneity ($⟨2S_{\xi }{\kappa }_m{⟩}_c$) remains negligible in the range considered in this study.

Figure 10. Profiles of the normalised mean values of the tangential strain rate $2S_d{\kappa }_m$ due to flame propagation and its components (i.e. $2S_r{\kappa }_m$, $2S_n{\kappa }_m$, $2S_t{\kappa }_m$ and $2S_{\xi }{\kappa }_m$, where appropriate) conditioned upon $c$ for all cases considered here.

Effective strain rates

The profiles of the normalized mean values of the normal ($a_N^{eff}$) and tangential ($a_T^{eff}$) effective strain rates conditioned upon $c$ are shown in Figure 11. A comparison between Figures 4, 9 and 10 shows that the mean values of the fluid-dynamic strain rates, both $⟨a_N{⟩}_c$ and $⟨a_T{⟩}_c$, are almost an order of magnitude smaller compared to those arising from flame propagation (i.e. $⟨N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j{⟩}_c$ and $⟨2S_d{\kappa }_m{⟩}_c$). Thus, the mean values of $a_N^{eff}$ and $a_T^{eff}$ are predominantly influenced by the flame propagation-induced strain rates. For all cases, $⟨a_N^{eff}{⟩}_c$ maintains positive values for all $c$ levels. Following eq. (13)(13) $\frac{1}{|\mathrm{\nabla }c|}\frac{d|\mathrm{\nabla }c|}{dt}=-\frac{1}{\mathrm{\Delta }{x}_N}\frac{d\mathrm{\Delta }{x}_N}{dt}=-a_N-N_j\frac{\mathrm{\partial }{S}_d}{\mathrm{\partial }{x}_j}=-a_N^{eff}$(13) , a positive $a_N^{eff}$ acts to increase the distance between two adjacent $c$ iso-surfaces and hence reduces the value of $|\mathrm{\nabla }c|$ and indicates flame thickening. On the other hand, the negative values of $a_T^{eff}$ indicates a fractional reduction in flame area, as can be seen from eq. (14)(14) $\frac{1}{\delta A}\frac{d(\delta A)}{dt}=a_T+2S_d{\kappa }_m=a_T^{eff}$(14) . This reduction in flame area could be partially attributed to strong flame self-interactions owing to flame front curvatures. Figure 11 shows that changing the dilution level and mixture inhomogeneity within the parameter range considered in this study do not have significant impacts on $⟨a_N^{eff}{⟩}_c$ and $⟨a_T^{eff}{⟩}_c$.

Figure 11. Profiles of the normalised mean values of the normal and tangential effective strain rates (i.e. $a_N^{eff}$ and $a_T^{eff}$, respectively) conditioned upon $c$ for all cases considered here.

The balance of convection, diffusion, and reaction rates

Minamoto et al. (2014a) defined the variable $\mathcal{B}$ (eq. (15))(15) $\mathcal{B}\equiv |\mathcal{C}-\mathcal{D}|-|\mathcal{R}|$(15) to assess the balance of convection, diffusion, and chemical reaction rates in the reaction progress variable equation (eq. (5)(5) $\rho \frac{\mathrm{\partial }c}{\mathrm{\partial }t}+\rho u_j\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\rho D_c\frac{\mathrm{\partial }c}{\mathrm{\partial }{x}_j}\right)+{\dot{w}}_c+A_{\xi }$(5) ). Minamoto et al. (2014a) concluded that the regions associated with $\mathcal{B}<0$ are dominated by autoignition and flame interactions which can lead to small levels of convection ($\mathcal{C}$) and diffusion ($\mathcal{D}$) and hence negative values for $\mathcal{B}$. These regions (where $\mathcal{B}<0$) are termed as the reaction-dominated regions by Minamoto et al. (2014a). The regions with $\mathcal{B}>0$ are dominated by propagating flames, particularly, due to the large convection contribution. Thus, two distinct regions were identified based on the value of $\mathcal{B}$: 1) reaction dominated at $\mathcal{B}<0$; 2) propagating-flame region at $\mathcal{B}>0$. Since both regions co-exist in MILD combustion, it is worthwhile to investigate the contribution of each region to the mean SDF, $S_d$ and the strain rates due to flame propagation ($N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$ and $2S_d{\kappa }_m$). Here, only the additional strain rates due to flame propagation are considered since it has been shown in the discussion of Figure 11 that these strain rates are an order of magnitude larger than the strain rates due to the background fluid motion and thus $⟨2S_d{\kappa }_m{⟩}_c$ and $⟨N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j{⟩}_c$ are the sole drivers of the behaviors of $⟨a_{{ }_T}^{eff}{⟩}_c$ and $⟨a_{{ }_N}^{eff}{⟩}_c$, respectively. Moreover, the effect of mixture inhomogeneity was not explicitly taken into consideration when defining $\mathcal{B}$ since it was shown to be negligible in the previous sections (see the discussions of Figures 9b and 10b).

Figure 12 shows the PDF of ${\mathcal{B}}^{+}\equiv \mathcal{B}\times ({\mathrm{\ell }}_{{ }_0}/{\rho }_{{ }_0}{u}_{{ }_0})$ for all cases considered here. It can be seen from Figure 12 that the PDF of ${\mathcal{B}}^{+}$ is broad with the peak $P({\mathcal{B}}^{+})$ value close to ${\mathcal{B}}^{+}\approx 0$. Figure 12 also shows that ${\mathcal{B}}^{+}$ is predominately positive indicating the prevalence of the propagating-flame-dominated regions throughout the domain. However, a negative branch of ${\mathcal{B}}^{+}$ persist in all cases indicating the presence of reaction-dominated regions. Figure 12 shows that increasing the dilution levels leads to a reduced probability of negative ${\mathcal{B}}^{+}$ compared to the low dilution case (HM-A) due to the reduced reaction rate (see Figure 6c). On the other hand, the mixture inhomogeneity gives rise to slightly larger negative region of ${\mathcal{B}}^{+}$ compared to the homogeneous case (HM-A) which can also be attributed to slightly increased reaction rate, as shown in Figure 6b.

Figure 12. Pdfs of ${\mathcal{B}}^{+}\equiv \mathcal{B}\times ({\mathrm{\ell }}_{{ }_0}/{\rho }_{{ }_0}{u}_{{ }_0})$ within the region given by $0.05<c<0.95$ for all cases considered here.

Figure 13. Iso-surfaces indicating the regions of negative ${\mathcal{B}}^{+}$ within the range given by $0.05<c<0.95$ for all cases considered here.

The spatial distribution of the regions where are shown in Figure 13. Figure 13 emphasizes the findings from Figure 12, where increasing the dilution level leads to reduction in the volume of negative ${\mathcal{B}}^{+}$ regions, while mixture inhomogeneity gives rise to larger ${\mathcal{B}}^{+}<0$ regions. It is also worth noting that the length-scale separation in the ${\mathcal{B}}^{+}$ field does not change significantly for the cases considered here, with the integral length scale calculated from the ${\mathcal{B}}^{+}$ field having the values $5.06\times {10}^{-4}$ m, $5.23\times {10}^{-4}$ m, and $2.0\times {10}^{-4}$ m for the HM-A, HM-B, and IM-A cases, respectively. Figure 13 also shows that ${\mathcal{B}}^{+}<0$ islands are rare for the conditions analyzed in this study, which prompted the use of multiple snapshots to gather the conditional statistics.

Figure 14 shows the profiles of normalized mean values of $|\mathrm{\nabla }c|$ conditioned upon $c$ for both the reaction dominated ($\mathcal{B}<0$) and propagating-flame ($\mathcal{B}>0$) regions. It can be seen from Figure 14 that, in the homogeneous cases, the reaction-dominated regions yield a noticeably reduced mean value for the SDF compared to that resulting in the propagating-flame region. Thus, one can conclude that on average the flame thickening is slightly more pronounced in the reaction-dominated regions. The occurrences of self interactions of flame fronts lead to small values of $|\mathrm{\nabla }c|$ in these regions (Griffiths et al. 2015; Trivedi, Nivarti, and Cant 2019) and thus the mean values of $|\mathrm{\nabla }c|$ conditional upon $c$ for $\mathcal{B}<0$ are smaller than those for $\mathcal{B}>0$. In the inhomogeneous case, the above trend is observed up to $c\approx 0.5$. For $c>0.5$ values, both the reaction-dominated and propagating-flame regions give rise to similar mean values of the SDF.

Figure 14. Profiles of the mean values of $|\mathrm{\nabla }c|\times {\delta }_{\mathrm{\ell }}$ at both the reaction dominated (i.e. $\mathcal{B}<0$) and propagating-flame (i.e. $\mathcal{B}>0$) regions conditioned upon $c$ for all cases considered here. The ’LAM’ curve refers to the diluted laminar flame used to create the initial MILD field following the methodology by Minamoto et al. (2013).

The profiles of normalized mean displacement speed conditioned upon $c$ in both the reaction-dominated and propagating-flame regions are shown in Figure 15. It can be seen from Figure 15 that, in general, the mean displacement speed in the reaction-dominated regions remain larger than that in the propagating-flame regions. However, the differences between the mean values of displacement speed in these two regions remain small in the homogeneous cases, particularly, in the high dilution case (HM-B) where this difference in the mean values of $S_d$ becomes almost negligible for $c⪆0.6$. This is not the case though in the inhomogeneous case (IM-A) where the mean values of $S_d$ in the reaction-dominated regions can be significantly larger than that in the propagating-flame regions. There is a variation of equivalence ratio within the reaction zone in the inhomogeneous mixture, and the presence of more reactive mixtures than the mixture with global mean equivalence ratio could potentially give rise to larger mean values of $S_d$ in these regions than in the propagating-flame regions where the greater reaction rate magnitudes for more reactive mixtures are partially countered by larger magnitudes of molecular diffusion rate effects. However, Figure 16 shows that the mean equivalence ratio $\varphi$ conditioned upon $c$ remains mostly below the stoichiometric level indicating a lean mixture over the most of $c$ range with the exception of small $c$ values ($c⪅0.2$). This remains true for both the propagating-flames regions ($\mathcal{B}>0$) and reaction-dominated regions $\mathcal{B}<0$ with only a slight increase in $⟨\varphi {⟩}_c$ observed for $c⪆0.4$ in the regions where $\mathcal{B}<0$.

Figure 15. Profiles of the mean values of $S_d$ at both the reaction-dominated (i.e. $\mathcal{B}<0$) and propagating-flame (i.e. $\mathcal{B}>0$) regions conditioned upon $c$ for all cases considered here.

Figure 16. Profiles of the mean equivalence ratio ϕ at both the reaction-dominated (i.e. $\mathcal{B}<0$) and propagating-flame (i.e. $\mathcal{B}>0$) regions conditioned upon $c$ for the inhomogeneous case (IM-A).

Figure 17 shows the normalized mean values of the additional flame normal strain rate due to flame propagation $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$ conditioned upon $c$ at both the reaction-dominated and propagating-flame regions. It can be seen from Figure 17 that the mean value of $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$ in the reaction-dominated regions (i.e. $\mathcal{B}<0$) remains greater than that in the propagating-flame regions (i.e. $\mathcal{B}>0$) in all cases for $c⪅0.6$. This trend continues in the high dilution case (HM-B) throughout the flame front. However, in the low dilution cases, the value of $⟨N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j{⟩}_c$ in the reaction-dominated region decreases at $c⪆0.6$ to become smaller than that in the propagating flame region. The increased levels in $⟨N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j{⟩}_c$ in the reaction-dominated regions across most of the flame front ($c⪅0.6$) can explain the corresponding reduced levels of the SDF at the same $c$ range. In the low dilution cases, the reduction in $⟨N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j{⟩}_c$ in the reaction-dominated regions at $c⪆0.6$ is also reflected with an increase in the SDF value in the corresponding $c$ range, as can be clearly seen in Figure 14b.

Figure 17. Profiles of the normalised mean values of normal strain rate due to flame propagation $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$ at both the reaction-dominated (i.e. $\mathcal{B}<0$) and propagating-flame (i.e. $\mathcal{B}>0$) regions conditioned upon $c$ for all cases considered here.

The profiles of the normalized mean values of additional tangential strain rate due to flame propagation $2S_d{\kappa }_m$ conditioned upon $c$ at both the reaction-dominated and propagating-flame regions are shown in Figure 18. In all cases, the mean value of $a_T^{eff}$ in the reaction-dominated region has a larger negative value than that in the propagating flame region across most of the flame front. However, in the low dilution cases (HM-A and IM-A), the mean value of $2S_d{\kappa }_m$ increases with the progress variable from $c\approx 0.75$ and even becomes positive at high $c$ values. This could indicate slightly increased levels of flame area destruction in the reaction-dominated region (compared to that in the propagating-flame region) across most of the flame front, but local flame area generation is possible at high $c$ values in the inhomogeneous case.

Figure 18. Profiles of the mean values of tangential strain rate due to flame propagation $2S_d{\kappa }_m$ at both the reaction-dominated (i.e. $\mathcal{B}<0$) and propagating-flame (i.e. $\mathcal{B}>0$) regions conditioned upon $c$ for all cases considered here.

Implications on modelling

As can be seen from eqs. (12)(12) $\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }t}+\frac{\mathrm{\partial }(u_j|\mathrm{\nabla }c|)}{\mathrm{\partial }{x}_j}=a_T|\mathrm{\nabla }c|-\frac{\mathrm{\partial }\left(S_d{N}_j|\mathrm{\nabla }c|\right)}{\mathrm{\partial }{x}_j}+2S_d{\kappa }_m|\mathrm{\nabla }c|$(12) to (14(14) $\frac{1}{\delta A}\frac{d(\delta A)}{dt}=a_T+2S_d{\kappa }_m=a_T^{eff}$(14) ), the statistical behaviors of the various strain rates ($a_N$, $a_T$, $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$, $2S_d{\kappa }_m$, $a_N^{eff}$, and $a_T^{eff}$) affect the evolution of SDF, and thus can have a significant influence on flamelet-type modeling in Large Eddy Simulation (LES) and Reynolds Averaged Navier–Stokes (RANS) simulations. This is particularly the case when using the generalized flame surface density (FSD) (i.e. ${\mathrm{\Sigma }}_{gen}=\overline{|\mathrm{\nabla }c|}$) and the Scalar Dissipation Rate (SDR) (i.e. $N_c=D_c|\mathrm{\nabla }c{|}^2$) based modeling approaches. To illustrate the impact of the various strain rates on the generalized FSD, it is helpful to rewrite eq. (13)(13) $\frac{1}{|\mathrm{\nabla }c|}\frac{d|\mathrm{\nabla }c|}{dt}=-\frac{1}{\mathrm{\Delta }{x}_N}\frac{d\mathrm{\Delta }{x}_N}{dt}=-a_N-N_j\frac{\mathrm{\partial }{S}_d}{\mathrm{\partial }{x}_j}=-a_N^{eff}$(13) as follows (Klein, Alwazzan, and Chakraborty 2018):

(18) $\begin{array}{cc}\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }t}+u_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}& =-a_N|\mathrm{\nabla }c|-N_j\frac{\mathrm{\partial }{S}_d}{\mathrm{\partial }{x}_j}|\mathrm{\nabla }c|-S_d{N}_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}\\ \hspace{0.17em}& =-a_N^{eff}|\mathrm{\nabla }c|-S_d{N}_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}\end{array}$(18)

The effect of $a_T$, $2S_d{\kappa }_m$ and $a_T^{eff}$ on $|\mathrm{\nabla }c|$ can also be shown by writing eq. (12)(12) $\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }t}+\frac{\mathrm{\partial }(u_j|\mathrm{\nabla }c|)}{\mathrm{\partial }{x}_j}=a_T|\mathrm{\nabla }c|-\frac{\mathrm{\partial }\left(S_d{N}_j|\mathrm{\nabla }c|\right)}{\mathrm{\partial }{x}_j}+2S_d{\kappa }_m|\mathrm{\nabla }c|$(12) as follows (Awad et al. 2021):

(19) $\begin{array}{cc}\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }t}+\frac{\mathrm{\partial }(u_j|\mathrm{\nabla }c|)}{\mathrm{\partial }{x}_j}& =a_T|\mathrm{\nabla }c|+2S_d{\kappa }_m|\mathrm{\nabla }c|-\frac{\mathrm{\partial }\left(S_d{N}_j|\mathrm{\nabla }c|\right)}{\mathrm{\partial }{x}_j}\\ \hspace{0.17em}& =a_T^{eff}|\mathrm{\nabla }c|-\frac{\mathrm{\partial }\left(S_d{N}_j|\mathrm{\nabla }c|\right)}{\mathrm{\partial }{x}_j}\end{array}$(19)

A transport equation for ${\mathrm{\Sigma }}_{gen}$ can then be obtained by Reynolds averaging/LES filtering eq. (18)(18) $\begin{array}{cc}\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }t}+u_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}& =-a_N|\mathrm{\nabla }c|-N_j\frac{\mathrm{\partial }{S}_d}{\mathrm{\partial }{x}_j}|\mathrm{\nabla }c|-S_d{N}_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}\\ \hspace{0.17em}& =-a_N^{eff}|\mathrm{\nabla }c|-S_d{N}_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}\end{array}$(18) or eq. (19)(19) $\begin{array}{cc}\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }t}+\frac{\mathrm{\partial }(u_j|\mathrm{\nabla }c|)}{\mathrm{\partial }{x}_j}& =a_T|\mathrm{\nabla }c|+2S_d{\kappa }_m|\mathrm{\nabla }c|-\frac{\mathrm{\partial }\left(S_d{N}_j|\mathrm{\nabla }c|\right)}{\mathrm{\partial }{x}_j}\\ \hspace{0.17em}& =a_T^{eff}|\mathrm{\nabla }c|-\frac{\mathrm{\partial }\left(S_d{N}_j|\mathrm{\nabla }c|\right)}{\mathrm{\partial }{x}_j}\end{array}$(19) , and the resulting FSD transport equation will include the various strain rates affecting the behavior of SDF. Furthermore, a transport equation for the SDR can be obtained by multiplying eq. (18)(18) $\begin{array}{cc}\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }t}+u_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}& =-a_N|\mathrm{\nabla }c|-N_j\frac{\mathrm{\partial }{S}_d}{\mathrm{\partial }{x}_j}|\mathrm{\nabla }c|-S_d{N}_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}\\ \hspace{0.17em}& =-a_N^{eff}|\mathrm{\nabla }c|-S_d{N}_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}\end{array}$(18) or eq. (19)(19) $\begin{array}{cc}\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }t}+\frac{\mathrm{\partial }(u_j|\mathrm{\nabla }c|)}{\mathrm{\partial }{x}_j}& =a_T|\mathrm{\nabla }c|+2S_d{\kappa }_m|\mathrm{\nabla }c|-\frac{\mathrm{\partial }\left(S_d{N}_j|\mathrm{\nabla }c|\right)}{\mathrm{\partial }{x}_j}\\ \hspace{0.17em}& =a_T^{eff}|\mathrm{\nabla }c|-\frac{\mathrm{\partial }\left(S_d{N}_j|\mathrm{\nabla }c|\right)}{\mathrm{\partial }{x}_j}\end{array}$(19) by $2|\mathrm{\nabla }c|$ and some algebraic manipulations. The resulting equation can be written as (Ahmed et al. 2020; Klein, Alwazzan, and Chakraborty 2018):

(20) $\begin{array}{cc}\rho \frac{\mathrm{\partial }{N}_c}{\mathrm{\partial }t}+\rho u_j\frac{\mathrm{\partial }{N}_c}{\mathrm{\partial }{x}_j}& =\stackrel{-2\rho a_N^{eff}{N}_c}{\overbrace{-2\rho a_N{N}_c-2\rho N_j\frac{\mathrm{\partial }{S}_d}{\mathrm{\partial }{x}_j}{N}_c}}-\rho S_d{N}_j\frac{\mathrm{\partial }{N}_c}{\mathrm{\partial }{x}_j}\\ \hspace{0.17em}& +\rho S_d{N}_j{N}_c\frac{1}{D_c}\frac{\mathrm{\partial }{D}_c}{\mathrm{\partial }{x}_j}+\frac{\rho N_c}{D_c}\left(\frac{\mathrm{\partial }{D}_c}{\mathrm{\partial }t}+u_j\frac{\mathrm{\partial }{D}_c}{\mathrm{\partial }{x}_j}\right)\end{array}$(20)

It is evident from eqs. (18)(18) $\begin{array}{cc}\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }t}+u_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}& =-a_N|\mathrm{\nabla }c|-N_j\frac{\mathrm{\partial }{S}_d}{\mathrm{\partial }{x}_j}|\mathrm{\nabla }c|-S_d{N}_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}\\ \hspace{0.17em}& =-a_N^{eff}|\mathrm{\nabla }c|-S_d{N}_j\frac{\mathrm{\partial }|\mathrm{\nabla }c|}{\mathrm{\partial }{x}_j}\end{array}$(18) to (20(20) $\begin{array}{cc}\rho \frac{\mathrm{\partial }{N}_c}{\mathrm{\partial }t}+\rho u_j\frac{\mathrm{\partial }{N}_c}{\mathrm{\partial }{x}_j}& =\stackrel{-2\rho a_N^{eff}{N}_c}{\overbrace{-2\rho a_N{N}_c-2\rho N_j\frac{\mathrm{\partial }{S}_d}{\mathrm{\partial }{x}_j}{N}_c}}-\rho S_d{N}_j\frac{\mathrm{\partial }{N}_c}{\mathrm{\partial }{x}_j}\\ \hspace{0.17em}& +\rho S_d{N}_j{N}_c\frac{1}{D_c}\frac{\mathrm{\partial }{D}_c}{\mathrm{\partial }{x}_j}+\frac{\rho N_c}{D_c}\left(\frac{\mathrm{\partial }{D}_c}{\mathrm{\partial }t}+u_j\frac{\mathrm{\partial }{D}_c}{\mathrm{\partial }{x}_j}\right)\end{array}$(20) ) that the strain rates $a_N$, $a_T$, $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$, $2S_d{\kappa }_m$, $a_N^{eff}$ and $a_T^{eff}$ have a significant influence on the evolutions of FSD and SDR. The present analysis showed that the effect of mixture inhomogeneity and dilution level on the normalized mean values of $a_N$, $a_T$, $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$, $2S_d{\kappa }_m$, $a_N^{eff}$ and $a_T^{eff}$ were small in MILD combustion within the ranges considered in this study. Thus, when extending the existing FSD and SDR closures for the modeling of the EGR-type MILD combustion, the effect of mixture inhomogeneity and dilution levels may not need to be explicitly accounted for within the parameter ranges considered in this study.

Conclusions

A comparison of the mean behaviors of the magnitude of the reaction progress variable gradient (surface density function, SDF), displacement speed, and flame normal and tangential strain rates in EGR-type homogeneous and inhomogeneous mixture MILD combustion has been conducted using three-dimensional DNS with reduced chemical mechanism. The effects of dilution level and mixture inhomogeneity have been considered and the statistical behaviors of the SDF and $S_d$ were also characterized in both the reaction-dominated and propagating-flame regions of the domain.

It is found that the mean values of SDF decrease significantly under turbulent conditions indicating flame thickening compared to the unstretched laminar flame with the degree of the SDF reduction being less pronounced in the high dilution case. It is also found that the effect of mixture inhomogeneity was small on the SDF evolution for the parameters considered here.

The mean displacement speed decreases with $c$ across most of the flame front, but then starts to increase at $c\approx 0.6$ in the low dilution cases. The increase of the mean value of $S_d$ toward the burned gas side is found to be driven by an increase in the reaction rate of the progress variable. In the high dilution case, the mean value of $S_d$ continues to decrease almost monotonically with $c$ since increasing the dilution level gives rise to reduced levels of the reaction rate of the progress variable and thereby offers smaller reaction rate contribution to the mean values of displacement speed $S_d$.

The mean dilatation rate is found to be negligible due to the small heat release associated with MILD combustion. The mean values of the effective flame normal and tangential strain rates were driven by the mean values of the additional flame normal and tangential strain rates due to flame propagation, respectively. A generally positive mean effective normal strain rate explains the decrease in the SDF, while the negative effective tangential strain rate indicates a reduction in flame area that could be attributed to flame self-interaction.

In the homogeneous cases, the reaction-dominated regions yield a noticeably reduced mean value of the SDF compared to that resulting in the propagating-flame region, while this is only true toward the unburned gas side in the inhomogeneous case. The reaction-dominated regions of the domain also yield higher mean values of $S_d$ in all cases and generally enhanced negative mean values for $2S_d{\kappa }_m$, while little difference is observed between the reaction-dominated and propagating-flame regions in terms of mean values of $N_j\mathrm{\partial }{S}_d/\mathrm{\partial }{x}_j$.

In summary, the present study shows that the statistical behavior of the SDF is mostly driven by the additional strain rates due to flame propagation, and that mixture inhomogeneity and dilution level have limited effects on these additional strain rates within the parameter range considered here. This indicates that the molecular diffusion contributions of displacement speed, especially the curvature contributions, have to be accurately accounted for when extending the FSD/SDR modeling to EGR-type MILD combustion. This also suggests that the modeling approaches applicable to homogeneous MILD combustion can potentially be extended for the inhomogeneous cases. These findings are consistent with the recent observations by Swaminathan (2019), which indicate that MILD combustion of homogeneous and inhomogeneous mixtures behave similarly for comparable values of thermochemical and turbulence parameters. However, the present analysis is carried out for a simple canonical configuration without any mean shear effects and moderate values of turbulent Reynolds number. Thus, further analyses at higher values of turbulent Reynolds number and more complex configurations will be needed to confirm whether the conclusions drawn in this analysis for the reactive scalar gradient statistics for MILD combustion hold in other configurations and flow conditions. These aspects will form the foundations of further investigations.

Acknowledgements

The authors are grateful to EPSRC (EP/S025154/1) for the financial support. The computational support was provided by ARCHER2 (EP/R029369/1) and the HPC facility at Newcastle University (Rocket).

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The work was supported by the Engineering and Physical Sciences Research Council [EP/R029369/1, EP/S025154/1]