Evaluation of Modeling Approaches for MILD Combustion Systems With Internal Recirculation

Numerical simulations employing two different modeling approaches are performed and validated against experimental results from a moderate or intense low-oxygen dilution (MILD) system with internal recirculation. The flamelet-generated manifold (FGM) and partially stirred reactor (PaSR) closures are employed in a Reynolds-averaged Navier–Stokes (RANS) framework to carry out the numerical simulations. The results show that the FGM model strongly overpredicts temperature profiles in the reactive region, while yielding better results along the central thermocouple. The PaSR closures based on a prescribed mixing time constant, Cmix, of 0.01, 0.1, and 0.5 are compared, showing that a Cmix value of 0.5 is the most appropriate choice for the cases under investigation. A PaSR formulation allowing local estimation of the Cmix value is found to provide improved results for both the lateral and central thermocouples. A flame index analysis, used to assess the ability of FGM and PaSR to capture intense mixing of the cyclonic burner, indicates how the FGM model predicts a typical non-premixed region after the injection zone, contrary to the experimental observation.


INTRODUCTION
The simultaneous requirements of thermal efficiency, fuel flexibility, and low pollutant emission are the main drivers for the development of innovative combustion-based systems. Several concepts were developed in the last decades, to achieve lower temperatures of the oxidation region and avoid the formation of several classes of pollutants. In particular, a number of combustion concepts based on the internal or external dilution and/or preheating of main reactants were proposed. They include flameless oxidation (Wünning and Wünning, 1997), colorless combustion (Khalil and Gupta, 2017), high-temperature air combustion (HiTAC) (Rafidi and Blasiak, 2006), and moderate or intense low-oxygen dilution (MILD). Examples of implementation of these technologies include furnaces, boilers, and gas turbines.
MILD combustion plays today a crucial role in the development of flexible, efficient, and nonpolluting technologies. Notably, its feeding conditions allow the establishment of a distributed reaction zone, with almost uniform temperatures in the combustion chamber, absence of a visible flame, very low levels of noise and fluctuations, absence of soot particles, and negligible NO x and CO emissions (Cavaliere and de Joannon, 2004). In the last decades, strong attention was given to the turbulent combustion modeling of such processes. In fact, the unique thermochemical features related to the turbulence-chemistry interaction have led to the development of new modeling paradigms or to the adaptation of the ones already used for conventional conditions (De and Dongre, 2015). In such systems, the flue gases' recirculation decreases the oxygen level of the reactant mixture, thus slowing down the kinetic rates in such a way that the characteristic chemical timescales approach the mixing timescales; that is, the Damköhler number is close to unity (Özdemir and Peters, 2001). Due to such modifications of the reaction zone, appropriate turbulence/chemistry interaction models need to be developed (Frassoldati et al., 2010). In particular, detailed kinetic mechanisms should be considered to correctly capture the thermochemical variable distributions (temperature and main species fields) and pollutant emissions in particular.
The complexity of MILD reaction structures has pushed the development of a number of approaches for Reynolds-averaged Navier-Stokes (RANS) (Galletti et al., 2007) and large-eddy simulations (LES) (Duwig et al., 2007), which can be divided into flamelet-like and finite-rate approaches, depending on the assumptions made with respect to the accessible reaction state space (Pope, 2013).
Flamelet-based models allow us to include detailed chemistry effects at an affordable computational cost, as the chemical state space is calculated in a preprocessing step, and it is then accessed during the actual simulation via a number of scalars (typically the mixture fraction and a progress variable) that parameterize the chemical state space Moin, 2001, 2004). In this context, tabulated chemistry techniques such as the flameletgenerated manifold (FGM) (Van Oijen and De Goey, 2000) and flamelet progress variable (FPV) Moin, 2001, 2004) techniques have acquired a key role in the combustion community, and their effectiveness has been demonstrated for a number of problems and configurations. Recently, Abtahizadeh et al. (2015) showed that the canonical structure used for tabulation strongly influences the representation of MILD combustion in jet in hot coflow (JHC) configurations. Such results were confirmed by Sorrentino et al. (2018a), for a cyclonic flow burner. A relevant challenge in FGM/FPV approaches is to include the effects of exhaust gas recirculation in the lookup tables, to properly represent the dilution and heat loss effects (Lamouroux et al., 2014). The inclusion of additional variables in the tabulation procedure is a very demanding task and can lead to expensive computational fluid dynamics (CFD) simulations due to the increased manifold dimensionality (from 3D or 4D to 6D).
Reaction rate-based approaches have been also proposed in the literature and applied to MILD conditions. Among them, it is worth mentioning the eddy dissipation concept (EDC) (Magnussen, 2005) and the partially stirred reactor (PaSR) model, originally proposed by Chomiak (1990). Both models divide each grid cell into two regions: the reacting structures, modeled as ideal reactors, and the surrounding fluid mixture. The main difference between the two approaches lies in the estimation of the reacting structures' volume fraction. In EDC, the latter is solely based on the turbulence properties of the flow, obtained by means of an energy cascade model (Magnussen, 2005). In PaSR, the reacting fraction explicitly depends on both the mixing and chemical timescales. Recently, the PaSR model was used to simulate the Adelaide JHC, using both RANS (Ferrarotti et al., 2019) and LES (Ihme and See, 2011), with very satisfactory results. Despite the encouraging progress in the use of reactor-based models for MILD combustion, there are still open questions related to the general applicability of the concept to different configurations. To this end, the availability of experimental data in MILD conditions is key.
Both the Adelaide and Delft JHC configurations (Dally et al., 2002;Oldenhof et al., 2011) mimic MILD conditions by feeding diluted and hot streams to the coflow, impacting directly the system characteristic chemical timescale. This simplified configuration makes it possible to perform optical diagnostics and deliver high-fidelity experimental data that are crucial for model development. However, in industrial MILD configurations, the distributed conditions are a consequence of the internal aerodynamic recirculation and the modification of the local mixing timescales. Closed furnaces operating in MILD conditions allow, in principle, challenging of the ability of turbulent combustion models to appropriately capture turbulence/chemistry interactions. For instance, Ferrarotti et al. (2018) applied the PaSR model for a quasi-industrial burner fed with natural gas, showing the relevance of both mixing and chemical timescales as well as of the comprehensiveness of the chemical mechanism.
The availability of experimental data in furnace configurations is generally limited due to the presence of an enclosure. In addition, available data in MILD conditions are limited to narrow operative conditions. Based on such considerations, the present study reports the investigation of a novel cyclonic burner (Sorrentino et al., 2018b). A methane-fired small-scale burner was employed, and detailed internal temperature measurements were used to characterize the local thermal field. Experimental temperature measurements were used to validate two different modeling paradigms, a tabulated chemistry approach, FGM, and a reactor-based model, PaSR. The comparison was carried out for two experimental configurations, involving standard air and diluted air as the oxidizer stream. The main aim of the paper is to identify current potential and limitations of turbulent combustion closures for the investigation of more complex MILD combustion configurations.

