Comparison of the Reactive Scalar Gradient Evolution between Homogeneous MILD Combustion and Premixed Turbulent Flames

: Moderate or intense low-oxygen dilution (MILD) combustion is a novel combustion technique that can simultaneously improve thermal efﬁciency and reduce emissions. This paper strain rate assumed positive mean values, but it was overcome by negative mean values of curvature stretch rate to yield negative mean values of stretch rate for both the premixed ﬂames and MILD combustion cases. This behaviour is explained in terms of the curvature dependence of displacement speed. These ﬁndings suggest that the curvature dependence of displacement speed and the scalar gradient alignment with local principal strain rate eigendirections need to be addressed for modelling EGR-type homogeneous-mixture MILD combustion.


Introduction
Fossil fuels are expected to supply 85% of the primary energy required for transportation until 2040 [1]. Thus, these fuels will continue to be a major energy source for applications involving high energy density, such as in the transport sector [1], at least for a foreseeable future. Fossil fuel combustion is a major contributor to pollutants and greenhouse gas emissions. Thus, it is imperative for the new generation of combustion devices to be both energy-efficient and environmentally friendly. Moderate or intense low-oxygen dilution (MILD) combustion is one of the modern concepts with the potential to deliver both high-energy efficiency and ultra-low emission [2][3][4]. This concept involves a fuel-air mixture preheated above its autoignition temperature [4], which can be achieved by recirculating the hot flue gas of furnaces and boilers [5]. Such gas recirculation inherently increases the thermal efficiency of furnaces and boilers, as the maximum temperature rise remains smaller than the autoignition temperature [4,5]. MILD combustion provides uniform temperature distribution, extremely low levels of NOx, CO and greenhouse gases, high thermal efficiency, and excellent combustion stability. All of these make this technique attractive for power generation and process applications. MILD combustion is characterised by (a) high level of dilution by burnt gases such that the temperature rise ∆T remains smaller than the autoignition temperature T ign (i.e., ∆T < T ign ), and (b) the reactants are significantly preheated such that the reactant temperature T 0 remains greater than T ign (i.e., T 0 > T ign ) [4]. This definition of MILD combustion (alternatively known as flameless combustion) follows from the well-stirred reactor (WSR) theory, which can be readily applied to turbulent premixed combustion. MILD combustion has been analysed extensively [4][5][6][7][8][9][10][11][12][13][14][15] by experimental means to understand the underlying combustion process. Plessing et al. [6] compared MILD combustion with conventional combustion by investigating both planar laser-induced fluorescence images of OH radicals (OH-PLIF) and Rayleigh thermometry. The presence of flame fronts in MILD combustion was observed using OH-PLIF but the temperature measurements revealed indications of distributed burning. A similar observation was also reported by Ozdemir and Peters [10] based on their experimental investigations. Several Reynolds averaged Navier-Stokes (RANS) and large eddy simulations (LES) investigations of MILD combustion using conventional flamelet and eddy dissipation approaches [8,9,13,16,17] have also been conducted which reported 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 (e.g., CO and OH) concentrations [8,12]. The advancements in high-performance computing have enabled direct numerical simulations (DNS) of MILD combustion [18][19][20][21][22][23][24][25][26][27]. Van Oijen [18] carried out DNS of autoigniting mixing layers of CH 4 /H 2 fuel-blends and provided important insights to the roles of non-unity Lewis number and O 2 dilution on autoignition delay under MILD combustion conditions. Minamoto, Swaminathan and co-workers [18][19][20][21][22] analysed the distribution of species, temperature, and reaction-rate fields [18,19], morphological and topological structures [20] and scalar gradient statistics in terms of temperature and equivalent OH-PLIF signal [21] constructed from DNS data based on three-dimensional simulations of the combustion an exhaust gas recirculation (EGR)-type homogeneous-mixture under MILD conditions. They concluded that the appearance of a distributed combustion is due to the interaction of thin reaction zones [19][20][21][22][23], which were also reported based on OH-PLIF visualisations [4,10,14,15]. These analyses have been further extended by Doan, Swaminathan and co-workers for stratified-mixture MILD combustion [24][25][26] and provided important insights into the flame structure [24,25], reaction-zone topology [24,25] and markers of combustion mode [26]. A review of the physical insights gained from DNS can be found in [27]. These DNS analyses [19][20][21][22][23][24][25][26] revealed that the convolution of the reaction zones in MILD combustion is not only affected by turbulence-chemistry interaction but the dilution factor also plays a significant role. Furthermore, recent studies [20,23] using DNS data have indicated that flamelet-based models can be applied to gaseous-fuel MILD combustion with some modifications, which is contrary to the conventional expectation that flamelet closures are unlikely to be valid for Energies 2021, 14, 7677 3 of 28 distributed reactions expected in MILD conditions. As flamelet models need the statistics of the reactive scalar gradient, it is important to investigate the statistical behaviour of the gradient of the reaction progress variable, c (surface density function (SDF) =|∇c| [28]), since it is closely related to the scalar dissipation rate (SDR) N c = D c | ∇c| 2 (with D c being the reaction progress variable diffusivity) [29,30] and generalised flame surface density (FSD) ( Σ gen = ∇c with the overbar, suggesting a Reynolds averaging/LES filtering operation, as appropriate) [31], which are often used for the modelling of the combustion of homogeneous-mixtures. In this respect, it is important to understand the influence of the strain rates induced by fluid motion and flame propagation on SDF evolution in homogeneous-mixture MILD combustion, and the differences in these mechanisms in comparison with conventional premixed combustion. Although the SDF evolution within the flame front and the strain rates affecting this process have been extensively analysed, in the case of conventional premixed flames [28,[32][33][34][35][36][37][38][39][40][41][42][43][44], this aspect has rarely been addressed in the case of EGR-type homogeneous-mixture MILD combustion [23]. The present analysis aims to address this gap in the existing literature by conducting three-dimensional DNS with a skeletal methane-air chemical mechanism [45] for EGR-type homogeneousmixture combustion under MILD conditions, in the case of an O 2 dilution percentage of 4.8% by volume, for an equivalence ratio of φ = 0.8 under different root-mean-square turbulent-velocity fluctuations u . The SDF statistics of these MILD combustion cases have been compared to the corresponding statistics obtained from DNS of undiluted CH4-air premixed flames of φ = 0.8 with an unburned gas temperature of T 0 = 300 K following the methodology adopted previously by Minamoto, Swaminathan and co-workers [19][20][21][22][23] to compare the flame statistics between EGR-type homogeneous-mixture MILD combustion and conventional premixed combustion. In this respect, the main objectives of the current analysis are: 1 A comparison between the SDF statistics in EGR-type homogeneous-mixture MILD combustion and conventional premixed combustion to understand how different strain rates influence the SDF evolution in MILD combustion under different turbulence intensities. 2 Modelling implication of the differences in the SDF statistics between MILD and conventional premixed turbulent combustion.
The rest of the paper will be organised as follows. The mathematical background pertaining to this analysis will be presented in the next section. This will be followed by a brief discussion of the numerical implementation. Following that, the results will be presented and subsequently discussed. Finally, key findings will be summarised, and conclusions will be drawn.

Mathematical Background
The extent of the completion of the chemical reaction can be quantified in terms of a reaction progress variable c, which increases monotonically from 0.0, in the unburned gas, to 1.0 in the fully burned products. For the present analysis, the reaction progress variable c is defined based on the fuel mass fraction Y F (i.e., CH 4 mass fraction) in the following manner: where subscripts R and P refer to values in unburned reactants and fully burned products. The reaction progress variable transport is governed by [46]: where, ρ is the gas density, D c is the mass diffusivity of the species based on which the reaction progress variable is defined, u j is the j th velocity component, and .
ω c is the reaction ω c is expressed as follows [46]: where . ω F is the fuel reaction rate. Equation (2) can be rewritten for a given c iso-surface in the following manner [46]: where S d is the displacement speed at a given c iso-surface and it is given by: The displacement speed arises because of the relative balance between its components, and can be split in the following manner [46][47][48]: Here, S r , S n and S t are the reaction, normal diffusion and tangential diffusion components of displacement speed, respectively, which are defined as [48][49][50]: where, → N = −∇c/|∇c| is the flame-normal vector and κ m = 0.5∇· → N is the mean flame curvature. According to these definitions, the flame-normal points towards the unburned gas and a reaction progress variable iso-surface has a positive (negative) curvature when it is convexly (concavely) curved towards the reactants.
The transport equation for the magnitude of the reaction progress variable gradient (i.e., surface density function (SDF) [28]) for a given c-isosurface takes the following form [31][32][33][34][35][36][37][38][39][40][41][42][43][44]: where, a T |∇c| is the SDF strain rate term with a T = δ ij − N i N j ∂u i /∂x j being the tangential strain rate, −∂ S d N j |∇c| /∂x j arises due to flame-normal propagation and thus referred to as the SDF propagation term, 2S d κ m |∇c| originates due to curvature and therefore is referred to as the SDF curvature term. Equation (8) can be rewritten as [37][38][39]: where, d(. . .)/dt = ∂(. . .)/∂t + V c j ∂(. . .)/∂x j is the total derivative in the reference frame attached to the flame, V c j = u j + S d N j is the jth component of flame propagation velocity, ∆x N is the distance between two adjacent c-isosurfaces, a N = N i N j ∂u i /∂x j is the fluid-dynamic normal strain rate, N j ∂S d /∂x j is the normal strain rate induced by flame propagation and a e f f N = a N + N j ∂S d /∂x j is the effective normal strain rate [37][38][39]. It is evident from Equation (9) that a positive value of a e f f N acts to increase the distance between two adjacent c-isosurfaces, which in turn decreases |∇c|. By contrast, a negative value of a e f f N decreases the distance between two adjacent c-isosurfaces, and therefore increases |∇c|.
As |∇c| is closely related to the surface area of the reaction zone, it is worthwhile to consider the evolution of the elemental reaction zone surface area δA [37][38][39]49,50]: Here, a e f f T is commonly referred to as the flame stretch rate or the effective tangential strain rate, where a T is the flow component of a e f f T and the component 2S d κ m arises due to flame propagation and therefore is commonly termed as the curvature stretch rate.
It is evident from the above discussion that a N , a T , and a e f f T determine the evolutions of the SDF and the area of the reaction zone and thus the statistical behaviours of these strain rates will be discussed in detail in Section 4 of this paper.

