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 (Pierce and Moin, 2001, 2004). In this context, tabulated chemistry techniques such as the flamelet-generated manifold (FGM) (Van Oijen and De Goey, 2000) and flamelet progress variable (FPV) (Pierce and 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.

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² 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², for the investigated conditions. In the second case, the air stream was diluted with 90% N₂ 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.

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 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″², 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, Pf=−2ρ¯uk″Z″˜(∂Z˜∂xk) is the production of scalar fluctuation, and Pk=−ρ¯uk″ui″˜(∂U˜i∂xk) is the production of turbulent kinetic energy. The set of coefficients C_P1, C_P2, C_D1, C_D2 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₂ and H₂O, to define the progress variable. Nevertheless, in this work, HO₂ was also included. Indeed, Medwell et al. (2008) stressed the role of precursor species like HO₂ in MILD combustion. Therefore, a combination of H₂O, CO₂, and HO₂ 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 α_H2O = 100/W_H2O, α_CO2 = 100/W_CO2, α_HO2 = 1,000/W_HO2, 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₁ and C₂ 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.

Case 1 (described in Table 1) is first considered. Figure 2 shows the temperature profiles along the lateral thermocouple 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, 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 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 Knudsen 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_mixeq for the dynamic case. The equivalent C_mixeq 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_mixeq 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.

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.

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.

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

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.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.