EXPERIMENTAL TEST CASE
The geometrical features and characteristics of the cyclonic burner used in this article to mimic MILD combustion conditions have been reported in previous works (de Joannon et al., 2017;Ferrarotti et al., 2019;Sabia et al., 2019;Sorrentino et al., 2019). In the following, the key features of the device are discussed.
It consists of a lab-scale combustor (20 * 20 cm 2 of square section and a height of 5 cm) where a cyclonic and toroidal flow field is established through two couples of coaxial oxidant/fuel inlet jets located in an antisymmetric configuration. In Figure 1, a sketch of the burner midplane, the feeding configuration, and the position of the thermocouples are shown. Oxidizer and fuel injection pipes are placed 2 and 4.5 cm from the lateral walls, respectively. The flue gas exit is placed at the top-central side. In the present manuscript, methane was used as a fuel and air as oxidizer.
MILD conditions are ensured thanks to the strong internal recirculation of burned products that stabilizes locally a distributed ignition. Such a microscopic feature of MILD combustion corresponds macroscopically to reduced temperature gradients (high levels of temperature uniformity) with very low pollutant emissions.
Two shielded N-type thermocouples can be moved along the y-direction. The lateral one is placed 0.02 m from the wall, whereas the central thermocouple is located 0.01 m from the confinement.
External ceramic preheating ovens were used to fix the boundary conditions in terms of external temperature, thus controlling the heat flux through the surroundings. The burner was built with a three-layer structure of alumina, a thick layer of heat-insulating material, and an outer shell of stainless steel 310s. The oxidant flow passes through heaters, to keep fixed the inlet temperature to the desired values (T in ).
In the present work, the two experimental configurations illustrated in Table 1 were investigated. In the first case, the fuel (pure methane) was injected at the environmental temperature (300 K) with a jet velocity of 15.8 m/s, while the oxidizer (preheated air) was fed at a fixed preheating level of 1,000 K (T in ) and with a jet velocity of 17.8 m/s. The inlet equivalence ratio was kept constant at the stoichiometric value. Moreover, a mean overall heat loss flux at walls (Q W ) was estimated by means of thermocouples placed on outer and inner surfaces of the walls to be equal to 13 kW/m 2 , for the investigated conditions. In the second case, the air stream was diluted with 90% N 2 by volume. The equivalence ratio was kept to stoichiometric value, and the inlet temperatures were as in case 1; the fuel jet velocity was 15.8 m/s, while the oxidizer jet was 54.17 m/s.