Numerical Implementation
A compressible DNS code SENGA2 [51] is used for the simulations conducted in the current study. In SENGA2, the standard conservation equations of mass, momentum, energy and species mass fractions for compressible reacting flows are solved in dimensional form. The thermophysical properties, such as viscosity, thermal conductivity, specific heat capacities and mass diffusivities, are taken to be temperature-dependent and are expressed in terms of polynomials similar in form to CHEMKIN polynomials [51]. A 10th-order central difference scheme is used 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 boundaries. A fourth-order explicit low-storage Runge-Kutta scheme is used for time advancement. The boundary conditions are specified according to the Naiver-Stokes characteristic boundary condition (NSCBC) methodology [52]. A skeletal chemical mechanism comprising 16 species and 25 reactions (including 10 reversible reactions) has been taken to represent the chemical kinetics of CH 4 / Air combustion [45]. A domain, of size 10 mm × 10 mm × 10 mm and discretized by a Cartesian grid of 252 × 252 × 252 with uniform grid spacing, has been used for the simulations of MILD CH 4 -air mixtures with an O 2 dilution percentage of 4.8% by volume for an equivalence ratio of φ = 0.8 and an initial temperature of about 1500 K (i.e., T 0 ≈ 1500 K). The unburned gas temperature of 1500 K is higher than the autoignition temperature (i.e., T ign = 1100 K based on wellstirred reactor calculations [19][20][21][22][23]) for CH 4 -air mixtures with an O 2 dilution percentage of 4.8% by volume for an equivalence ratio of φ = 0.8. The simulation domain consists of a turbulent inflow with specified density, velocity and species at the left-hand x 1 −boundary and a partially non-reflective outflow boundary condition is specified at the right-hand x 1 − boundary. All other boundaries are taken to be periodic.
For the purpose of comparison, DNS calculations of statistically planar conventional turbulent premixed flames of φ = 0.8 with an unburned gas temperature of 300 K have been conducted for a domain size of 20 mm × 10 mm × 10 mm, which was discretized by a uniform Cartesian grid of 504 × 252 × 252. A longer domain in the streamwise flow direction was taken for the turbulent premixed flame simulations to allow for flame propagation in the mean direction of flame propagation (which is taken to align with the x 1 − direction) without the influence of computational boundaries. In the premixed flame simulations, the domain boundaries in the direction of the mean flame propagation were taken to be partially non-reflecting, and transverse boundaries were taken to be periodic.
The grid spacing for both the MILD and premixed combustion simulations ensures that the thermal flame thickness δ th = (T P − T 0 )/max|∇T| L (with T, T p and T 0 being the instantaneous, product and reactant temperatures, respectively and subscript L refers to the values in the corresponding 1D unstretched premixed flame) is resolved using at least 12 grid points. The grid spacing also ensured the resolution of the Kolmogorov length scale, η, by at least 1.5 grid points for all turbulent cases.
The thermochemical conditions in the unburned gas, including the mole fractions of the reactants, and unburned gas temperature for the MILD and the lean premixed cases are shown in Table 1. Under the conditions summarised in Table 1 Karlovitz number Ka = (u /S L ) 3/2 (l/δ th ) −1/2 are summarised in Table 2. The turbulence levels reported in Table 2 matches and expands upon those used by Minamoto et al. [19][20][21][22] and represent a jet in hot co-flow (JHC)-type flame configuration. This provides a base line to validate the current DNS dataset against, while also expanding the investigation to include the effect of turbulence intensities on homogeneous MILD combustion. It is worth noting the identical values of u /S L and l/δ th between the MILD and premixed cases ensure that they are nominally in the same location on the Borghi-Peters diagram [53] and all the cases considered here nominally represent combustion within the thin reaction zones regime [53]. It is important to note that the turbulent Reynolds number Re t = u l/ν with ν being the kinematic viscosity scales as: Re t ∼ (u /S L )(l/δ th ). Thus, the identical values of u /S L and l/δ th between the MILD and premixed cases also ensure that they have comparable values of Re t . The MILD combustion simulations are split into two stages, which are summarised in Figure 1. The first stage corresponds to the mixing of the products of combustion with the reactant mixture. The second stage simulates the combustion process [19][20][21][22][23]. The species field obtained from the first stage acts as the initial, as well as the inflow fields for the second stage. The first stage mimics what happens inside the recirculation zone of a real MILD combustor. The mixture at stage 1 is prepared according to the following steps, as reported by Minamoto, Swaminathan and co-workers [19][20][21][22][23]: • A homogenous isotropic field, generated using a well-known pseudo-spectral method [54], is used to initialize turbulent velocity fluctuations following the Batchelor-Townsend spectrum [55].

