DNS of MILD combustion with mixture fraction variations

Direct numerical simulations of Moderate or Intense Low-oxygen Dilution combustion inside a cubical domain are performed. The computational domain is specified with inflow and outflow boundary conditions in one direction and periodic conditions in the other two directions. The inflowing mixture is constructed carefully in a preprocessing step and has spatially varying mixture fraction and reaction progress variable fields. Thus, this mixture includes a range of thermo-chemical states for a given mixture fraction value. The combustion kinetics is modelled using a 58-step skeletal mechanism including a chemiluminescent species, OH ∗, for methane–air combustion. The study of reaction zone structures in the physical and mixture fraction spaces shows the presence of ignition fronts, lean and rich premixed flames and non-premixed combustion. These three modes of combustion are observed without the typical triple-flame structure and this results from the spatio-temporally varying mixture fraction field undergoing turbulent mixing and reaction. The flame index and its pdf are analysed to estimate the fractional contributions from these combustion modes to the total heat release rate. The lean premixed mode is observed to be quite dominant and contribution of non-premixed mode increased from about 11% to 20% when the mean oxygen mole fraction in the inflowing mixture is reduced from about 2.7% to 1.6%. Also, the non-premixed contribution increases if one decreases the integral length scale of the mixture fraction field. All of these results and observations are explained on physical basis. © 2017 The Authors. Published by Elsevier Inc. on behalf of The Combustion Institute. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Improvements in energy efficiency and pollutants emission are ever increasing demands on combustion devices and these demands could be met using Moderate or Intense, Low-oxygen Dilution (MILD) combustion concept [1,2] . The key aspect of this concept involves preheating and dilution of reactant mixture which are achieved by recirculating combustion products. This recirculation can either be external or internal to combustion chamber and this arrangement typically gives reactants (or air) temperature, T r , higher than the autoignition temperature, T ign , for a given thermochemical condition. Also the temperature rise from combustion, T = T p − T r , is lower than T ign [2] . As a result of these characteristics, energy efficiency is improved by recovering exhaust heat and the emission is reduced substantially because of homogeneous heat release with low peak temperature resulting from dilution [1][2][3] . * Corresponding author.
Simplified configurations have been used in past studies to improve our understanding of MILD combustion [1][2][3] . The homogeneous temperature field combined with low gradients of species and temperature observed in MILD combustion experiments [4,5] led Cavaliere and co-workers [2,6,7] to use the perfectly stirred reactor (PSR) as a representative model for MILD combustion. Other studies used the counterflow configuration to investigate characteristics of this combustion in laminar flow conditions [8][9][10] . The studies using the PSR configuration suggested that these combustion attributes were compatible with a twostage oxidation process and it occurred for a particular combination of inlet temperature, dilution level and mixture residence time highlighting a need to understand autoignition of MILD mixtures, which was explored in later studies using laminar counterflow configurations [8][9][10] . The suitability of these configurations to turbulent MILD combustion is, however, an open question. Indeed, MILD combustion has been noted to be a novel combustion mode which cannot be categorised using conventional classification such as premixed and non-premixed modes [2] . It has indeed been shown that turbulent MILD combustion has distinctive features of its own while sharing some common attributes of conventional combustion [11,12] .
Direct Numerical Simulation (DNS) has been used in the past to address some of the above questions and to shed light on fundamental processes of MILD combustion. Past DNS studies of an auto-igniting mixing layer between methane and hot diluted coflow highlighted the role of the most reactive mixture fraction and preferential diffusion effects in the formation of reaction zones [13,14] . The broadening and frequent interactions of reaction zones in MILD combustion established using recirculation (internal to combustor) of exhaust gas were highlighted in other DNS studies [11,12,15,16] . Furthermore, the suitability of PSR based tabulated chemistry approach to model MILD combustion was highlighted through detailed analysis of the MILD combustion DNS [17,18] . The flamelet methodology may also become suitable if the effects of flame-flame interactions are included in this approach [15] . DNS studies of Minamoto et al., which are the first of their kind, considered a single mixture fraction while actual MILD combustion would involve spatio-temporally varying mixture fraction. Thus, one would expect turbulent mixing of combustible MILD mixtures, combustion, and their mutual influences to play a vital role in the establishment of MILD combustion. These effects resulting from spatio-temporal variations of mixture fraction on the structure of MILD reaction zones are yet to be explored, which provide a motivation for this work.
The specific objectives for this study are (i) to conduct DNS of turbulent combustion of MILD mixtures with spatio-temporally varying mixture fraction using complex chemistry, including the kinetics for a chemiluminescent species like OH * ; (ii) to gather physical insights through detailed interrogation of the numerical data and (iii) thereby shed light on the influence of mixture fraction variation on the structure of turbulent MILD reaction zones and the physical mechanisms involved.
This paper is organised in the following way. The methodology to conduct the DNS is described in Section 2 along with the chemical kinetic mechanism used in this work. Results are presented and discussed in Section 3 . Conclusions are summarised in the final section.