MODELING TOOLS
Numerical simulations were carried out using RANS, with the commercial code Ansys Fluent 19.3. The grid was generated with the software Ansys Workbench. A grid independence study was performed using a polyhedral grid with a number of cells ranging from 100 to 400 k. The selected grid consisted of 450 k polyhedral

Case
Fuel elements. Reynolds stresses were solved using the RNG kε turbulence model with swirl-dominated flow corrections to take into account the high swirl level in the combustor. The discrete ordinate (DO) radiation model was used, while the radiation properties of the reacting mixture are accounted for with the weighted-sum-of-gray-gases (WSGG) model, using the coefficients proposed by Smith et al. (1982).
For the turbulence-chemistry interaction, PaSR and FGM models were adopted. As with other reactor-based models, the computational cell in PaSR is split into two zones, the reaction structures and the surrounding fluid. The overall reaction progress is determined by the exchange between the two regions. The mass fraction of the reaction zone is estimated, considering both the characteristic chemical and mixing timescales: A first set of simulations was performed to determine the influence of the kinetic mechanism on the results. The detailed mechanisms GRI2.11 (Bilger et al., 1990;Bowman et al., 1996) and GRI3.0 (Smith et al., 2011), along with the reduced mechanism KEE (Bilger et al., 1990), were used. Then, a sensitivity to the PaSR model parameters was carried out, with special emphasis on the determination of the mixing scale. First, τ mix was estimated as in Ferrarotti et al. (2018), considering a fraction of the integral timescale κ/ǫ by means of a constant C mix : Values of C mix ranging 0.01, 0.1, and 0.5 were considered. Then, a dynamic approach was used, as proposed in Ferrarotti et al. (2019), based on local properties of the flow field. τ mix,d is estimated with the ratio of the mixture fraction variance, Z ′′ 2 , to the mixture fraction dissipation rate, χ: The determination of τ mix,d requires the solution of transport equations for Z ′′ 2 andχ , expressed as where Z is the mixture fraction, D t is the turbulent diffusivity, ) is the production of scalar fluctuation, ) is the production of turbulent kinetic energy. The set of coefficients C P 1 , C P 2 , C D 1 , C D 2 used is the one proposed in Keehan Ye (2011). The resulting system is solved using an unsteady RANS (URANS) approach.
The PaSR model was benchmarked to the FGM model (Van Oijen and De Goey, 2000). In FGM, 1D flame configurations were employed by means of igniting mixing layer (IML) (Abtahizadeh et al., 2015) and premixed flamelet equations (Van Oijen et al., 2016). The thermochemical quantities were tabulated as functions of the mixture fraction and progress variable. At the CFD runtime, only the transport equations for the mean values and the variances of the controlling variables (in addition to continuity, momentum, and energy) were solved, and all the required parameters were retrieved. Specifically, the main assumptions behind the FGM method are that a turbulent flame can be represented by an ensemble of laminar flames, and the n-dimensional composition space can be replaced by a lower-dimensional manifold (Chen et al., 2018;Sorrentino et al., 2018a). The flamelet equations were computed through the specialized solver Chem1D, developed at the University of Technology of Eindhoven (TU/e) (Somers, 1994). A unity Lewis number assumption was chosen, whereas the GRI3.0 chemical mechanism was used.
Two controlling variables were chosen, a progress variable Y and a mixture fraction Z, to represent the mixing between fuel and oxidizer. Usually, under traditional combustion conditions, the standard FGM model adopts combustion products, such as CO 2 and H 2 O, to define the progress variable. Nevertheless, in this work, HO 2 was also included. Indeed, Medwell et al. (2008) stressed the role of precursor species like HO 2 in MILD combustion. Therefore, a combination of H 2 O, CO 2 , and HO 2 was selected as progress variable to include the effects of both preignition and oxidation chemistry. The general form of the progress variable is where X i is the molar fractions and α i the weight factors of the generic species i, which were optimized to yield a smooth mapping between the state space and the transported variables. In this work, the coefficients α were chosen as α H 2 O = 100/W H 2 O , α CO 2 = 100/W CO 2 , α HO 2 = 1,000/W HO 2 , where W i denotes the molar mass of species i. α i = 0 for all other species. Moreover, the progress variable was scaled between 0 and 1, to make it independent of the mixture fraction. The controlling variables and the lookup tables were defined by means of Matlab R2019a (The Mathworks Inc., 2019) scripts. Thus, all the thermochemical parameters calculated from flamelet equations were tabulated as a function of Y and Z.
Once defined, the laminar and adiabatic 2D tables were imported on Fluent using the FGM dedicated routine. The effect of the turbulence was considered incrementing the manifold dimension from 2 to 4, using a presumed β -PDF approach (Fiorina et al., 2005). In this way, both the mean values (100 table points) and the variances (11 table points) of Y and Z are considered as controlling variables.
During a simulation, the turbulent transport equations for the mean and variance of the mixture fraction (Z and Z υ , respectively) and for the mean and variance of the progress variable (Ỹ and Y υ , respectively) are solved in addition to the Favre-averaged mass, momentum, and enthalpy conservation equations. The four steady-state transport equations read as follows: In the equations above,ρ is the mean density,ũ j is the mean velocity vector, D is the molecular diffusion coefficient, µ t is the turbulent viscosity, Sc t is the turbulent Schmidt number, ω Y is the mean chemical source term of the progress variable, C 1 and C 2 are the modeling constants,k is the mean turbulent kinetic energy, andε is the mean dissipation rate. Since the lookup tables were created from adiabatic flamelets, heat losses were not included in the manifold and, therefore, only considered in the CFD simulations. The accuracy of this approach is then related to the accuracy associated with the inclusion of the enthalpy in Fluent once the adiabatic tables are imported. In particular, the software includes enthalpy by freezing the species in the adiabatic flamelets and recalculating the temperatures and the progress variable source terms with the enthalpy.