•
The thermochemical conditions shown in Table 1 are used to simulate a 1D laminar premixed flame for φ = 0.8. An initial bimodal distribution of c Y = 1 − Y F /Y FR with prescribed integral length scale l c = 1.50 mm is then developed using the methodology by Eswaran and Pope [56]. The 1D laminar profiles of a species mass fraction for φ = 0.8 are specified as a function of the progress variable c Y . These functions, alongside the bimodal distribution of c Y , are used to initialise the density and species mass fractions corresponding to an atmospheric pressure and an unburned gas temperature of 1500 K (i.e., T 0 = 1500 K).

•
The generated bimodal scalar fields in the previous step are then allowed to evolve with turbulence for about 1 turnover time (l/u ) in a periodic domain mimicking the EGR in MILD combustion. At the end of this step the mean and variance of the preprocessed c field are c ≈ 0.50 and c 2 ≈ 0.09 respectively where Q is the The prepared mixture at the end of previous stage will then be fed into the domain through the inlet with a mean velocity U in = 20 m/ sec by scanning a plane through it, as used previously by Minamoto, Swaminathan and co-workers [19][20][21][22][23].
the EGR in MILD combustion. At the end of this step the mean and variance of the preprocessed field are 〈 〉 ≈ 0.50 and 〈 ′ 〉 ≈ 0.09 respectively where 〈 〉 is the mean value evaluated over the whole domain. The temperature in the preprocessed mixture has a variation of about ±0.3% of its mean value.  The prepared mixture at the end of previous stage will then be fed into the domain through the inlet with a mean velocity = 20 m/sec by scanning a plane through it, as used previously by Minamoto, Swaminathan and co-workers [19][20][21][22][23].
The species mass fractions, temperature and density for turbulent premixed flame simulations are initialised by the steady unstretched 1 D laminar premixed flame with = 0.8. A pseudo-spectral method [54] is used to initialise turbulent velocity fluctuations by a homogeneous isotropic field with prescribed values of ′ and , according to Batchelor-Townsend spectrum [55].  The species mass fractions, temperature and density for turbulent premixed flame simulations are initialised by the steady unstretched 1 D laminar premixed flame with φ = 0.8. A pseudo-spectral method [54] is used to initialise turbulent velocity fluctuations by a homogeneous isotropic field with prescribed values of u and l, according to Batchelor-Townsend spectrum [55].

Flame Morphology
The instantaneous views of the c = 0.8 isosurface, for both statistically planar turbulent premixed flames and homogeneous MILD combustion cases and for different turbulence intensities are shown in Figure 2. For methane-air combustion, the maximum heat release takes place close to c = 0.8 for the unstretched laminar premixed flame, and thus the c = 0.8 isosurface can be taken as a marker of the flame surface. It can be seen, from Figure 2, that the flame wrinkling increases with increasing u /S L for statistically planar turbulent premixed flames and, in these cases, the flame surface cleanly separates the reactants from the products, but the increased flame wrinkling, with an increase in u /S L , may give rise to self-interactions between flame surface wrinkles that can lead to localised pockets. Figure 2 reveals that there are significant differences in flame morphologies between the homogeneous MILD combustion and premixed flame cases. There is a significant amount of self-interaction of the c = 0.8 isosurface for all values of turbulence intensity in the MILD combustion cases and the extent of self-interaction increases with increasing turbulence intensity. The contours of temperature, T, in the central midplane for both the MILD and premixed cases considered here are shown in Figure 3, which reveals that the conventional premixed flame cases show a steep temperature gradient (from a temperature of 300 K in the reactants to 1960 K in the products) within the flame and the temperature contours become increasingly wrinkled with an increase in u /S L . By contrast, the MILD combustion cases show a small change in temperature (∆T ≤ 150 K) and thus the temperature rise is much more gradual in these cases than in the conventional premixed flames. The equivalent OH-PLIF signal constructed from DNS data can be evaluated in the following form for 1000 K ≤ T ≤ 1800 K [57][58][59][60][61]: Here the value of β is taken to be 0.0 (i.e., β = 0.0) following previous studies [22]. The contours of S OH /[S OH ] max, L at the central midplane for the cases considered here are shown in Figure 4, where [S OH ] max, L is the maximum value of the OH-PLIF signal in the corresponding unstretched laminar premixed flame. It can be seen from Figure 4 that the OH-PLIF signal reveals a much thicker reaction zone in MILD combustion cases than in the corresponding conventional premixed flames, but there is still a clearly defined interface separating reactants and products for turbulent MILD combustion. The observations made from temperature and the equivalent OH-PLIF signal from the DNS data are consistent with previous experimental [6][7][8][9][10][11][12][13][14][15] and numerical [19][20][21][22][23] findings.

Mean Behaviour of the SDF
The mean values of the normalised SDF (|∇c| × δ th ) conditioned upon c for both MILD and conventional premixed combustion cases, are shown in Figure 5, where the standard deviations of |∇c| × δ th conditioned upon c are also shown. In contrast to premixed combustion, which shows that the peak mean value of the SDF is skewed towards the burned gas side, the MILD combustion cases show a |∇c| × δ th profile skewed towards the unburned gas side of the flame front. This is a consequence of the differences in reaction and diffusion contributions to the SDF evolution between MILD combustion cases and conventional premixed flames considered here. This is contingent on the maximum heatrelease zone within the flame and highlights the difference between conventional premixed cases, where the maximum heat release occurs towards the burned gas side, in the range c = 0.7-0.8, while in MILD combustion, the reaction contribution remains active for a broad range of reaction progress, variable across the flame front (see later in Figure 9). Moreover, the differences in molecular diffusion contribution (especially the flame tangential diffusion rate) between the premixed flames and MILD combustion cases also contribute to this (see Figure 9 later in this paper).

Mean Behaviour of Fluid-Dynamic Strain Rates
The mean values of the normalised SDF (|∇ | × ) conditioned upon for bot MILD and conventional premixed combustion cases, are shown in Figure 5, where th standard deviations of |∇ | × conditioned upon are also shown. In contrast to pre mixed combustion, which shows that the peak mean value of the SDF is skewed toward the burned gas side, the MILD combustion cases show a |∇ | × profile skewed to wards the unburned gas side of the flame front. This is a consequence of the difference in reaction and diffusion contributions to the SDF evolution between MILD combustio cases and conventional premixed flames considered here. This is contingent on the max mum heat-release zone within the flame and highlights the difference between conven tional premixed cases, where the maximum heat release occurs towards the burned ga side, in the range = 0.7 − 0.8, while in MILD combustion, the reaction contribution re mains active for a broad range of reaction progress, variable across the flame front (se later in Figure 9 ). Moreover, the differences in molecular diffusion contribution (espe cially the flame tangential diffusion rate) between the premixed flames and MILD com bustion cases also contribute to this (see Figure 9 later in this paper).
The inverse of the peak value of the mean value of |∇ | conditioned upon pro  both the laminar and turbulent conditions in conventional premixed flame cases, suggest that the flame thickness does not change significantly under turbulent conditions. However, an increase in ⁄ increases the standard deviation of |∇ |, since the range of fluctuation increases with increasing turbulence intensity. It is worthwhile to analyse the strain rates, which affect the evolution of the SDF (see Equation (9)), to explain the observed behaviours of |∇ | in both premixed flame and MILD combustion cases.