Conservation equations
The conservation equations solved, using SENGA2 [19] , are for mass, momentum, total internal energy, E , and species mass fractions, Y α , for (N − 1) chemical species in fully compressible form.
These equations are written using standard notations as mass : ∂ρ ∂t momentum : The mass fraction of the N th species is obtained using The various transport properties of the fluid mixture are temperature dependent and Fickian diffusion is used for molecular diffusive flux of species. The viscous stress tensor, τ ij , and the heat flux vector, q j , are defined as where μ and λ are, respectively, the mixture dynamic viscosity and thermal conductivity, h α is the enthalpy of species α. The temperature dependence of the transport coefficients are approximated using fifth order polynomials. The diffusion coefficient of species α is estimated using non-unity constant Lewis number for the species involved in the chemical mechanism. The thermal equation of state for the mixture is where R is the universal gas constant and W α is the molar mass of species α.
The above equations are discretised on a uniform grid detailed in Section 2.5 and the spatial derivatives are obtained using a tenth order central difference scheme which reduces to fourth order near the boundaries. The discretised equations are advanced in time using a low-storage third-order Runge-Kutta scheme although a fourth-order scheme can be used in SENGA2. This code has been used in many previous studies of MILD [11,12,15,16] , premixed [20][21][22] and droplet [23][24][25] combustion.
The above equations are to be specified with initial and boundary conditions before numerical solution can be started. These conditions depend on the flow configuration, which is described next, followed by a description of initial and boundary conditions used here. The choice of chemical mechanism and combustion conditions are described towards the end of this section.

Numerical flow configuration
MILD combustion of mixtures having spatio-temporal variation of mixture fraction with recirculated exhaust gases is of interest in the present work. However, typical MILD combustion systems are of a few tens of centimeters in size involving EGR or either spatially or temporally staged fuel injection techniques. DNS of such system is not feasible at this time because of the heavy computational costs involved and thus a simplified two-stage approach is used here following earlier studies [12] . The first stage mimics the mixing of fuel, air and exhaust gases yielding inhomogeneous mixture field which subsequently undergoes combustion in the second stage. This two-stage process is illustrated schematically in Fig. 1 . The left box titled "Preprocessing" represents the first stage elaborated in Section 2.3 . The symbol ˙ Q out represents the cooling of exhaust gases.
The computational domain for MILD combustion (second stage) is a cuboid with periodic boundary conditions in the y and z directions. Reflecting inflow and non-reflecting outflow are specified for the x -direction using the Navier-Stokes characteristics boundary condition (NSCBC) approach [26] . Spatially and temporally varying mixture of vitiated air and fuel (generated in the first stage) is fed at an average velocity of U in through the inlet located at x = 0 of the computational domain. The turbulent velocity, u , species mass fractions, Y i , of the partially premixed mixture and temperature, T , of the inflowing mixture are specified as where u , Y i and T are the turbulent velocity, species mass fractions and temperature obtained during the (first) preprocessing stage to be explained in Section 2.3 . The x -location of the scanning plane, x ( t ), moves with a velocity of U in = 20 m / s inside the preprocessed field.