RESULTS AND DISCUSSION
Case 1 (described in Table 1) is first considered. Figure 2 shows the temperature profiles along the lateral thermocouple Frontiers in Mechanical Engineering | www.frontiersin.org obtained with the PaSR model in combination with three chemical mechanisms, the GRI3.0, GRI2.11, and KEE, using a C mix of 0.1. It can be observed that the temperature profiles are almost superimposed, with differences never exceeding 4%. Based on this observation, it was decided to keep the KEE mechanism for the subsequent simulations, to reduce the associated computational time.
In Figures 3A,B, the numerical predictions were compared to the experiments in terms of temperature profiles along the central and lateral thermocouples. It can be observed how the FGM method shows a good agreement with the experiments for the central thermocouple ( Figure 3A). In particular, the temperature profile predicted by FGM is quite homogeneous, and it well captures the temperature peaks in the near-wall regions, only slightly overestimating it by about 50 K. However, when looking at the lateral position, FGM strongly overpredicts the experimental trend of about 200 • , with a very localized peak close to the inlet zone (x < 0.06 m). Such a strong difference between the central and lateral positions suggests that the model predicts a very early ignition after injection, contrary to the experimental observations suggesting distributed reaction conditions. The fact that the reaction and heat release region predicted by FGM are very localized also explains why FGM well captures the central thermal field, the progress variable being unitary in this area and heat transfer dominating the flow features. The inability of FGM to capture the initiation and extension of the reaction region can be related to the absence of a controlling variable, taking into account the internal dilution by the combustion products, in the construction of the lookup table. Figures 3A,B also show the PaSR predictions for different values of the C mix constant. As a remainder, increasing C mix increases the mixing timescale value. The temperature profiles obtained with C mix values of 0.01, 0.1, and 0.5 clearly show the strong sensitivity of the model to the estimation of the mixing scale. In particular, decreasing C mix impacts the mean reaction rate in two ways, via the reacting structure fraction κ and the residence time in the reacting structures (taken equal to the mixing time). When C mix is reduced, representative of a condition of intense mixing, the reacting fraction κ (Equation 1) is increased, indicating that the whole computational cells evolve toward a perfectly stirred reactor. At the same time, the denominator in Equation (2) decreases, thus impacting directly the mean reaction rate, which is increased. A close observation of Figures 3A,B suggests that a C mix value of 0.5 is the most appropriate choice for the case under investigation, in agreement with previous results (Ferrarotti et al., 2018(Ferrarotti et al., , 2019. Nevertheless, when looking at the centerline thermocouple (Figure 3A), all PaSR simulations underpredict the temperature near the outlet zone (0.06 < x < 0.14 m), while overestimating the temperature Frontiers in Mechanical Engineering | www.frontiersin.org level on the sides. The shapes of the PaSR temperature profile clearly indicate that the model can capture the existence of a distributed reaction region, but not to the extent suggested by the experimental measurements, showing a very uniform profile. This is further confirmed by the analysis of Figure 3B, which shows how the PaSR results show a delayed ignition with respect to the experimental data, as the temperature peak reduces and moves downstream of the fuel injection. Increasing C mix helps to reduce the temperature overpredictions along the thermocouples, although the central thermocouple profile is still not quite well represented, as the maximum temperature is about 150 K higher than the experimental values. Opposed to experimental data, the PaSR simulation still retains two defined reaction zones, which are the cause of the temperature spikes along the central thermocouple profile (Figure 3A).
Interestingly, as C mix is decreased, FGM and PaSR become more similar. This is explained by the direct impact that the mixing timescale has on the mean reaction rate (Equation 2). As discussed above, a reduction in τ mix directly impactsω and increases it. This is further confirmed by the analysis of the Damköhler number maps presented in Figure 4, showing how the Damköhler in the domain decreases by one order of magnitude as C mix is decreased from 0.5 to 0.01. As Dα decreases, the reactive zone mass fraction κ (Equation 1) increases, thus increasing the mean reaction rate as well. The combined effects of reducing τ * and increasing κ in Equation 2 results in the above-mentioned temperature spikes in the C mix = 0.01 case.
One tool that can be helpful to understand the differences between the two models results is the so-called normalized flame index, introduced by Yamashita et al. (1996) and defined by Bray et al. (2005), Domingo et al. (2002), and Pitsch (2009, 2012) as By definition, the index takes the value of ξ = +1 in premixed flames and ξ = −1 in non-premixed ones. It can be then useful to assess the ability of the model to capture the premixed nature  of the combustion process under MILD conditions. Figure 5 shows the contours of the flame index obtained using both the PaSR model, using C mix = 0.5, and FGM. It can be observed that the two models yield quite different results. With the PaSR case, the non-premixed region is almost exclusively restricted to the injection zone, between x = 0 m and x = 0.05 m; in the rest of the chamber mixing, it occurs effectively before reaction, and the flame index is +1, indicating a premixed region. The same cannot be said about FGM: beside the injection zone, the flame index contour presents four additional regions close to the chamber walls, where the flame index value is −1, thus indicating a non-premixed flame structure. This can be ascribed to the very assumptions of the FGM formulation used in the present work and specifically to the lack of a dilution parameter to describe the mixing between fresh reactants and combustion products before reaction, thus resulting in the temperature overprediction shown in Figure 3B.
The PaSR results presented above clearly show that the model can potentially capture the distributed features of MILD combustion. At the same time, the estimation of the mixing timescale has a key role in the structure of the reaction zone and the temperature levels. In particular, it is clear that the use of a simple approach for τ mix , based on a fraction of the integral mixing scale, is not appropriate to describe the complexity of turbulence-chemistry interactions in the system under investigation. Therefore, the dynamic model (Equations 5, 6) was used for the same configuration, to assess whether the scalar time of a (non)reactive scalar could represent an appropriate choice. It should be stressed here that a large uncertainty is associated to the determination of this scalar dissipation rate via Equation 6 (Knudsen and Pitsch, 2009): the coefficients used in such an equation are the result of a calibration on simple-flow configurations (Ferrarotti et al., 2019) and might be inadequate for the complex flow under investigation. The results of the dynamic model using the coefficient by Keehan Ye (2011) show results very close to the C mix = 0.5 case. As shown in Figures 6A,B, the temperature profiles along the lateral thermocouple are almost identical, with the exception of the zone between x = 0.07 m and x = 0.15 m, where the dynamic model yields slightly higher temperatures and a flatter profile. The dynamic model yields slightly lower maximum temperatures in the reaction zone along the central thermocouple and lower temperatures around x = 0.1 m. Figure 7 shows the contour of the equivalent C mix eq for the dynamic case. The equivalent C mix eq is calculated as follows (Ferrarotti et al., 2019): where τ mix,d is the mixing timescale for the dynamic approach defined in section Modeling Tools. It is possible to see in Figure 7 how the average value of C mix eq in the reaction zone is close to 0.5, thus explaining the agreement of the static and dynamic cases. However, the dynamic model still cannot capture the initial ignition of the system. The reason for such deficiency can be explored by means of sophisticated models for turbulent mixing, using for instance detached eddy simulation (DES) and LES. The FGM and static PaSR with C mix = 0.5 were then tested for the experimental conditions of case 2 ( Table 1). With respect to case 1, both FGM and PaSR better capture the experimental trends. This is probably due to the dilution of the oxidizer stream, which causes the smoothing of the reaction process, resulting in a homogeneous temperature distribution. Figure 8 shows that, while FGM better predicts the temperature peak along the central thermocouple, the PaSR model enhances the temperature predictions in the ignition region along the lateral thermocouple, although it overpredicts the peak temperature by about 100 K. With respect to the central thermocouple, the PaSR model slightly overpredicts the peak temperature, while FGM underpredicts the observed values in the central region, that is, close to the outlet.

CONCLUSIONS
Numerical simulations of a lab-scale MILD cyclonic burner were carried out with two different turbulence-chemistry interaction closures, FGM and PaSR, in a RANS framework. The numerical results were compared to the available experimental temperature measurements, in the ignition zone and in the burner midplane.
The main findings can be summarized as follows: • The thermocouple profiles overall are better represented by the PaSR model results; while the FGM model yields slightly better results along the central thermocouple, the PaSR model better captures the typical slow ignition of MILD systems in both operative conditions. • A C mix of 0.5 seems to be the better fit for these combinations of fuel and operative conditions. The a posteriori analysis of the dynamic PaSR model equivalent C mix eq shows that for the static model, C mix = 0.5 is the optimal value for the present cases. • The analysis of the flame index reveals that the two models determine quite different mixing behaviors: FGM results show that the reaction regions all exhibit non-premixed behavior, while on the other hand the PaSR model determines mostly a premixed environment. As can be expected, the two different behaviors reflect on the temperature profiles, as the different ways mixing is handled between the two models strongly affect the computation of the reaction rates. Further experimental investigations of the cyclonic burner could then be key to determining which of the two numerical approaches would be the better option. • With the increase of the dilution level in case 2, the differences in the predictive capabilities of the two models lower; the dilution slows down the reaction process, yielding a more homogeneous temperature field, thus explaining the lack of a strong temperature overprediction in the reactive regions by the FGM model as in case 1.
Even though PaSR results are overall better than those of FGM, both closures display some limitations which can be ascribed to the turbulent mixing models employed in RANS computations. For this reason, LES analysis of the same burner will be carried out, to explore the effect of turbulent mixing modeling accuracy on the results.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
GC and GS provided experimental results and FGM results. RA and MF provided PaSR results and wrote first draft of the manuscript. GC wrote sections of the manuscript. All authors contributed to manuscript revision and approved the submitted version and contributed conception and design of the study.