Mean Behaviour of Fluid-Dynamic Strain Rates
The profiles of the mean values of { , and ∇. ⃗} × / conditioned upon are presented in Figure 6. The dilatation rate ∇. ⃗ assumes positive values because of the The inverse of the peak value of the mean value of |∇c| conditioned upon c provides the measure of the flame thickness (i.e., δ ∼ 1/max |∇c| c where . . . c is the mean value conditioned upon c) [32]. A comparison between the profiles of the MILD combustion and premixed flame cases reveals that the mean values of |∇c| × δ th in MILD combustion cases remain smaller than those in the premixed flame cases, which is indicative of the thickened reaction zone separating reactants and products in contrast to a thin reaction zone in the premixed flame cases. This is also consistent with the observations made from Figures 2-4. It can also be seen from Figure 5 that the maximum value of |∇c| c in the MILD combustion cases increases slightly with increasing u /S L but the differences between the mean values at different turbulence intensities remain much smaller than the standard deviation suggesting that differences in |∇c| c are not significant enough to imply any change in the thickness of the reaction zone under turbulent conditions. However, it can be seen from Figure 5 that the peak mean value of |∇c| c in the MILD combustion cases remains significantly smaller than the corresponding value in the 1D premixed flame case, indicating that the flame thickens considerably under turbulent conditions for the MILD combustion cases. The interactions of the flame surfaces in the MILD combustion cases, as shown in Figures 2-4, contribute to this thickening of the flame under turbulent conditions. By contrast, the peak mean values of |∇c| c for both laminar and turbulent conditions are almost identical to each other for the conventional premixed flame cases. Figure 5 also shows that the peak value of |∇c| c for the turbulent premixed flame cases do not change significantly with a change in u /S L , but a small decrease in the peak value of |∇c| c is discerned with an increase in u /S L . However, the standard deviations of |∇c| in the premixed flame cases also remain greater than the differences in mean values and, therefore, the differences in the mean values of |∇c| do not signify a change in the mean flame thickness in response to the variation in turbulence intensity within the parameter range considered here. Moreover, comparable peak mean values of |∇c| c , for both the laminar and turbulent conditions in conventional premixed flame cases, suggest that the flame thickness does not change significantly under turbulent conditions. However, an increase in u /S L increases the standard deviation of |∇c|, since the range of fluctuation increases with increasing turbulence intensity. It is worthwhile to analyse the strain rates, which affect the evolution of the SDF (see Equation (9)), to explain the observed behaviours of |∇c| in both premixed flame and MILD combustion cases.

Mean Behaviour of Fluid-Dynamic Strain Rates
The profiles of the mean values of {a N , a T and ∇· → u } × δ th /S L conditioned upon c are presented in Figure 6. The dilatation rate ∇· → u assumes positive values because of the Energies 2021, 14, 7677 13 of 28 thermal expansion due to heat release for both the premixed and MILD flame cases but the mean values of ∇· → u × δ th /S L , in the MILD combustion cases, remain negligibly small in comparison with the premixed flame cases. This suggests that the thermal expansion effects are relatively weak in the MILD combustion cases compared with conventional premixed flames. In the premixed flames, the peak mean value of ∇· → u × δ th /S L is obtained in the middle of the flame with a skewness towards the burned gas side (i.e., c ≈ 0.65), whereas the mean value of ∇· → u × δ th /S L remains vanishingly small for all c values in the MILD combustion cases.  It can be seen from Figure 6 that the mean value of a N remains negative towards the unburned gas side of the flame, while positive values are obtained for larger values of c in the premixed flame cases. The normal strain rate a N can be expressed as [41][42][43][44]62]: a N = (e α cos 2 θ α + e β cos 2 θ β + e γ cos 2 θ γ ) (12) where e α , e β and e γ are the most extensive, intermediate and the most compressive principal strain rates, respectively and θ α , θ β and θ γ are the angles between ∇c with eigendirections corresponding to e α , e β and e γ , respectively. The PDFs of |cos θ α |, cos θ β and |cos θ γ | for different c− isolevels across the flame for both MILD combustion and premixed flame cases considered here are shown in Figure 7. A high probability of |cos θ i | = 1.0 (for i = α, β, γ) indicates a preferential collinear alignment between ∇c with the eigenvector associated with e i (for i = α, β, γ). Figure 7 shows that the reaction progress variable gradient, ∇c, preferentially aligns collinearly with the most compressive principal strain rate eigendirection throughout the flame in the MILD combustion cases, similarly to passive scalar mixing [63][64][65][66][67][68][69][70][71][72], and a preferential collinear alignment between ∇c and the eigenvector associated with e γ , in the premixed flame cases, are observed towards the unburned and burned gas sides of the flame where the effects of heat release are weak. However, it can be seen from Figure 7 that ∇c preferentially collinearly aligns with the most extensive principal strain rate eigendirection in the region of the flame, where the effects of thermal expansion are strong (i.e., where dilatation rate assumes large values) in the premixed flame cases. In the case of passive scalar mixing, ∇c preferentially collinearly aligns with the most compressive principal strain rate eigendirection [63][64][65][66][67][68][69][70][71][72], whereas, in premixed flames, ∇c may align preferentially with the most extensive principal strain rate eigendirection if the strain rate induced by flame normal acceleration overcomes turbulent straining [62,72]. The strain rate induced by flame normal acceleration a chem can be taken to scale with τS L /δ th (where τ = ρ u /ρ b − 1 is the heat release parameter with ρ u and ρ b are the unburned gas density and the burned gas density, respectively) [62], whereas large-scale turbulent straining is given by a turb ∼ u /l. This suggests that ∇c aligns preferentially with the most extensive principal strain rate eigendirection for τS L /δ th u /l, which is characterised by τDa 1. By contrast, ∇c aligns preferentially with the most compressive principal strain rate eigendirection for τS L /δ th u /l, which is characterised by τDa 1. For MILD combustion cases, τDa remains small due to small values of τ and thus a preferential alignment of ∇c with the most compressive principal strain rate eigendirection is observed in these cases, leading to predominantly negative values of a N , according to Equation (12). This tendency strengthens with an increase in u /S L , as a result of the reduction in Da, and therefore, the mean value of a N assumes negative values of higher magnitudes for higher values of u /S L (see Figure 6). For the premixed flame cases considered here, remains greater than 1.0 (i.e., > 1) and thus, a preferential alignment of ∇ with the most extensive principal strain rate eigendirection is observed within a major part of the flame in these cases, leading to predominantly positive values of , according to Equation (12). This is valid in the re-  For the premixed flame cases considered here, τDa remains greater than 1.0 (i.e., τDa > 1) and thus, a preferential alignment of ∇c with the most extensive principal strain rate eigendirection is observed within a major part of the flame in these cases, leading to predominantly positive values of a N , according to Equation (12). This is valid in the regions of heat release, but an increase in u /S L leads to a drop in Da, which reduces the extent of the preferential alignment of ∇c with the most extensive principal strain rate eigendirection leading to a marginal decrease in the positive mean values of a N with an increase in u /S L in the premixed flame cases.
As the dilatation rate ∇· → u is vanishingly small in the MILD combustion cases, the predominantly negative mean values of a N lead to positive mean values of tangential strain rate a T = ∇· → u − a N in these cases. In the premixed flames considered here, the mean value of ∇· → u remains greater than that of a N , which leads to positive mean values of a T throughout the flame. As the mean value of ∇· → u is primarily determined by thermochemistry, the decrease in mean values of a N with increasing u /S L leads to an increase in the mean value of tangential strain rate a T = ∇· → u − a N . The positive mean value of a T is indicative of flame area generation under fluid-dynamic straining in all the cases considered here.