Generation of initial and inflow fields
The flow and mass fraction fields for the combustion DNS are generated by following the methodology described in [11,12,15] . This preprocessing step mimics the mixing of vitiated air with fuel (first stage of the two-stage approach) before auto-ignition occurs. Earlier studies [11,12,15] considered premixed MILD combustion and thus there was no variation of mixture fraction. Since the influence of mixture fraction variation on the MILD reaction zone structure is of interest for this study, appropriate modifications are introduced to generate fields which are representative of non-premixed MILD combustion system. The steps described below are followed to obtain the desired fields of u and Y i .
Step 1: A turbulence field is generated first by conducting a DNS of freely decaying, homogeneous isotropic turbulence inside a periodic cube. The initial turbulence field is specified using a prescribed turbulent kinetic energy spectrum as described in [27] . This simulation is run until the turbulence is fully developed, which is recognised when the velocity derivative skewness reaches a value of about −0 . 5 .
Step 2: A one-dimensional freely propagating premixed laminar flame for a wide range of mixture fraction (or equivalence ratio) are simulated using the same kinetic mechanism to be employed for the combustion DNS. The reactants considered are a mixture of methane and diluted-air (without intermediate species) with desired dilution levels (specified in Section 2.5 ). Based on the earlier studies [11,12,15] , the reactant temperature is set at T r = 1500 K , which is higher than the reference auto-ignition temperature. This step yields a database of one-dimensional laminar premixed flames having thermochemical conditions of interest and these flames are parametrised using the mixture fraction (or equivalence ratio) and progress variable.
Step 3: Following the method described in [28] , two scalar fields bounded between 0 and 1 are generated by specifying scalar-energy spectrum functions having prescribed length scales, c and Z , and means, c and Z . These two fields are taken to be the initial reaction progress variable, c Y , and mixture fraction, Z , fields. The progress variable is defined as where the subscripts r and p in the methane mass fraction, Y CH 4 , denote, respectively, the reactant and product mixtures of a laminar flame computed in Step 2. The mixture fraction is defined using Bilger's definition [29] : (12) where the subscripts f and ox denote fuel and vitiated air streams respectively (see Fig. 1 ), and β = 2 with the elemental mass fractions and atomic masses for the elements carbon, hydrogen and oxygen. Typical variations of these two fields generated thus are shown in Fig. 2 for the mid x − y plane of the computational domain. As expected c Y presents variations bounded between 0 and 1 by construction with a prescribed lengthscale. On the other hand, Z has a smaller range of variation with a typically larger lengthscale. The probability density functions (pdfs) of these two fields are also shown in this figure. The pdfs of c Y for the various cases show a clearly bimodal behaviour for the choices of c and c . On the other hand, the pdfs of Z vary according to the case considered. For cases AZ1 and BZ1, the pdfs have a peak value centred around the specified mean value of Z while for case AZ2, the pdf does not present a peak value. It should be noted that there are no pockets of pure fuel in the present domain as a result of the small value for Z and the limited variations for Z . Further details about various cases considered are provided at the end of this section and in Section 2.5 .
Step 4: The species mass fractions obtained in Step 2 are used to initialise the scalar mass fraction fields depending on the local values of Z and c Y . If the local value of Z is beyond the flammability limit then the mixing line solution in the mixture fraction space is used to obtain the local species mass fraction values. The temperature is set to 1500K . The scalar field obtained thus has partially premixed mixture (varying equivalence ratio) containing reactants, partially burnt mixtures and products. These conditions would be quite similar to those in the recirculation zones of a MILD combustion burner investigated experimentaly, for example in [5] . The fluctuations in the scalar mass fraction fields obtained in this step do not have any correlation with the turbulence field obtained in Step 1. This lack of correlation may lead to some non-physical transients, which are overcome using the method described in the next step.
Step 5: The scalars and the turbulence from Step 1 are allowed to interact and evolve inside the periodic domain without any chemical reactions for a duration of about one large eddy timescale, 0 / u where 0 and u are the integral length scale and root-mean-square (rms) value of the turbulence field. Based on past studies, the duration of this simulation was kept deliberately below the ignition delay time to avoid occurrence of combustion but sufficiently long for the development of correlation between scalar and velocity fields through their interactions. The scalar fields obtained at the end of this step are shown in Fig. 3 and it is clear that unburnt ( c Y = 0 ), partially burnt (intermediate values of c Y ), and fully burnt ( c Y = 1 ) mixtures are present. Furthermore, the equivalence ratio, ranges from 0 to 10 inside the computational domain. Intermediate species and radicals are also present as seen in Fig. 3 d. These fields serve as the initial and inflowing fields for the combustion DNS.
A desired level of dilution is obtained by changing the reactant mixture used for the laminar calculations in Step 2. The length  scales c and Z denote typical eddy sizes for the progress variable and mixture fraction fields respectively and they will influence the spatio-temporal evolution of mixing and reactions. The chemical (flame) length scales are usually smaller than the mixing length scale and so the case of c > Z is excluded. Hence, three cases, which are discussed in Section 2.5 , are considered here. The first two cases, noted AZ1 and AZ2, have almost the same dilution level but different c / Z values. The third case, called BZ1, has the same length scale ratio as the first case but substantially higher dilution level. The PDFs of c Y and Z for the initial and inflowing fields of these three cases are shown in Fig. 4 . The maximum Z observed in the domain is about 0.1 which corresponds to φ ≈ 11 , but the pdf for Z > 0 . 04 is small and thus is not shown in Fig. 4 . The value of Z = 0 . 04 corresponds to φ = 4 . 1 for case AZ1 and AZ2, and φ = 7 . 1 for case BZ1, implying a substantial variation in equivalence ratio. Compared to the initial pdfs of Z shown in Fig. 2 , the pdfs of Z after the mixing in step 5 have not changed much since this mixing simulation time is kept to be about one eddy turnover time and the peak of the Z pdf is already close to Z . On the other hand, the pdfs of c Y presents much more variation showing the presence of unburnt, intermediate and burnt mixtures. The mixing process has also created high probability for intermediate values of c .
It should be noted that for the Step 2 described above, the choice of using the laminar premixed flames to construct the initial scalar field was made for the following reasons. The variations of Y i with Z and c are observed (not shown here) to be very similar for premixed, PSR and counterflow non-premixed configurations under fully reacting conditions. Such a variation is observed for most of the species except for some minor radical species. Furthermore, Z and c vary in the local normal direction across the reaction zone and thus the reaction zone structure in Z − c space do not depend on the flow configuration. The reaction zone structure depends on the scalar dissipation rate of the mixture fraction acting on the stoichiometric point in the non-premixed condition and thus the dissipation rate is a natural parameter, rather than c , to characterise the structure. It is also possible to map the non-premixed reaction zone structure and mass fractions fields using Z and c using a number of non-premixed flames generated by varying the mixture fraction dissipation rate at the stoichiometry [30] . However, the variations of Y i ( Z, c ) are very similar for these three configurations as noted earlier and thus the DNS solution and the results reported here are not influenced by the choice of the laminar premixed flame configuration.

Combustion mechanism
The fuel considered for this study is methane as in previous studies [11,12,15] with similar turbulence and dilution conditions. This helps us to contrast the difference between premixed and non-premixed turbulent MILD combustion. The skeletal mechanism of Smooke and Giovangigli [31] is selected following earlier studies. Since the chemiluminescent species such as OH * and CH * are also of interest from the perspective of experimental investigation [32,33] , the Smooke and Giovangigli mechanism is modified to include precursor species required for OH * chemistry taken from Kathrotia et al. [34] . The elementary reactions required for the precursor species are taken from Bilger et al. [29] . The kinetic mechanism assembled thus involves 19 species and 58 reactions in total and this mechanism is listed in Appendix A along with sources for the respective rate constants. This modified mechanism is referred as MS-58 and that from Bilger et al. [29] is called KEE-58 in the discussion below. Figures 5 and 6 compare the measured and computed laminar flame speeds and ignition delay times respectively. The laminar flame speeds are computed using Cantera [37] with 4 different chemical mechanism, viz., GRI-3, Smooke and Giovangigli, MS-58 and KEE-58. The experimental results are taken from [35] and [36] . It is clear that the MS-58 mechanism assembled using reactions from [31] and [29] predicts the laminar flame speed, s L , variation with equivalence ratio, φ, for methane-air mixture very well.
These results are comparable to GRI-3 values and MS-58 slightly under-predicts s L for φ ≥ 1.4.
The computed ignition delay times for a stoichiometric methane-air mixture at 1 and 4 atm over a wide range of temper- ature is compared to measurements in Fig. 6 . As one would expect, the GRI-3 estimates the ignition delay time accurately and, Smooke and Giovangigli and MS-58 mechanisms overestimate the ignition delay time for both pressures but the level of overestimation, in comparison to GRI-3, remains almost constant over the temperature range. The KEE-58 mechanism substantially overestimates the ignition delay time as shown in Fig. 6 for the case at 4 atm. At 1 atm, it performs better than the other two mechanisms for relatively high temperature, but substantially over predicts the ignition delay time for the temperature of interest here, T = 1500 K. It seems that the mechanism MS-58 is very good for flame propagation, quite good for ignition delay time estimates and also includes OH * precursors, and thus it is adequate to address the objectives of this study. However, the over estimate of the ignition delay time must be borne in mind while investigating the simulation results. As expected, the reactions involving OH * do not affect these estimates while allowing us to track OH * in the DNS. It is also worth remarking here that CH * and other chemiluminescent species are not included because C 2 and higher species acting as precursors for them are not available in MS-58 mechanism. The accuracy for the ignition delay times can be improved by including C 2 chemistry. However, including a comprehensive mechanism such as GRI-3 with C 2 and higher species is prohibitively expensive for DNS at this time. Also, there is no chemical mechanism with OH * precursor, which is suitable for DNS and so MS-58 mechanism is assembled for this study.