Mean Behaviour of Flame Displacement Speed
It is evident from Equations (9) and (10)   It can be seen from Figure 9 that the mean value of . ω c × δ th /ρ u S L peaks towards the burned gas side of the flame for both MILD combustion and conventional premixed flames. However, it can be seen from Figures 8 and 9 that the mean value of . ω c × δ th /ρ u S L , and, accordingly, of S r /S L , assume non-negligible values in MILD combustion for a broader range of c than in the conventional premixed flame cases. The unburned gas temperature in MILD combustion cases is above the autoignition temperature, which gives rise to non-negligible reaction rates . ω c for the unburned gas side of the flame-front, whereas the   The mean value of the flame normal diffusion rate ⃗ • ∇ ⃗ • ∇ is positive towards the unburned gas side but becomes negative on the burned gas side of the flame for both the MILD combustion and conventional premixed flame cases, and therefore the mean normal diffusion component of displacement speed, , also assumes positive (negative) mean values for the unburned (burned) gas side of the flame. The magnitude of the mean contribution of −2 |∇ | remains small in comparison to the leading order mean contributions of ̇ and ⃗ • ∇ ⃗ • ∇ for all cases considered in this analysis (see Figure 9). This is due to the statistically planar nature of these cases. Accordingly, the mean value of remains negligible in comparison to the mean values of and for both premixed and MILD combustion cases. Thus, and are the major contributors to the mean values of for all cases considered here. However, there are differences in the mean behaviours of and between premixed flame and MILD combustion cases. N·∇c) for all cases considered in this analysis (see Figure 9). This is due to the statistically planar nature of these cases. Accordingly, the mean value of S t remains negligible in comparison to the mean values of S r and S n for both premixed and MILD combustion cases. Thus, S r and S n are the major contributors to the mean values of S d for all cases considered here. However, there are differences in the mean behaviours of S r and S n between premixed flame and MILD combustion cases. Figure 9 indicates the mean values of . ω c and → N·∇(ρD → N·∇c) (and also ∇(ρD∇c)) remain of the same order of magnitude but with different signs in the reaction zone, and, accordingly, the mean values of S r and S n are of the same order in the reaction zone both in the premixed flame and MILD combustion cases. It is important to note that δ th /ρ u S L , in the MILD combustion cases, is a factor of 2.0 smaller than the premixed flame cases considered here. Therefore, the magnitudes of . ω c × δ th /ρ u S L and → N·∇(ρD → N·∇c) × δ th /ρ u S L ) between MILD combustion and conventional premixed flame cases in Figure 9. However, the relative magnitudes of mean S r and S n are such that the mean S d increases with increasing c in turbulent premixed flames, whereas the positive contribution of . ω c is partially nullified by the negative contribution of → N·∇(ρD → N·∇c) in the MILD combustion cases in such a manner that the mean displacement speed S d decreases from the unburned gas side to the middle of the flame and subsequently increases mildly with c for c > 0.6 (see Figure 8). Moreover, it can be seen from Figure 8 that the mean S d is not affected by turbulence intensity in the premixed flame cases. The statistical behaviour of S d and the differences in its behaviour between premixed flames and MILD combustion cases, in turn, affect the strain rates induced by flame propagation. Although the mean values of S d /S L in the MILD combustion cases are comparable for the conventional premixed flame and MILD combustion cases, the higher value of S L in the MILD combustion cases implies higher magnitude of S d in these cases in comparison to that in the conventional premixed flame cases.