MILD combustion conditions
Two dilution levels are considered, the reference case has 3.5% oxygen by volume in the oxidiser mixture, which is the same as for the Case B of Minamoto et al. [12] in terms of dilution conditions. The earlier study considered the premixed case but this study is with non-premixed counterpart which includes spatio-temporally varying mixture fraction. The oxygen level is reduced to 2% by volume for the case BZ1. The mole fractions of major species in the diluted oxidiser mixture for these two conditions are listed in Table 1 . The sum of the mole fractions must be equal to 1 and so the rest of the mixture is N 2 .
The turbulence and thermochemical conditions for the combustion DNS are summarised in Table 2  There are about 7 integral eddies inside the computational domain of size L x = L y = L z = 10 mm . Since there are spatial variations of mixture fraction and progress variable (unburnt and burnt mixtures) fields, the ratio of their integral length scales is an important parameter. The influence of this parameter is also investigated by changing the value of c / Z from 0.77 to about 1 while keeping the dilution level the same. This case is denoted as AZ2 in Table 2 . The condition of c / Z > 1 is seen to be unphysical because the thickness of either a propagating flame or an ignition kernel at a reactant temperature of about 1500 K is expected to be thinner than a  [38] and [39] . Table 2 DNS initial conditions for MILD combustion. representative mixture fraction length scale resulting from turbulent mixing of fuel and vitiated air. From Table 2 , it can be observed that all cases correspond, on average, to lean condition with an equivalence ratio of approximately 0.8. The variance of the mixture fraction in case AZ2 is higher than for case AZ1 as that case was initialised with a smaller Z which allows for more variation of Z . Furthermore, the maximum O 2 concentration found in the computational domain for all cases correspond to the concentration from the oxidiser composition (see Table 1 ) meaning that pure oxidiser (vitiated air) is present in the domain. All cases share similar statistics, c and σ 2 c , for the progress variable field which is to be expected because similar initial c field was used. An average value of c equal to 0.56 signifies the presence of both burnt and unburnt mixture.
The computational domain is discretised using 512 × 512 × 512 uniform grid points which ensures that there are about 30 grid points inside the smallest chemical thickness for the mechanism MS-58 used here. This thickness is related to OH * species reaction rate. The timestep used is 1 ns which is smaller than the shortest chemical timescale imposed by OH * chemistry.
Each of these simulations has been run for about 1.5 flowthrough time, which is defined as τ f = L x /U in . Statistics were collected after the first flow-through time to ensure that the transients from the initial conditions had left the domain before data sampling. These simulations have been run on ARCHER, a Cray XC30 systems, and each simulation took about 550 hours on a wall-clock using 4096 cores.

Results and discussion
The overall influence of combustion can be investigated using the variations of averaged, in the y − z plane and in time, quantities such as temperature and various scalar mass fraction values with x . Analysis of these quantities shows that the averaged temperature increases progressively from 1500 K to about 1600 K and the reactant species such as fuel and oxygen mass fractions decrease with x while the product mass fractions increase gradually. Intermediate species like OH or CH 2 O show a small decrease from their inlet values because of their consumption to initiate combustion reactions and then increases as they are produced by combustion. These are expected behaviours and so they are not shown here.

Instantaneous reaction rate
The typical iso-surface of normalised heat release rate, ˙ Q + = ˙ Q δ th / (ρ r s L C p (T b − T r )) is shown in Fig. 7 for ˙ Q + = 2 . 0 from the case AZ1. The normalisation is done using the thermo-chemical quantities, appearing in the above expression, for the local mixture fraction. One can also use these quantities for Z st in the normalisation but this would mask regions with relatively low heat release rate per unit volume, for example lean regions. It can be observed that reactions occur over a very large portion of the computational domain with an extremely convoluted aspect. This increases the possibility for interactions of reaction zones and clearly differentiates MILD combustion from conventional combustion where a clear flame front with localised heat release rate is usually present. Furthermore, reactions occur near the entrance of the computational domain, shown by the presence of the iso-surfaces there (see Fig. 8 ). This is due to the elevated temperature and the presence of radicals that initiate reactions. The variation of this isosurface in the other two cases and also at other times is similar to that shown in Fig. 7 . Movies of these iso-surfaces and of the mid x − y cuts are available for all cases as supplementary materials.
For a close analysis of these reaction zones, 2D slices in the mid x − y plane are shown in Fig. 8 for the three cases. One observes both thin and thick heat releasing zones in Fig. 8 a and b. The thickening of the reaction zones originates from interactions of reaction zones yielding larger regions of chemical activity and heat release. These regions can then occupy large patches in the computational domain. The typical thickness of these reaction zones present a large variation as pictured in Fig. 8 a and  of about 2% shows chemical activity throughout the domain. In addition, the degree of convolution of the reaction zones between cases AZ1 and AZ2 appears to be similar despite the different initial mixture fraction lengthscales (see Table 2 ). Comparing the cases AZ1 and BZ1, shown respectively in Fig. 8 a and c, the case BZ1 has much thicker reaction zones nearly over the entire computational domain and no region could be identified as a thin reaction zone as seen for AZ1 and AZ2. This finding is in agreement with those observed in premixed MILD combustion investigations [12] suggesting that a higher dilution level led to thicker reaction zones with more frequent interactions of reaction zones.
The evolution of species in these reaction zones is studied in Fig. 9 to investigate the physical space structure of these zones. The spatial evolutions of ˙ and c T at the locations marked in Fig. 8 a are shown. The * superscript denotes quantities normalised by the maximum value in the plane. The distance d in Fig. 9 is along the local normal n = −∇ c T / |∇ c T | and d = 0 corresponds to the point along n with peak heat release rate per unit volume. The variations of the above quantities with d are shown for a non-interacting reaction zone (marked using filled circles in Fig. 8 a) and an interacting one (open circles) in Fig. 9 b and d, respectively. Figure 9 a and c shows these variations for the same locations at an earlier time, t = 1 . 46 τ f . For the non-interacting zone ( Fig. 9 a and b), a clear localised heat release profile with a progressive increase of c T is observed and this variation is typical of a flame. The variations of normalised methane and oxygen mass fractions shown in the figure suggest that this reaction zone is of premixed type. However, there is a substantial variation of mixture fraction across this reaction zone (0.49 ≤ φ ≤ 1.15) and this variation is well within the flammability limits for the conditions analysed here. The propagation of this reaction zone is tracked through detailed interrogation of the time evolution of ˙ Q + shown in Fig. 8 . This analysis and the results in Fig. 9 show that there is a propagation of this flame in the downstream direction (compare Fig. 9 a  and b).  Figure 9 c is shown for t = 1 . 46 τ f depicting two distinct peaks for the heat release rate. The variations of other quantities shown clearly suggest the presence of two distinct flames, which start to interact as they are propagating towards each other. This interaction leads to a broader peak for the heat release rate as seen in Fig. 9 d shown for t = 1 . 5 τ f . This is confirmed further by the variation of OH along n showing a clear discrepancy with the variation of heat release as observed in earlier studies on interactions of laminar flames [40] . Conversely, OH * seems to adequately follow the variation of heat release. This suggests that OH may not be a suitable surrogate to identify reaction zones in nonpremixed MILD combustion as also observed for premixed MILD combustion [11] .