Mean Behaviour of Strain Rates Induced by Flame Propagation
The profiles of the normalised mean values of the additional flame normal strain rate due to flame propagation (N j ∂S d /∂x j × δ th / S L ) and its components (i.e., N j ∂(S r + S n )/∂x j × δ th / S L and N j ∂S t /∂x j × δ th / S L ) conditioned upon c are shown in Figure 10. It is worth noting that δ th values for the conventional premixed flame and MILD combustion cases remain comparable [20,21] but S L values are significantly different, so the chemical timescale δ th / S L in the MILD combustion cases remains almost 1/10 of the corresponding value in the premixed flame cases. This suggests that N j ∂S d /∂x j and its components assume larger magnitudes in the MILD combustion cases than in the conventional premixed flames. The higher magnitude of displacement speed in the MILD combustion cases and comparable values of flame thickness contribute to larger magnitudes of mean values of N j ∂S d /∂x j and its components in the MILD combustion cases than those in the premixed flame cases. It can further be seen from Figure 10 that the mean value of N j ∂S d /∂x j remains positive towards the unburned gas side but assumes negative values for a major part of the flame front in the premixed flame cases because the mean value of S d increases with increasing c, which is consistent with previous studies [41][42][43]. The mean contribution of N j ∂S d /∂x j assumes positive values throughout the flame in the MILD combustion cases because of the qualitatively different variation of S d in comparison to the corresponding premixed flame cases (see Figure 5). It can also be seen from Figure 10 that the mean value of N j ∂S d /∂x j is principally driven by the mean contribution of N j ∂(S r + S n )/∂x j in the premixed flame cases. However, it is worth noting that the local curvature variations give rise to non-negligible positive mean value of N j ∂S t /∂x j for the premixed flame cases. In MILD combustion cases, the contribution of N j ∂S t /∂x j assumes positive values across the whole range of c values and is of the same order of magnitude as N j ∂(S r + S n )/∂x j . As the extent of flame curvature variation increases with increasing turbulence intensity, the magnitude of N j ∂S t /∂x j = −2N j ∂(D c κ m )/∂x j increases with increasing u /S L . The diffusivity, D c , increase from the unburned to the burned gas side of premixed flames acts to promote negative values of N j ∂(D c κ m )/∂x j , which eventually leads to positive mean values of N j ∂S t /∂x j = −2N j ∂(D c κ m )/∂x j . The effects of the increase in D c remain weak in the MILD combustion cases and a decrease in flame curvature towards the burned gas side, due to the weakening of turbulence, induces predominantly positive values of N j ∂(D c κ m )/∂x j , which eventually leads to a decrease in the mean values of N j ∂S t /∂x j = −2N j ∂(D c κ m )/∂x j towards the burned gas side of these cases. The mean values of N j ∂(S r + S n )/∂x j for the turbulent premixed flame cases remain comparable for all the turbulence intensities considered in this analysis. However, the higher positive mean value of N j ∂S t /∂x j in the premixed flame case with initial u /S L = 8.0 gives rise to a smaller negative mean value of N j ∂S d /∂x j than in cases with initial u /S L = 4.0 and 6.0. As the mean value of S d is significantly affected by turbulence intensity in MILD combustion cases, the mean value of N j ∂S d /∂x j also shows significant u /S L dependence and the mean value of positive contribution of N j ∂S d /∂x j increases with increasing u /S L . According to Equation (9), a negative (positive) N j ∂S d /∂x j tends to promote flame thinning (flame thickening) and therefore the mean contribution of the additional normal strain rate N j ∂S d /∂x j to the SDF evolution is different in the MILD combustion cases in comparison with the conventional premixed flame. The relative magnitudes of the mean values of a N and N j ∂S d /∂x j determine the mean behaviour of the effective normal strain rate a / tends to promote flame thinning (flame thickening) and therefore the mean contribution of the additional normal strain rate / to the SDF evolution is different in the MILD combustion cases in comparison with the conventional premixed flame. The relative magnitudes of the mean values of and / determine the mean behaviour of the effective normal strain rate . The profiles of the normalised mean values of the additional flame tangential strain rate due to flame propagation (i.e., alternatively curvature stretch rate given by 2 × / ) and its components conditioned upon are shown in Figure 11. It can be seen from Figure 11 that the mean value of 2 remains negative for both conventional premixed flames, and MILD combustion cases considered here. In all cases the negative definite value of 2 = −4 is the major contributor to the mean value of curvature stretch 2 . This is consistent with the scaling arguments by Peters [53], suggesting that the stretch rate induced by the tangential diffusion component of displacement speed is the dominant component of 2 for ≫ 1 cases, as in the cases considered here. The mean value of 2( + ) increases from the unburned gas side to the burned gas side for both conventional premixed flames and MILD combustion cases considered here, but this contribution is dominated by the negative contribution of 2 = −4 . The negative mean values of 2 indicate that the curvature stretch in the cases considered here act to reduce the flame area generation under turbulence for all cases considered here. However, there are qualitative differences in the variations of the mean values of Figure 10. Profiles of the normalised mean values of the additional flame normal strain rate due to flame propagation (N j ∂S d /∂x j × δ th / S L ) and its components (i.e., N j ∂(S r + S n )/∂x j × δ th / S L and N j ∂S t /∂x j × δ th / S L ) conditioned upon c for u /S L = 4.0 (black) 6.0 (red) and 8.0 (blue) for statistically planar turbulent premixed flames (first column) and homogeneous MILD combustion cases (second column). Note the difference in scales between premixed flame and MILD combustion cases.
The profiles of the normalised mean values of the additional flame tangential strain rate due to flame propagation (i.e., alternatively curvature stretch rate given by 2S d κ m × δ th / S L ) and its components conditioned upon c are shown in Figure 11. It can be seen from Figure 11 that the mean value of 2S d κ m remains negative for both conventional premixed flames, and MILD combustion cases considered here. In all cases the negative definite value of 2S t κ m = −4D c κ 2 m is the major contributor to the mean value of curvature stretch 2S d κ m . This is consistent with the scaling arguments by Peters [53], suggesting that the stretch rate induced by the tangential diffusion component of displacement speed is the dominant component of 2S d κ m for Ka 1 cases, as in the cases considered here. The mean value of 2(S r + S n )κ m increases from the unburned gas side to the burned gas side for both conventional premixed flames and MILD combustion cases considered here, but this contribution is dominated by the negative contribution of 2S t κ m = −4D c κ 2 m . The negative mean values of 2S d κ m indicate that the curvature stretch in the cases considered here act to reduce the flame area generation under turbulence for all cases considered here. However, there are qualitative differences in the variations of the mean values of 2S d κ m between the MILD combustion and conventional premixed flame cases. In the conventional premixed flame cases, the mean value of 2S d κ m does not alter significantly across a major part of the flame-front, whereas the magnitude of the negative mean value of 2S d κ m in the MILD combustion cases decreases from the unburned to the burned gas side. As the extent of flame wrinkling increases with increasing turbulence intensity, a larger range of curvature is obtained for higher values of turbulence intensity, which acts to increase the magnitude of 2S t κ m = −4D c κ 2 m with an increase in u /S L . As δ th / S L in the MILD combustion cases remains almost 1/10 of the corresponding value in the premixed flame cases, the magnitudes of 2S t κ m and 2S d κ m , in the MILD combustion cases, are greater than those in the conventional premixed flames. The higher magnitude of displacement speed and greater extent of flame wrinkling (i.e., larger spread of curvature values) in the MILD combustion cases contribute to larger magnitudes of 2S t κ m and 2S d κ m in these cases than in the premixed flame cases. mixed flame cases, the magnitudes of 2 and 2 , in the MILD combustion cases, are greater than those in the conventional premixed flames. The higher magnitude of displacement speed and greater extent of flame wrinkling (i.e., larger spread of curvature values) in the MILD combustion cases contribute to larger magnitudes of 2 and 2 in these cases than in the premixed flame cases. Figure 11. Profiles of the normalised mean values of the additional flame tangential strain rate due to flame propagation (2 × / ) and its components (i.e., 2( + ) × / and 2 × / ) conditioned upon for ⁄ = 4.0 (black) 6.0 (red) and 8.0 (blue) for statistically planar turbulent premixed flames (first column) and homogeneous MILD combustion cases (second column).

Mean Behaviour of Effective Strain Rates
The profiles of the normalised mean values of the effective normal rate × / conditioned upon are shown in Figure 12. According to Equation (9), a positive value of at low values of indicates that there is a propensity for a decrease in |∇ |, while a negative effective normal strain rate at high values of gives rise to an increase in |∇ |. For the premixed flames considered here, the mean values of remain negative for 0.3 < < 0.9 for initial ⁄ = 4.0 and 6.0, whereas, in the case of initial ⁄ = 8.0, the mean value of remains mostly positive, but its magnitude decreases in the region where ̇ assumes high values (see Figures 9 and 12). An increase in / increases the positive mean value of towards the preheat zone for the premixed flame cases, which acts to reduce the SDF in this region (see Figure 5) and therefore acts to thicken the preheat zone in comparison with the unstretched laminar flame. The positive mean values of , in the premixed flame cases, are the consequence of the mean value of dominating over the negative mean values of / . The smaller magnitudes of negative mean values of / , for the initial ⁄ = 8.0 premixed flame case, dominates over the smaller positive mean values of to yield higher mean values of than in the premixed flames with initial ⁄ = 4.0 and 6.0. Figure 11. Profiles of the normalised mean values of the additional flame tangential strain rate due to flame propagation (2S d κ m × δ th / S L ) and its components (i.e., 2(S r + S n )κ m × δ th / S L and 2S t κ m × δ th / S L ) conditioned upon c for u /S L = 4.0 (black) 6.0 (red) and 8.0 (blue) for statistically planar turbulent premixed flames (first column) and homogeneous MILD combustion cases (second column).