Reaction zone structure in mixture fraction space
As mentioned in Section 1 , one of the main objectives of this study is to analyse the effects of spatial and temporal variation of the mixture fraction on MILD combustion. The variations of typical quantities such as temperature, heat release rate, major and minor species with mixture fraction are shown in Figs. 10-14 for case AZ1. The results are shown for two streamwise ( x ) locations and two instances which are marked in the respective figures. The conditional averages for the quantities are shown using lines with symbols with an error bar of width representing their 2 standard deviation σ from the conditional mean value. The dash-dotted lines represent the variation in a laminar counterflow flame with a very low ( 50 s −1 ) strain rate computed using Cantera [37] . This laminar flame has pure fuel as one stream and the vitiated air (see Table 1 ) used for case AZ1 as the other stream. This strain rate has been chosen to ensure a fully reacting solution. The choice of this value is, however, not critical for MILD combustion as there are no The mixing line for the oxidiser at 1500 K and the fuel at 300 K is also shown as dashed line in Fig. 10 which shows the variation of T with Z . The time variation of the incoming mixture fraction is observed by comparing Fig. 10 a and c. The maximum Z seen at t = 1 . 5 τ f is nearly 30% of the value at t = τ f . For the location close to the inlet, most mixture is at about 1500 K represented by the conditionally averaged temperature irrespective of the time of analysis, see Fig. 10 a and c. However, there is a relatively larger scatter for the lean mixture because of its lower ignition delay compared to rich ones [14] allowing some of these lean mixtures to have already reacted. As one moves downstream ( Fig. 10 b and  d), the conditionally averaged temperature increases for all mixtures due to heat release. The peak temperature in the turbulent case is observed to occur for lean mixture and not for a slightly rich mixture as in the counterflow flame. This is due to interactions (mixing) between the various mixture fractions present in the flow. The stoichiometric and slightly rich mixtures which are hotter diffuse heat towards the leaner (slightly cooler) mixture. As a result, the conditional averaged temperature of these hotter mixture decreases while the one of the leaner mixture increases. This yields a peak in temperature located around the average mixture fraction inflowing in the domain (of Z ≈ 0.008 for case AZ1). Furthermore, there are ignition fronts as one shall see shortly.
The peak of conditional and also unconditional temperature in the turbulent case is smaller than the temperature observed for the counterflow flame because the inflowing mixture at 1500 K in the DNS contains pockets of burnt and unburnt mixtures with the same equivalence ratio (see Fig. 3 ) which is not the case for the steady counterflow flame shown in Fig. 10 . The spatial variation of the conditional temperature at t = τ f can be understood by comparing Fig. 10 a and b. The time evolution at a given location can be observed by comparing the results shown in the top and bottom rows of Fig. 10 . The spatiotemporal variation of the mixture fraction can also be observed in these figures. The range of Z decreased nearly by half from about 0.06 to 0.03 when comparing Fig. 10 a to d. This change is because of turbulent mixing of the inhomogeneous incoming Z field and when this mixing is complete the mixture fraction will be equal to the mean value of 0.008 for the case AZ1 (see Table 2 ).    Fig. 10 a and d. Figure 10 b and d also shows the evolution of conditionally averaged T for case BZ1. Similar observations can also be made for this case.
The temperature evolution shown in Fig. 10 clearly suggests that there are ignition fronts moving in space and time. This front evolution in space can be seen clearly in Fig. 11 showing the evolution of conditional temperature in the mixture fraction space for five streamwise locations. One can clearly see the propagation of ignition front in space. If one plots this variation at a given location for various instances then one observes a similar behaviours (see Fig. 10 b and d for example). The mixtures ignited at an earlier location and time yield propagating flames subsequently. Figure 12 shows the variations of ˙ Q with Z . It is observed that the maximum heat release rate occurs for a slightly lean mixture for both the counterflow flame and turbulent case, and not for a slightly rich mixture as observed for a conventional flame. This behaviour is seen for all locations and times in the turbulent case. Furthermore, the conditional average of ˙ Q in the turbulent case is higher than the laminar counterflow flame. This is because the turbulent case contains mixture of different thermo-chemical state for a given mixture fraction value, leading to a large fluctuation in ˙ Q as shown by the vertical bars (representing the conditional rms values). It is also to be noted that the incoming mixture contains intermediates and radicals, which is not so for the counterflow case.
Further insights into these structures can be gathered by studying the mixture fraction space variations of major and minor species. This variation for major species, CH 4 and O 2 , is shown in Fig. 13 . The mass fraction of CH 4 follows closely the laminar case for the lean mixture, ie., Z < Z st = 0 . 01 for the locations and instances shown in the figure. For other mixture fraction values, the mass fractions are substantially larger compared to the fully reacting laminar case. This behaviour is typical of ignition phenomena. For example, if one compares the temperature, ˙ Q and these mass fraction values from Figs. 10 a, 12 a and 13 a, respectively, for Z = 0 . 03 then one can see that this mixture has not ignited yet. The results in these figures show that this specific mixture has not reacted even at x / L x ≈ 0.67. This is because this mixture at this specific location has emerged from a richer mixture through turbulent mixing which would drive the mean Z value towards Z = 0 . 008 .   This physics is reflected by the abrupt ending of the various curves in Fig. 11 at some intermediate values of Z . The variations of minor species for the same locations and times are shown in Fig. 14 . It is clear that these species variations do not follow those in counterflow configuration. For x / L x ≈ 0.03, which is close to the inlet, the reactions are just beginning to occur through ignition by consuming the incoming OH (see Fig. 3 ) and thus OH levels are substantially lower compared to the laminar case. The formaldehyde level is higher because it is also present in the incoming mixture. This gives substantially larger level even for unreacted mixture fraction value of Z ≥ 0.004. As the ignition front moves through the mixture fraction space, OH and CH 2 O mass fractions move towards counterflow flame values. The substantially larger values seen for Z ≥ 0.01 is because the ignition front has not moved to these mixtures yet. Thus, these variations are consistent with ignition front behaviour seen in earlier figures. For cases AZ2 and BZ1, the mixture fraction space variations of species and their spatio-temporal evolution are observed to be similar to that shown here for AZ1.
For example see the variations of conditionally averaged temperature with Z shown in Fig. 10 for the case BZ1. However, it is worth to note that the effects of c / Z and the dilution will be felt in the physical space but not in the mixture fraction space. The physical space structure shown in Figs. 8 and 9 , and the mixture fraction space evolution shown in Figs. 10 -14 suggest that nonpremixed MILD combustion involves ignition fronts, propagating premixed and non-premixed reaction zones. The later two aspects are analysed further in the next section. The presence of ignition fronts and propagating premixed flames was also observed for premixed MILD combustion [16] . However, their relative roles, contributions, and mutual influences are subject matters for a future study.

Premixed and non-premixed behaviour of reaction zones
Since the mixture temperature is higher than the ignition temperature for methane, autoignition events occur which can subsequently lead to either premixed or non-premixed combustion depending on the local mixture condition. The premixed characteristics arise because both oxidiser and fuel are present in the incoming partially premixed mixture. The spatio-temporal variation of mixture fraction yields non-premixed combustion features. This dual behaviour can be investigated using the flame index introduced by Yamashita et al. [41] which was modified by Briones et al. [42] . The modified flame index, FI, is given by: This FI provides information about the gradient (molecular mixing) direction and thus the combustion mode when computed in heat releasing regions. It takes a value close to 0 for non-premixed combustion, a value close to −1 for lean and +1 for rich premixed combustion. The spatial variation of FI in regions with substantial heat release rate per unit volume ( ˙ Q + ≥ 1 ) is shown in Fig. 15 for all three cases. This result is shown in the mid x − y plane of the computational domain.
This figure shows islands of non-premixed combustion surrounded by regions of premixed combustion (few regions showing these characteristics are marked in the insets of Fig. 15 ). The regions of non-premixed combustion originates from the flux of fuel coming from a neighbouring richer mixture and the flux of oxygen coming from an adjacent leaner mixture. As a result of this, the gradients of fuel and oxygen mass fractions have opposite direction giving FI = 0 (see Eq. (13) ). This situation is illustrated in the insets of Fig. 15 a-c where there is a richer mixture on one side and a leaner mixture on the other side. For the inset of Fig. 15 a,  Fig. 15. Flame Index shown in regions with ˙ Q + ≥ 1 . 0 for cases (a) AZ1, (b) AZ2, (c) BZ1 at t = 1 . 5 τ f . The results are shown in the mid x − y plane. despite a clear localised peak of heat release rate, similar to those found in premixed flame, a further analysis to be presented later in Fig. 17 shows that this heat release still comes from non-premixed mode of combustion. The inset in Fig. 15 b shows a non-premixed island which has just started to react as indicated by the relatively low ˙ Q * and low c T (the superscript * implies that the value is normalised using the respective maximum value seen in the corresponding x − y plane shown). As it has just begin to ignite, the levels of OH and OH * are low. Furthermore, by comparing the case AZ1 to AZ2, it is observed that the proportion of non-premixed combustion is higher in the case AZ2. This is linked to the relatively larger value of c / Z for  the case AZ2 compared to the case AZ1 (see Table 2 ). This creates steeper gradients of mixture fraction in case AZ2 leading to stronger and more frequent fluxes of CH 4 and O 2 coming from opposite directions. The flux of CH 4 comes from richer and that of O 2 from leaner mixtures.
The results depicted in Fig. 15 show the complex nature of MILD combustion and the coexistence of premixed and nonpremixed modes. It is to be noted that FI does not help to identify autoignition. The preponderance of these modes and their relative contribution can be investigated by studying the pdf of FI, p ( ψ).
This pdf constructed using 40 bins is shown in Fig. 16 for the three cases. The result is shown for t = 1 . 5 τ f and it is similar for other times. Note that the y -axis is shown in logarithmic scale to clearly see the pdf variation for intermediate values of FI. Let us remind ourselves here that FI = −1 corresponds to lean premixed while +1 corresponds to rich premixed and 0 implies non-premixed combustion modes. It is quite clear that the lean premixed mode is dominant in all three cases. The rich premixed mode is the next dominant mode in the case AZ1, which is reduced in the cases AZ2 and BZ1. This drop in the case BZ1 is because of high dilution level. The stronger mixture fraction gradient, resulting from smaller Z , Table 3 Volume fractions of regions having non-premixed, rich premixed, lean premixed and mixed-mode combustion at t = 1 . 5 τ f . Their fractional contributions to the total heat release rate in the domain are also listed. in AZ2 compared to AZ1 decreases the rich premixed mode. However, the non-premixed mode increases in AZ2 and BZ1 compared to AZ1. Even the mixed mode, having 0.1 ≤ |FI| ≤ 0.8, increases in these two cases. These changes can be quantified by calculating the fraction of the computational volume with these combustion modes.
These fractions can be obtained by simply using: for non-premixed, lean premixed and rich premixed combustion modes respectively while the fraction of the mixed mode is com- The fractions obtained thus are summarised in Table 3 . The symbol ψ is the sampling space variable for FI and, ξ 1 = 0 . 1 and ξ 2 = 0 . 8 . These limits are chosen carefully after extensive testing on the sensitivity of the respective integrals to these parameters. It is clear from those values that premixed combustion is predominant with the lean condition being the main contributor in the present configuration. This is to be expected since all cases have lean mixture with φ = 0 . 8 on average. Furthermore, in the configuration considered here, fuel and oxidiser are introduced as a partially premixed mixture which naturally yields a dominant premixed combustion mode. However, the volume fraction having non-premixed combustion is non-negligible, which increases substantially for the most diluted case BZ1. As noted earlier, the increased mixture fraction gradient resulting from larger c / Z in case AZ2 yields larger V NP . Nearly 70-80% of the domain experiences these classical modes of combustion and the remaining 20-30% has mixed mode combustion, ie., with 0.1 ≤ |FI| ≤ 0.8. The pdf of FI constructed using planar data for various streamwise ( x ) locations has also been studied for all cases. These pdfs are similar to that shown in Fig. 16 and it is observed that the three combustion modes exist at all locations but with varying level of their relative importance. In the early stage of the domain, the lean premixed mode is dominant with nearly nonexistent rich premixed and non-premixed modes. This originates from the early ignition of lean mixtures. Further downstream, non-premixed mode starts to appear with the flux of O 2 coming from the excess oxygen in the burnt lean mixture and CH 4 coming from unreacted or partially reacted rich mixtures. Rich premixed mode also appears only in the downstream region as the ignition delay time for rich mixtures is larger than for lean ones.
Following the analysis of non-premixed combustion by Bilger [43] , the heat release rate per unit volume in non-premixed mode can be written as: where h s is the sensible enthalpy of the mixture, through the use of mixture fraction theory. Similar non-premixed contribution was shown to exist by Bray et al. [44] using Z and c to investigate partially premixed combustion. This contribution can be estimated directly using the DNS data. Denoting the total heat release rate per unit volume at a given point as ˙ where h 0 f i is the enthalpy of formation of species i , one expects that ˙ Q NP / ˙ Q tot should be close to unity in regions with FI ≈ 0. This is simply another way to verify that the combustion in those regions is through non-premixed mode or not. The spatial variation of this ratio in the mid x − y plane for the case AZ1 is shown in Fig. 17 . The colour map shows this ratio and the contour lines of | FI | = 0 . 1 are used to mark regions with |FI| < 0.1. This inequality is used to include the contributions coming from the slightly longer negative tail in the pdf of FI shown in Fig. 16 . A good agreement between the regions having FI ≈ 0 and ˙ Q NP / ˙ Q tot ≈ 1 is observed confirming the nonpremixed nature of combustion in these regions (see Fig. 15 for FI).
The fractional contribution, Q , of these various modes of combustion to the total heat release rate for the entire domain is given in Table 3 . These fractions are obtained by computing the ratio of the heat release rate per unit volume, ˙ Q tot , integrated over these regions to the heat release rate integrated over the entire computational domain. It is observed that nearly 50% comes from lean premixed combustion and the contribution from rich premixed mode drops below 3% as dilution or c / Z is increased. However, the nonpremixed contribution still amounts to a non-negligible portion of the total heat release rate. The contribution from mixed mode is about 25%. Can this be seen to come from ignition fronts or not?
This is an open question at this time and will be answered in future work.