Mean Behaviour of Effective Strain Rates
The profiles of the normalised mean values of the effective normal rate a e f f N × δ th /S L conditioned upon c are shown in Figure 12. According to Equation (9) Figure 12. It can be seen from Figure 12 that the mean values of a e f f N |∇c| mildly increases with increasing u /S L for turbulent premixed flames considered here, whereas MILD combustion cases the mean values of a e f f N |∇c| increases more clearly with increasing u /S L . As δ th /S L in the MILD combustion cases is smaller than that in the premixed cases considered here, Figure 12 indicates that the mean values of a e f f N are higher in magnitude in the MILD combustion cases than in the premixed combustion cases. The The mean values of for the MILD combustion cases remain positive throughout the flame and its value increases with increasing ⁄ . The high mean positive value of is indicative of flame thickening under turbulence in comparison to the corresponding laminar flame in the MILD combustion cases, as shown in Figure 5. It is useful to consider that the left-hand side of Equation (9) indicates the fractional change in the SDF Equation (9) can be rewritten as: |∇ |/ = − |∇ |. The mean values of |∇ | conditional upon for both the premixed flame and MILD combustion cases considered here are also presented in Figure 12. It can be seen from Figure 12 that the mean values of |∇ | mildly increases with increasing ⁄ for turbulent premixed flames considered here, whereas MILD combustion cases the mean values of |∇ | increases more clearly with increasing ⁄ . As / in the MILD combustion cases is smaller than that in the premixed cases considered here, Figure 12 indicates that the mean values of are higher in magnitude in the MILD combustion cases than in the premixed combustion cases. The higher magnitude of / , in the MILD combustion cases contribute to the larger mean values of in these cases than those in the premixed flame cases. By the same token, the mean values of |∇ | in the MILD combustion cases are greater than those for turbulent premixed flames considered in this analysis because values are comparable for both MILD combustion and premixed flame cases [21].
The total derivative |∇ | | | ⁄ = |∇ | [ |∇ | ⁄ + ( + ) |∇ | ⁄ ] provides the temporal change of the SDF in a reference frame attached to the flame. However, the flames considered here are not stationary in an instantaneous sense (i.e., ( + ) ≠ 0) and thus |∇ | ⁄ ≠ |∇ | ⁄ . This can be substantiated from the non-zero mean values of |∇ | ( + ) |∇ | ⁄ conditioned upon , which are presented for both premixed and MILD combustion cases in Figure 13. Figure 13 shows that the mean value of |∇ | ( + ) |∇ | ⁄ remains negative towards the unburned gas side but becomes positive towards the burned gas side of the flame front for all the premixed flame cases considered here. The magnitudes of the mean value of |∇ | ( + ) |∇ | ⁄ remain comparable for all turbulence intensities for the premixed combustion cases. By contrast, the mean value of |∇ | ( + ) |∇ | ⁄ remains negative throughout the The total derivative |∇c| −1 d|∇c|/dt = |∇c| −1 [∂|∇c|/∂t + (u j + S d N j )∂|∇c|/∂x j ] provides the temporal change of the SDF in a reference frame attached to the flame. However, the flames considered here are not stationary in an instantaneous sense (i.e., (u j + S d N j ) = 0) and thus d|∇c|/dt = ∂|∇c|/∂t. This can be substantiated from the nonzero mean values of |∇c| −1 (u j + S d N j )∂|∇c|/∂x j conditioned upon c, which are presented for both premixed and MILD combustion cases in Figure 13. Figure 13 shows that the mean value of |∇c| −1 (u j + S d N j )∂|∇c|/∂x j remains negative towards the unburned gas side but becomes positive towards the burned gas side of the flame front for all the premixed flame cases considered here. The magnitudes of the mean value of |∇c| −1 (u j + S d N j )∂|∇c|/∂x j remain comparable for all turbulence intensities for the premixed combustion cases. By contrast, the mean value of |∇c| −1 (u j + S d N j )∂|∇c|/∂x j remains negative throughout the flame front but the magnitude increases with increasing u /S L . It has been demonstrated by Dopazo et al. [42] that the relative magnitudes of |∇c| −1 (u j + S d N j )∂|∇c|/∂x j and a e f f N determine the behaviour of |∇c| −1 ∂|∇c|/∂t = −|∇c| −1 (u j + S d N j )∂|∇c|/∂x j − a e f f N . The mean and standard deviations of |∇c| −1 ∂|∇c|/∂t × δ th /S L conditioned upon c are shown in Figure 14 for both premixed flame and MILD combustion cases. It can be seen from Figure 14 that the magnitude of the mean value of |∇c| −1 ∂|∇c|/∂t remains much smaller than the magnitudes of mean values of |∇c| −1 (u j + S d N j )∂|∇c|/∂x j and a e f f N for both conventional premixed combustion and MILD combustion cases. Moreover, the standard deviations of |∇c| −1 ∂|∇c|/∂t remain much greater than the differences in the mean values for different turbulence intensities. Furthermore, it has been found that the mean value of |∇c| −1 ∂|∇c|/∂t vanishes at the value of c where the peak value of |∇c| c is obtained for both conventional premixed flame and MILD combustion cases for all turbulence intensities considered here. This is consistent with the findings from Figure 5, which suggest that the peak value of |∇c| c is not significantly affected by the variation of u /S L for both conventional turbulent premixed flame and turbulent MILD combustion cases.
Energies 2021, 14, 7677 21 of 28 mean value of |∇ | |∇ | ⁄ vanishes at the value of where the peak value of 〈|∇ |〉 is obtained for both conventional premixed flame and MILD combustion cases for all turbulence intensities considered here. This is consistent with the findings from Figure 5, which suggest that the peak value of 〈|∇ |〉 is not significantly affected by the variation of / for both conventional turbulent premixed flame and turbulent MILD combustion cases.  mean values for different turbulence intensities. Furthermore, it has been found that the mean value of |∇ | |∇ | ⁄ vanishes at the value of where the peak value of 〈|∇ |〉 is obtained for both conventional premixed flame and MILD combustion cases for all turbulence intensities considered here. This is consistent with the findings from Figure 5, which suggest that the peak value of 〈|∇ |〉 is not significantly affected by the variation of / for both conventional turbulent premixed flame and turbulent MILD combustion cases. The mean and standard deviations of |∇ | ⁄ × / conditioned upon are also shown in Figure 14 for both premixed flame and MILD combustion cases, which reveal that the mean value of |∇ | ⁄ remains negative for the MILD combustion cases, whereas the mean value of this quantity remains weakly positive for the major part of the The mean and standard deviations of ∂|∇c|/∂t × δ 2 th /S L conditioned upon c are also shown in Figure 14 for both premixed flame and MILD combustion cases, which reveal that the mean value of ∂|∇c|/∂t remains negative for the MILD combustion cases, whereas the mean value of this quantity remains weakly positive for the major part of the flame front before assuming negative values towards the burned gas side for the premixed flame cases. Once again it has been found that the standard deviations of ∂|∇c|/∂t remained much greater than the differences in the mean values for different turbulence intensities and, thus, the change in turbulence intensity does not alter the mean behaviour of ∂|∇c|/∂t significantly for the range of u /S L considered here. However, the standard deviation of ∂|∇c|/∂t increases with increasing u /S L but the negative mean value of ∂|∇c|/∂t in the turbulent MILD combustion cases are consistent with the reduction in the peak mean value of |∇c| in these cases in comparison with the 1D unstretched laminar flame results. It is worth noting that the normalization factor δ 2 th /S L for ∂|∇c|/∂t for the MILD cases is 0.1 times of that of the conventional premixed flame cases and therefore the magnitudes of ∂|∇c|/∂t in the premixed flame cases are smaller than in the corresponding MILD combustion cases. These small magnitudes of ∂|∇c|/∂t, in the premixed flame cases, are consistent with the finding that the peak value of |∇c| c for turbulent premixed flames is comparable to the corresponding 1D unstretched laminar premixed flame solution.
The profiles of the normalised mean values of the effective tangential rate a e f f T × δ th /S L conditioned upon c are shown in Figure 15. It can be seen from Figure 15 Figure 15. Profiles of the normalised mean values and standard deviations of the effective tangential strain rate ( × / ) conditioned upon for ⁄ = 4.0 (black) 6.0 (red) and 8.0 (blue) for statistically planar turbulent premixed flames (first column) and homogeneous MILD combustion cases (second column). Note the difference in scales between the premixed flame and MILD combustion cases.
It has already been mentioned that / , in the MILD combustion cases, is smaller than that in the premixed cases considered here, and thus the mean values of × / presented in Figure 15 imply that the mean values of assume negative values of larger magnitude in the MILD combustion cases than in the premixed combustion cases. This is a consequence of a larger mean contribution of 2 = −4 in the MILD combustion cases than in the premixed cases considered here (see Figure 11). It can also be seen from Figure 15 that the standard deviations of are much greater than the differences in mean values for different initial ⁄ , but the magnitude of the standard deviation increases with an increase in ⁄ . However, the negative mean values of It has already been mentioned that δ th /S L , in the MILD combustion cases, is smaller than that in the premixed cases considered here, and thus the mean values of a e f f T × δ th /S L presented in Figure 15 imply that the mean values of a e f f T assume negative values of larger magnitude in the MILD combustion cases than in the premixed combustion cases. This is a consequence of a larger mean contribution of 2S t κ m = −4D c κ 2 m in the MILD combustion cases than in the premixed cases considered here (see Figure 11). It can also be seen from Figure 15 that the standard deviations of a e f f T are much greater than the differences in mean values for different initial u /S L , but the magnitude of the standard deviation increases with an increase in u /S L . However, the negative mean values of a e f f T imply that the flame area destruction due to flame self-interaction is more prevalent in the MILD combustion cases than in the premixed flame cases.