Summary and conclusions
DNS of MILD combustion with spatio-temporal variations of mixture fraction, Z , is performed to investigate the MILD reaction zone structures. The conservation equations are solved using SENGA2, which is based on fully compressible formulation. The computational domain is a cube with inflow and non-reflecting outflow boundary conditions in the streamwise ( x ) direction and periodic conditions in the other two directions. The required initial conditions are constructed carefully in a preprocessing step as described in Section 2.3 , which allowed us to include unburnt, partially burnt and fully burnt mixtures for a given mixture fraction value. Thus, the inflowing mixture has spatially and also temporally varying Z and progress variable fields. The preprocessing is used to mimic the mixing of fuel with vitiated air coming from the recirculated exhaust gases in a MILD combustion burner. The local thermo-chemical state of the inflowing mixture is determined using mass fractions of various species computed for laminar flames with the range of mixture fraction observed in the initial Z field. The temperature of the inflowing mixture is set to be 1500 K.
The combustion kinetics is modelled using a skeletal mechanism for methane-air combustion and this mechanism is assembled using the reactions in [31] , [29] and [34] in order to include chemiluminescent species like OH * . The assembled mechanism is tested for its ability to predict the laminar flame speeds for the range of equivalence ratio and the ignition delay times. These testing showed that the assembled mechanism is good with some overprediction of ignition delay time.
A close scrutiny of the DNS results offered the following insights. The chemical reactions of MILD mixtures occur over a wider region of the computational domain unlike the conventional combustion typically showing chemical reactions confined to narrower regions. As noted in earlier studies [12,15,16] , the MILD reaction zones are observed to have their salient features, such as strongly convoluted zones with frequent interactions while retaining conventional flame characteristics. However, allowing the mixture fraction to vary introduces further complexities. Spatial and temporal ignition fronts are observed, leaner mixtures with shorter ignition delay time are observed to ignite first and these kernels evolved into propagating premixed flames. These flames remain as lean or become rich premixed flames depending on the mixture fraction value along their propagation paths produced by turbulent mixing of Z field. The rich mixtures are observed to ignite later at downstream locations of the computational domain. These features are observed clearly by analysing the spatio-temporal evolution of the reactions zones and their structures in the physical and mixture fraction spaces. More importantly, the presence of non-premixed combustion mode is also observed through an analysis of flame index, FI, defined in Eq. (13) . A value of −1 and 1 for FI implies lean and rich premixed combustion respectively and FI = 0 means non-premixed combustion as noted in [42] . The flame index does not distinguish ignition fronts. The presence of non-premixed mode was not observed in the previous MILD combustion DNS studies [12,15,16] as the mixture fraction variations were excluded. The non-premixed mode combustion results from the fluxes of methane and oxygen coming from neighbouring pockets of richer and leaner mixtures and the premixed and non-premixed zones are observed to coexist along with ignition fronts. The fractions of computational volume having lean and rich premixed, and non-premixed combustion modes are computed using the pdf of flame index and the corresponding fractional contributions of heat release rate are calculated. The results suggest that nearly 50% of the contribution to the total heat release rate comes from the lean premixed mode. The contribution from rich premixed mode drops from about 13% in the reference case AZ1 to about 3% in the diluted case BZ1 when the mean oxygen mole fraction in the inflowing mixture is reduced from 2.7% to 1.6% (see Table 2 ). Whereas the nonpremixed contribution increases substantially and this increase is also observed when the integral length scale of the mixture fraction field is reduced as noted in Table 3 . However, the sum of these three fractional contributions to the total heat release rate is noted to be about 75% only and the remaining 25% is found to come from regions with 0.1 ≤ |FI| ≤ 0.8. Whether this contribution can be attributed to ignition fronts or not is an open question. Finding a meaningful answer to this question requires further detailed investigation. To conclude, non-premixed turbulent MILD combustion is a complex process with spatial and temporal ignitions fronts, propagating flames that are strongly convoluted and interacting, non-premixed combustion and interaction among these. Future work will investigate the coexistence of these various phenomena to provide physical insights. Also, a skeletal mechanism with C 2 -chemistry is required to accommodate rich premixed and non-premixed combustion, and for improved prediction of autoignition delay times. Such a mechanism of reasonable size is of great interest for MILD combustion DNS investigations as long as the computational burden for chemistry is not too onerous.

Table A4
Methane-air combustion mechanism used with OH * chemistry. Rate coefficients are

Supplementary material
Supplementary material associated with this article can be found, in the online version, at 10.1016/j.combustflame.2017.10.030.