Modelling Implications
In order to explain the relevance of the strain rate statistics discussed in previous sections on turbulent combustion modelling, it is useful to write the transport equation of |∇c| by expanding Equation (9) in the following manner: which can alternatively be written as: ∂|∇c| ∂t The transport equation of the generalised FSD (i.e., Σ gen = |∇c|) can be obtained by Reynolds averaging/LES filtering Equation (14).
Multiplying Equation (13a) by 2 |∇c| yields [39,41,42]: Algebraic manipulation of Equation (15) provides the transport equation of scalar dissipation rate (SDR) N c = D c |∇c| 2 [41,43,44]: Equations (13), (14) and (16) show that a N , a T , N j ∂S d /∂x j , 2S d κ m , a e f f N and a e f f T play key roles in both the FSD and SDR evolutions. The correlation coefficients of |∇c| and S d with tangential strain rate a T and curvature κ m are listed in Table 3 for both conventional premixed flame and MILD combustion cases considered here. It can be seen from Table 3 that the correlation coefficients between |∇c| and a T , and between S d and κ m have the same sign for both conventional premixed flame and MILD combustion cases but the correlations are weaker for the MILD combustion cases. Moreover, the correlations between S d and a T , and between |∇c| and κ m in the MILD combustion cases are significantly weaker than the corresponding values obtained for the premixed flame cases. Although on some occasions, the signs of S d − a T and |∇c| − κ m in the MILD combustion cases are different to the premixed flame cases, these correlation coefficient magnitudes are so small for the MILD combustion cases that these differences in signs of correlation coefficients are of limited significance. The qualitative similarities in the SDF statistics and the strain rates affecting its evolution, and local curvature and strain rate dependences of S d and |∇c|, respectively between premixed flame and MILD combustion cases indicate that the models for FSD [74][75][76] and SDR [29,[77][78][79][80] transports, which were proposed for the thin-reaction-Energies 2021, 14, 7677 24 of 28 zones regime combustion, may have the potential to be extended for the modelling of EGR-type MILD combustion of homogeneous-mixtures.

Conclusions
The statistical behaviours of the magnitude of the reaction progress variable gradient, known as the surface density function (SDF), between premixed flame and EGR-type homogeneous-mixture MILD methane-air combustion cases under different turbulence intensities have been compared using three-dimensional DNS, employing a skeletal chemical mechanism consisting of 16 species and 25 reactions proposed by Smooke and Giovangigli [45]. The turbulent MILD combustion cases showed significantly smaller mean values of the SDF than that in the corresponding unstretched laminar flame condition, whereas the mean SDF values for turbulent premixed flames were found to be comparable to the corresponding unstretched 1D laminar flame values. It was observed that the peak mean value of the SDF did not change significantly in response to the changes in turbulence intensity for both MILD and premixed flame cases in the range considered in this study. This behaviour has been explained in terms of the strain rates induced by fluid motion and the ones arising from flame propagation. It was found that the effects of dilatation rate were relatively weaker under MILD combustion conditions and ∇c exhibited greater propensity to align with the most compressive principal strain rate eigendirection, whereas, in the case of the premixed flames, ∇c tended to align with the most extensive strain rate eigenvector in the regions of high dilatation rate caused by intense heat release. As a consequence of the variations in the alignment of ∇c with the strain rate eigenvectors, a predominantly positive mean value of normal strain rate was observed in the premixed flames, which decreased in magnitude with increasing turbulence intensity. By contrast, the mean normal strain rate remained negative in the case of MILD combustion, and its magnitude increased with increasing turbulence intensity.
The mean behaviour of the displacement speed and its components have been investigated in detail, and it was found that the reaction contribution assumes non-negligible values in MILD combustion cases for a broader range of c compared with the conventional premixed flames. Moreover, the mean displacement speed increased from the unburned gas side to the burned gas side in conventional premixed flames, whereas, in the MILD cases, the mean displacement speed was found to decrease from the unburned gas side to the middle of the flame before increasing weakly towards the burned gas side. These differences in the mean displacement speed and its components between the premixed flames and MILD combustion cases gave rise to significant differences in the mean behaviour of the normal strain rate induced by the flame propagation and the effective normal strain rate. The aforementioned statistics have been utilized to explain the observed behaviour in the SDF evolution and its response to the variation of turbulence intensity in conventional premixed flames and MILD combustion cases. It was found that the tangential fluid-dynamic strain rate assumed positive mean values but was overcome by the negative mean values of curvature stretch rate for both the premixed flames and MILD combustion cases considered here. This yielded negative mean values of stretch rate (or effective tangential strain rate) for both the premixed flame and MILD combustion cases and its magnitude increased with increasing turbulence intensity. This behaviour is explained in terms of the variation in the displacement speed and its dependence on curvature. The local curvature and strain rate dependences of the SDF and displacement speed were found to be mostly qualitatively similar for both the premixed and MILD combustion cases. Overall, the findings from this work suggest that the curvature dependence of displacement speed and the scalar gradient alignment with local principal strain rate eigendirections are important in controlling the behaviour of SDF which need to be explicitly accounted for in the extended FSD/SDR based modelling methodology for the EGR-type homogeneous-mixture MILD combustion. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.