Modeling Pollutant Emissions of Flameless Combustion With a Joint CFD and Chemical Reactor Network Approach

The Flameless Combustion (FC) regime has been pointed out as a promising combustion technique to lower the emissions of nitrogen oxides (NO x ) while maintaining low CO and soot emissions, as well as high efﬁciencies. However, its accurate modeling remains a challenge. The prediction of pollutant species, especially NO x , is affected by the usually low total values that require higher precision from computational tools, as well as the incorporation of relevant formation pathways within the overall reaction mechanism that are usually neglected. The present work explores a multiple step modeling approach to tackle these issues. Initially, a CFD solution with simpliﬁed chemistry is generated [both the Eddy Dissipation Model (EDM) as well as the Flamelet Generated Manifolds (FGM) approach are employed]. Subsequently, its computational cells are clustered to form ideal reactors by user-deﬁned criteria, and the resulting Chemical Reactor Network (CRN) is subsequently solved with a detailed chemical reaction mechanism. The capabilities of the clustering and CRN solving computational tool (AGNES—Automatic Generation of Networks for Emission Simulation) are explored with a test case related to FC. The test case is non-premixed burner based on jet mixing and fueled with CH 4 tested for various equivalence ratios. Results show that the prediction of CO emissions was improved signiﬁcantly with respect to the CFD solution and are in good agreement with the experimental data. As for the NO x emissions, the CRN results were capable of reproducing the non-monotonic behavior with equivalence ratio, which the CFD simulations could not capture. However, the agreement between experimental values and those predicted by CRN for NO x is not fully satisfactory. The clustering criteria employed to generate the CRNs from the CFD solutions were shown to affect the results to a great extent, pointing to future opportunities in improving the multi-step procedure and its application.


INTRODUCTION
Combustion of fossil fuels and industrial processes were estimated to contribute with around 65% of all anthropogenic greenhouse gases emissions in 2010, while being responsible for ∼85% of anthropogenic CO 2 emissions (IPCC, 2014). While there have been efforts to reduce emissions, the global CO 2 emissions have been increasing every year (International Energy Agency, 2018). Most scenarios developed for the future energy supply largely involve the combustion of biofuels, biomass, synthetic fuels, or hydrogen (Ellabban et al., 2014;Scarlat et al., 2015;Nastasi and Basso, 2016). Therefore, progress in combustion technology is necessary for the energy transition and also for the long term solutions that are being considered (Paltsev et al., 2018).

Motivation
Even though the major species of combustion are primarily dictated by the fuel that is being used, the emission of minor pollutant species is dictated by the combustion process. The emission of minor pollutants (CO, NO x , soot, and unburnt hydrocarbons) have a significant influence on local as well as the global environment. Carbon monoxide (CO) is a harmful substance usually present in flue gases when hydrocarbon or alcohol fuels are employed. Other typical pollutants are nitrogen oxides (NO x ). When air is used as an oxidizer, the formation of NO x tends to occur during combustion as the N 2 contained in air dissociates in the combustion reactions (Glarborg et al., 2018). The NO x emission from aircraft cruising in the upper layers of the troposphere or in the lower layers of the stratosphere is responsible for the formation of ozone (Grewe et al., 2012). At lower altitudes, the NO x emission is directly linked to respiratory diseases (World Health Organization, 2013), and is responsible for the acidification of water and rain that harms vegetation and wildlife (Camargo and Alonso, 2006). Therefore, there are tight regulations for the emission of CO and NO x and future restrictions on these emissions are being tightened, thereby posing challenges to designers of combustion systems.
One of the promising combustion technologies to lower pollutant emissions is the Flameless Combustion (FC) regime. Despite the fact that the FC regime does not have a clear and universally accepted definition, or that its physics are not yet fully understood (Perpignan et al., 2018a), it is clear that the regime has potential to substantially decrease emissions. Usually characterized by distributed reaction zones with lower temperature and species gradients, FC was shown to yield low NO x emissions while maintaining low CO (Kruse et al., 2015).
The modeling of FC is still challenging. Particularly, modeling tools able to predict pollutant emissions are required in order to assess, design and improve devices reliant on the regime. The use of Computational Fluid Dynamics (CFD) for modeling combustion largely relies on the use of simplified chemistry, as detailed chemical reaction mechanisms composed of hundreds of species and reactions cannot be accommodated at reasonable computational costs. Apart from the sheer number of reaction, solving detailed chemistry also adds the complexity of dealing with the wide range time-scales of the various chemical reactions, which causes the systems of equations to be stiff. Usually, the focus and strength of CFD has been on predicting reaction zone structures, temperatures and velocities, and not pollutant emissions.

CRNs for Emission Prediction
The lack of detailed chemical schemes has been pointed out as the main cause of discrepancy often seen between CFD and experimental results in emissions (Lyra and Cant, 2013). Because of that, the use of Chemical Reactor Networks (CRNs) to predict pollutant emissions has gained attention. A CRN is a set of ideal chemical reactors, in which detailed chemical reaction mechanisms can be applied at relatively low computational costs. On the other hand, most approaches developed for CRNs neglect turbulence-chemistry interaction and complex flow structures.
The design of CRNs manually has provided valuable results. The designs are usually based on experimental results or CFD simulations. The work of Lebedev et al. (2009) presented a CRN composed of six reactors to model a gas turbine combustor. The authors defined the CRN based on the mixture fraction field predicted via CFD. Likewise, utilizing velocity, temperature and species fields coming from CFD simulations, Park et al. (2013) developed a CRN to model a lean-premixed gas turbine combustor and obtained good matches with respect to experimental data.
More recently, Prakash et al. (2018) studied the effect of exhaust gas recirculation (EGR) on the emissions of a lean premixed gas turbine combustor using a manually designed CRN. The adoption of EGR might cause the required conditions for FC to occur, as O 2 concentrations drop and initial temperatures rise. Studying FC, the work of Perpignan et al. (2018b) showed the results of a CRN manually built based on CFD made with the Flamelet Generated Manifolds (FGM) approach, which is also employed in the present work. The authors achieved better agreement with CO and NO x experimental data using the CRN than with the reacting CFD. The focus was on a gas turbine combustor designed to operate under the FC regime. The detailed chemistry enabled the analysis of the NO x formation pathways, showing that the thermal pathway is not as important as the prompt, NNH and N 2 O pathways. This showcases the capability of in-depth chemistry analysis when utilizing CRNs.
The application of CRNs to predict emissions from FC systems has, in principle, advantages when compared to conventional combustion systems (Perpignan et al., 2018a). The highly distributed reaction zones, as well as lower gradients of species and temperatures, are better represented by ideal reactors (Cavaliere and de Joannon, 2004). However, this observation is valid for the largest scales, as there is evidence FC might be composed highly interactive flamelets at the smaller scales (Minamoto et al., 2013. The manual generation of CRNs is, nonetheless, particularly reliant on the experience of the designer and suffers from the lack of repeatability as a significant amount of trial and error approach is involved. Additionally, creating networks manually impedes the reproducibility of results and hampers systematic studies. As an attempt to avoid these issues, strategies to generate CRNs automatically based on CFD solutions have been developed. Early applications of automatic CRN generation were employed in simulations for analyzing furnaces (Benedetto et al., 2000;Faravelli et al., 2001;Falcitelli et al., 2002). Promising results were achieved as 3D CFD simulations were post-processed. The works utilized temperature and stoichiometry to define the ideal reactors, which could be PSRs or PFRs, based on the angle of the velocity vector. The underlying assumption was that the simplified chemistry models employed in the CFD were enough to predict temperature, velocity and major species. Moreover, these works provided insights into important variables to be taken into account when clustering CFD cells: temperature, a variable capturing the flow direction, and at least one measure of composition.
The same computational tool (with the same clustering criteria) was employed by Frassoldati et al. (2005) in the case of the swirling flame of the TECFLAM burner test case. The authors showed the effect of varying the number of reactors on the result of NO x emissions. For that specific case, having 300 reactors proved to be enough to guarantee no significant variation with a further increase in the number of reactors. Good agreement with the outlet value of NO x was shown, as only one operating condition was simulated. Fichet et al. (2010) presented another strategy to use a CFD solution as an input to build a CRN. The authors utilized a gas turbine combustor to showcase the capabilities of the approach. They opted for discarding the temperatures calculated with CFD and recalculated the temperatures based on the heat release of the chemical reactions predicted by the CRN. There was no systematic study on the advantages and disadvantages of this approach. The authors reported good agreement between measured and simulated NO x emissions. However, only one operating condition was available.
Similarly, Monaghan et al. (2012) utilized the Sandia Flame D as a test case and post-processed a CFD solution to build CRNs. The authors reported improved results for minor species (OH and NO) with respect to the initial CFD solution. The strategy to cluster CFD computational cells into reactors involved using mixture fraction, temperature and the axial coordinate.
In yet another example, Cuoci et al. (2013) utilized a few test cases to explore their computational tool, including the one utilized in the present work, presented in section The Test Case. The authors analyzed one of the operating conditions to which they reported improved NO x predictions, while CO was overpredicted to some extent. According to the authors, the source of deviation in NO was the incorrect temperatures predicted by the CFD simulations.
Despite the valuable developments, the full potential and limitations of using both CFD and CRN are not fully known. For example, the effect of different clustering criteria on the final solution is not clear, and that is one of the objectives of the present work. Additionally, there are other unknowns in the approach, as the effect of neglecting turbulence chemistryinteraction or taking it into account (via a PaSR-Partially Stirred Reactor approach, for example), or the pros and cons of solving the energy equation in the CRN step.
For these reasons, Automatic Generation of Networks for Emission Simulation (AGNES) (Sampat, 2018) was developed at the Delft University of Technology. The computational tool uses Cantera (Goodwin et al., 2018), an open-source software dedicated to chemical kinetics as a framework. Yousefian et al. (2017) reviewed the available hybrid computational tools based on CRNs for emissions predictions. In their evaluation of the different available solvers for CRNs, the authors highlighted two of the characteristics of Cantera that made it attractive for the development of AGNES: the capability of solving the energy equation and the fact that it is a free and open source solver, as further discussed in section Chemical Reactor Networks. Being able to modify and control the code was essential for the development of AGNES.
In the present paper, AGNES is utilized to simulate a test case related to the FC regime, in order to showcase its capabilities, limitations, and improvement opportunities. Additionally, the analysis of the results aims to clarify the reasons for the emission behavior in the chosen test case. The unique contribution of the present research is the assessment of different clustering criteria on the outlet emissions of an FC system that has a complex NO x behavior with varying equivalence ratio, as shown in section The Test Case.

THE TEST CASE
In order to evaluate the performance of AGNES, a test case was chosen based on the available information and the relevance to FC. Data availability on emission characteristics of the combustor for various operating conditions was a stringent requirement in making the selection. The test case described by Veríssimo et al. (2011) was chosen as a test case. The combustor was developed to study FC in a non-premixed combustion mode. The combustor consists of a cylindrical combustion chamber with a central air jet surrounded by 16 fuel jets in the burner head. This configuration is a variation of the most common jet-induced recirculation geometry in which a central fuel jet is surrounded by air jets (Hosseini and Wahid, 2013). The inlet air was preheated to 673.15 K, while the fuel (pure methane) was at room temperature.
The equivalence ratio was varied by maintaining the fuel mass flow (heat input of 10 kW) and altering the air mass flow. Data on OH * chemiluminescence, flue-gas temperatures and, more importantly, emissions of CO and NO x were acquired for all operating conditions. The equivalence ratio was varied from φ = 0.455-0.909 and the behavior of CO and NO x emissions was found to be non-monotonic. The authors reported a peak in NO x around φ = 0.53, while CO was practically undetectable for the runs at ∼φ = 0.53 and 0.59 (see Figure 1).
This behavior was explained by the authors based on a combined effect of the global equivalence ratio and the mixing characteristics of the combustor. Starting from condition a in Figure 1 and going toward leaner conditions, more air flow was added. The higher central air jet momentum was responsible, according to the authors, for quicker and stronger entrainment of the fuel, which had weaker jets. The stronger mixing of fresh reactants caused more intense reaction zones (as shown in Figure 1), while not allowing as much mixing of combustion products prior to the reactions. The result was then higher NO x and lower CO. This trend was dominant up to conditions e or f, in which the usual pattern of having lower NO x and higher CO as going leaner becomes dominant over the mixing characteristics.
In a more recent study, Zhou et al. (2017) performed OH and CH 2 O PLIF on the same setup, and referred to operating conditions a, c, and e as flameless, transition and conventional modes, respectively. They support that conditions e and f behave like conventional diffusion flames, as mixing with combustion products is apparently minimal and combustion is dictated by the mixing between fuel and air. As mixing with combustion products is allowed by the weaker entrainment of fuel by the air jet, the combustion regime transitions to FC. A summary of the operating conditions with φ, air mass flow rate and outlet O 2 concentration is shown in Table 1.
For conditions b and d, Veríssimo et al. (2011) provided pointmeasurements at various radial stations. Temperatures, major species (CO 2 and O 2 ) and pollutants (CO, UHC, and NO x ) were measured at 10 radial positions at each of the 10 stations.
The test case was simulated before, by Cuoci et al. (2013) and Lamouroux et al. (2014), for example. However, the focus of these previous works was on local values and a single operating condition for the former, and on the effect of including heat losses on local values for the latter. The focus of the present work is

COMPUTATIONAL MODELING
In an attempt to overcome the challenges in predicting emissions at affordable computational costs, a three-step approach is adopted in the current research: (1) solution of the flow-field with CFD using simplified chemistry and a turbulence-chemistry interaction model, (2) clustering of computational cells into ideal reactors based on criteria imposed to the CFD solution, and (3) solution of the generated CRN with detailed chemical reaction mechanisms. In the following subsections, the details of these steps are presented. The specific objectives of the performed modeling are: i. Evaluating the performance of the chosen CFD modeling ii. Evaluating the performance of the developed computational tool iii. Comparing the results obtained with CFD and CRNs iv. Obtaining minor species concentration more accurately than those obtained with the CFD simulations used as input to the CRN simulations v. Analyzing the NO x formation pathways for different operating conditions of the combustor

Computational Fluid Dynamics
The simulations were performed in order to assess the performance of CFD in predicting pollutant emissions and, more importantly, to generate the inputs for the subsequent modeling steps. The RANS (Reynolds Averaged Navier-Stokes) approach was adopted along two different turbulence-chemistry interaction models: EDM (Eddy Dissipation Model) (Magnussen and Hjertager, 1976) and FGM (van Oijen and De Goey, 2000). All simulations were performed in ANSYS Fluent R . The EDM was chosen to perform a comparison with FGM.
The assumption behind the model is that reactions are chemically fast and are controlled by turbulent mixing. Therefore, it was not expected that the EDM would perform accurately, at least not for all the operating conditions of the chosen test case. The EDM is possibly the simplest turbulence-chemistry interaction model and the objective of utilizing it was to assess how the quality of the CFD simulation utilized as input affects the results obtained solving the resulting CRNs.
A two-step reaction mechanism was adopted for CH 4 combustion, as shown by Equations (1) and (2). The EDM determines the reaction rates based on large turbulence time scales (k/ε), as shown in Equation (3) (Magnussen and Hjertager, 1976).
N prod. C prod.,r M prod. (3) The choice for FGM was based on its relatively low computational cost with respect to other models (Eddy Dissipation Concept, Conditional Source-term Estimation or transported-PDF, for example), and on its performance for various combustion systems (Verhoeven et al., 2012;van Oijen, 2018). This approach has shown to be promising for the modeling of FC, provided the progress and control variables are adequately chosen (Perpignan et al., 2018a). The approach allows the use of detailed chemistry in the pre-calculation generation of flamelets. Non-premixed flamelets were chosen due to the nature of the analyzed burner, and were solved using the GRI 3.0 chemical reaction mechanism (Smith et al.). Tests were performed utilizing the GRI 2.11 (Bowman et al.) and the POLIMI C1-C3 (Ranzi et al., 2012) mechanisms, but no significant differences were observed. The underlying assumption of the FGM model, or any flamelet-based approach for turbulent combustion, is that a turbulent flame can be represented as an ensemble of laminar flames. The flamelets were calculated in the mixture fraction space according to Equation (4), for species, and Equation (5), for temperature, according to the formulation of van Oijen and De Goey (2000). The pre-calculated flamelet quantities were then tabulated based on mixture fraction, a predefined progress variable and enthalpy, which are the variables that require transport equations. A presumed β shape PDF was employed to account for the turbulence-chemistry interaction, which was a function of the mean and variance of the mixture fraction.
The adopted progress variable as assumed to be dependent on the mass fraction of CO and CO 2 . were considered. Reaction rates of NO x species accounted for the temperature fluctuations by means of a β-PDF.
The closure for the RANS equations was achieved using the k-ε turbulence model. Since the burner has round air and fuel jets, the C ε1 constant was adjusted to 1.6 to correct for the wellknown round jet anomaly (Pope, 1978;Shih et al., 1995). Tests using a Reynolds Stress model did not provide superior results.
The modeling of heat loss is extremely important for the prediction of the temperature and the resulting emissions. The modeling of radiation was performed with the Discrete Ordinates model along with the weighted-sum-of-gray-gases approach to determine the required fluid properties. Heat conduction through the walls was imposed via wall temperature profiles. The profiles were determined partially based on the reported temperature values 5 mm from the walls for conditions b and d, as well as the outlet temperature. The profiles were first estimated by extrapolating the available data at the radial locations to the wall location. Subsequently, the resulting profile at the wall was multiplied by a factor to meet the outlet temperatures obtained in the experiments, thereby resulting in the same overall heat loss. For conditions in which local measurements were not available, linear interpolations and extrapolations of the profiles were performed, based on the equivalence ratio.
The computational mesh used was fully hexahedral and a 45 • sector was simulated, which included two fuel ports, as shown in Figure 2. The angle was adopted to guarantee a good mesh quality in terms of skewness. The mesh refinement was defined based on monitoring the outlet values of chemical species, as well as mid-plane averaged quantities, to guarantee no significant changes were attained with further refinement. The initial mesh size was ∼1 million elements and other four mesh sizes were tested until the difference in outlet species and averaged quantities was negligible. The final mesh was composed of ∼2.5 million elements.
Periodicity was imposed on the lateral boundaries. The first set of simulations was performed with mass flow inlets of air and fuel with uniform profiles. As the effect of having a developed flow velocity profile imposed as the boundary was shown to influence the results, this was adopted for the simulations herein reported. The developed flow velocity profile was assumed to follow Equation (6), a power-law velocity profile. The turbulence intensity was assumed to be 5% of the mean velocity. Tests with increased intensity did not significantly affect the results. Turbulence length scale at the boundaries was estimated based on the pipe diameters for both air and fuel, as it was imposed to be 7% of these dimensions. The outlet boundary condition was imposed to have zero static gauge pressure.

Chemical Reactor Networks
In order to simulate combustion with detailed chemistry at affordable computational costs, AGNES was employed. Having the results from the CFD simulations, the tool first clusters computational cells into reactors based on user-defined criteria. Clustering is executed with a Breadth First Search (BFS) algorithm used to traverse the computational domain. Such a domain traversal algorithm was chosen to ensure that the clustered cells are indeed connected to each other in the mesh, forming a continuous domain, and preventing the clustering of discontinuous pockets of cells satisfying the criteria. The clustering process proceeds until the total number of reactors (clusters) reaches the set-point imposed by the user, or until the specified maximum tolerance is attained. The quantity range value α is defined as shown in Equation (7), and it is dependent on a certain tolerance δ. The new local value, calculated by averaging the clustered cells, is allowed to deviate from the original local value as expressed by Equation (8). Several variables can be used simultaneously as criteria and every criterion must be satisfied for a cell to be included in a cluster.
α ≥ |c cluster − c cell | The resulting reactors have their properties assigned based on the properties of the CFD cells that compose the reactors. When maintaining temperatures as obtained from the CFD solution, which is the case for the present paper, averaged static temperatures are assigned, although total values are used during clustering to guarantee total enthalpy conservation. The mass flow exchanged between reactors is calculated based on the mass flow between cells in the CFD solution. There is an inherent mass imbalance due to the precision of the CFD solution which needs to be corrected to ensure consistent mass conservation. This is done by accounting for the mass flow between any two reactors as a fraction of the total outflow from the source reactor and assembling a system of equations. The boundary conditions of total inflow and outflow are accounted for in the system. The system of equations (mass, species, and energy conservation) is solved for, giving a vector of total outflow from each reactor and the "corrected" mass flow between reactors is recalculated using this vector and the matrix of outflow fractions. The mass exchange between reactors is also stored and maintained when solving the CRN simulation.
The system of ODEs is characterized by non-linearity and stiffness. The stiffness of the system is attributed to the wide range of time scales encountered for the reactions involved. Some reactions, such as those responsible for heat release, are much faster than those responsible for the formation of minor species such as that of NO x . This makes it difficult to integrate as the fast reactions require a small time step to be captured whereas the slow reactions need to be integrated over a larger time scale. Therefore, two levels of solvers are implemented, a local and a global. The local solver treats each reactor individually, whereas the global solver treats all the reactors simultaneously as a system. Cantera is used in the following ways: i. As a chemistry book keeping tool, in which Cantera ensures a consistent physical and chemical state of the reactors in the form of their temperature, density, and species mass fraction and storing the chemical reaction mechanism with all its thermodynamic properties. ii. It can also be used to solve a reactor network, in which Cantera calls SUNDIALS to perform ODE integration of a stiff system of equations, while treating it as a dense matrix.
At the local solver level, both features are used, whereas at the global solver level Cantera is used only as a chemistry book keeping tool, as shown in Figure 3. This is due to Cantera's solver limitations of being unable to simultaneously handle large number of reactors, hence the governing equations are written explicitly and solved using SciPy, a Python-based ecosystem for mathematics, science, and engineering, to solve a sparse matrix equation for the entire network of reactors. As previously mentioned, both a local and a global solver are employed. Locally, the solver available in Cantera is used to advance each reactor individually, changing its state and that of its connected reservoirs. This solver employs a different timestep to each reactor, based on its residence time. The global solver employs Cantera only to maintain the consistency of the chemical states, while all reactors are solved simultaneously with the sparse solver.
All reactors were considered to be PSRs. No turbulence fluctuation was taken into account in the CRNs. They were neglected because approaches with and without the inclusion of fluctuations are still being successfully employed (Yousefian et al., 2017) and a logical first step in the development of AGNES was to start without them. Moreover, this first application of AGNES is aimed at investigating the effect of clustering criteria on the solution. The inclusion of fluctuations may be important in some cases, as shown by Cuoci et al. (2013). For this reason AGNES will be modified to be capable of including fluctuations in the future.
The results herein presented were obtained by imposing a 3% tolerance. This resulted in different numbers of reactors for each test condition and imposed clustering criteria. Simulations were performed with both GRI 3.0 and GRI 2.11 mechanisms, to investigate the effect on NO x formation, as it has been shown that there are relevant differences in the NO x formation pathways under FC conditions (Perpignan et al., 2018b).
Several tests with various clustering criteria were performed in order to evaluate their efficacy to the present test case (Sampat, 2018). Apart from the obvious choice of temperature as a criterion, it was observed that the inclusion of velocity direction FIGURE 4 | Experimental and CFD temperature results along the centerline of the combustor (Above) and temperature contour plots (Below) for condition b (φ ∼ = 0.77) using FGM and EDM.
as a criterion was fundamental to capture the recirculation within the combustor, a key for capturing the chemistry within FC. Additionally, a variable acting as a tracer of the fuel, as well as a variable indicating the progress of reactions were required. In the present work, these are the Y CH4 and Y H2O , respectively. Two sets of criteria are explored, one with the aforementioned variables, and one with the inclusion of Y O2 .

RESULTS
In this section, the results obtained from CFD simulations are presented and discussed. The extent to which the CFD modeling was able to replicate the experimental data is shown, in terms of temperatures, major species and pollutant emissions. Subsequently, the results obtained with AGNES are thoroughly analyzed.

CFD Results
As expected, the FGM simulations showed better results when compared to the EDM, as shown in Figure 4. The temperature rise occurs earlier in FGM, better representing the experimental data points. It is worth noting the discrepancy between experimental and simulation values in the data point (z/D = 0.11). This difference might be attributed to the heating of the burner head and its piping, which might have caused the preheating of air to a temperature higher than that reported by the authors. Another hypothesis would be the effect of radiation on the thermocouple, although the authors assessed this effect. The error in temperature measurements was estimated by Veríssimo Although the centerline values of temperature shown in Figure 4 depict that FGM was superior, the extent of that is not well-represented. The peak temperatures attained in the combustor are not located at the centerline, but at the region where fuel and air mix. In this region, EDM peak temperatures are much higher than those predicted by FGM, as seen in the contours of Figure 4. Additionally, looking only at the centerline values shown in Figure 4 it may seem that reactions with FGM occur faster than with EDM, but that is not the case.
A major difficulty in simulating the case at hand is the uncertainties related to heat losses. Radiative heat losses play a role in the combustor, as well as conduction through the walls. The total heat loss can be derived by the reported values of outlet temperatures, although the exact location where the temperatures were measured is not fully clear. Therefore, a wall temperature profile was imposed based on the results available at 5 mm from the wall (r/D = 0.45), attempting to maintain them as close as possible. Temperature at this position is, however, not only dependent on the wall temperatures, but also on the reaction rates, recirculation, and radiation. The extent to which the simulations are able to reproduce these values is shown in Figure 5.
The centerline values of O 2 and CO 2 shown in Figure 6 indicate that predictions for condition d were more accurate than for the condition b. The reasons for this may be related to the fact that condition d is, according to Veríssimo et al. (2011) and Zhou et al. (2017), closer to a conventional diffusion flame, while condition b would be more representative of FC, making it more difficult to model, as a control variable representing vitiated recirculated gases may be required for the representation of FC with FGM (Huang et al., 2017). However, the values of CO 2 and O 2 close to the outlet should not be dependent on the combustion regime or the type of modeling. As also shown in Figure 7, simulations for condition b have higher CO 2 and lower O 2 values at the outlet. The discrepancy between simulations and experimental values on the O 2 values close to the outlet for condition b has been shown before in the results of Cuoci et al. (2013). The experimental uncertainty of 10% in the composition data does not justify the difference alone. Therefore, the uncertainty of the reported mass flows is probably causing the discrepancies.
The overall agreement between experimental data and the FGM simulations is good (Figures 7, 8). The largest deviations for temperatures occur closer to the burner head (z/D = 0.11), especially for condition b. Further downstream in the combustion chamber, the agreement is better. Simulations presented the overall characteristic of retaining lower temperatures near the centerline for longer axial distances, while combustion was more progressed (i.e., lower O 2 and higher CO 2 concentrations). This can be explained by a possible imperfect prediction of the mixing between burnt gases and reactants, as well as the effect that burnt gases have on the incoming reactants. Moreover, the aforementioned apparent inconsistency in temperature modeling (or measurements) close to the burner head, as well as the uncertainty related to mass flows might also have contributed.
The radial profiles (Figures 7, 8) show that the effects of the central jet spreading rate seem to be overpredicted in the simulations. This can be noticed by the fact that the gradients of the profiles tend to be lower than those obtained in the experiments for intermediate axial positions (z/D = 1.13 and 1.81). This was also the case in the reported previous computational works on the same test case (Cuoci et al., 2013;Lamouroux et al., 2014). However, the discrepancies could also be a result of poor prediction of reaction rates or heat transfer, instead of an artifact of the jet spreading prediction.
The predictions of pollutant emissions from the CFD modeling is poor, both in terms of the absolute values as well as in the trends (Figure 9). The emission of CO is highly overpredicted, in some conditions up to two orders of magnitude. The NO prediction is, on the other hand, underpredicted. The trend of NO is also not captured, as values increase monotonically with increasing φ. These results were to a certain extent expected, as the accurate prediction of CO with FGM is challenging (Ramaekers et al., 2010) and the NO x species modeling was not based on a detailed chemical reaction mechanism. Finally, the CFD results based on FGM performed reasonably well in order to be used as an input to the CRN simulations. Most of the variables have a good level of agreement (apart from pollutants). Better results could be expected by utilizing an FGM approach that is more suitable for modeling FC, such as the diluted-FGM, developed by Huang et al. (2017). Perhaps the mixing in the combustor has unsteady characteristics which can be captured adequately only by LES. However, previous attempts of simulating the test case with LES along with tabulated chemistry did not provide improved results (Lamouroux et al., 2014).

AGNES Results
The results obtained from the CRN simulations performed with AGNES are discussed herein with respect to the combustion FIGURE 9 | Outlet CO and NO results from CFD performed with FGM compared with experimental values. model obtained from the CFD input (FGM or EDM), to the clustering criteria, and chemical reaction mechanism.
The CO predictions are improved significantly with respect to the CFD simulations performed with FGM, regardless of the clustering criteria (Figure 10). The same is valid for simulations based on the EDM, in which the main difference with respect to FGM is in the predictions of CO for conditions closer to stoichiometry. All results predict the trend of CO and have fairly similar values. Additionally, no significant difference between GRI 3.0 and GRI 2.11 can be noticed.  The calculated NO x values based on the EDM simulations are overpredicted, as expected (Figure 11). The high peak temperatures attained with the modeling are responsible for the results being up to one order of magnitude higher. The results obtained from FGM are more complex to analyze (Figure 12). When temperature, velocity direction, Y CH4 and Y H2O are employed as clustering criteria, the NO x trend obtained with GRI 2.11 reaction mechanism is monotonic and its slope is the opposite of the experimental data: condition a had the highest NO x value. For GRI 3.0, overall values are higher than those of GRI 2.11 (a behavior previously reported by Perpignan et al., 2018b), but condition a has a lower value than condition b, making the overall trend non-monotonic. The values, however, are not in good agreement with experiments for both chemical reaction mechanism.
The inclusion of Y O2 as a clustering criterion is fundamental to attain a trend closer to experimental data with respect to the non-monotonic behavior of the emissions. The results using GRI 2.11 have values in the same order of magnitude as the experimental values (below 10 ppm). Such difference highlights that the choice of clustering criteria has an effect on the accuracy of the solution. The computational time can also be affected by the choice of criteria, as the required number of reactors changes. Additionally, the variables available from the CFD input can potentially influence the choice for a given CFD modeling approach, provided it has advantageous variables to be employed as clustering criteria.
One should bear in mind that all experimental values were below 10 ppm, and such low values pose a significant challenge for computational simulations. The improvement with respect to the values predicted directly via CFD is remarkable, and that enables the use of AGNES to aid the design of FC combustors. The peak of NO x emissions is, however, not predicted by any simulation. While experiments reported the highest values of NO x for the condition e, simulations had maximum values for the condition b, which are also closer to experimental values for the case of GRI 2.11 that included Y O2 as a clustering criterion.
The analysis of local NO x values for condition b presented in Figure 13 shows that there are significant discrepancies. Similarly to the results from CFD (Figure 7), the radial profiles have larger differences between central (r/D ∼ 0) and peripheral (r/D ∼ 0.5) locations than the experimental profiles for locations close to the burner head. Further downstream (z/D = 1.81 and 2.80), the opposite is true, as computations have flatter profiles than experiments. On the other hand, the centerline profile shows that the results are satisfactory and certainly represent an improvement with respect to the original CFD calculations. The overprediction of temperatures in locations close to the burner head (as shown in Figure 7) also play a role, as NO x values are increased at z/D = 0.11. The fact that NO x formation takes place in a relatively short axial region (around z/D= 1.13) and then stays relatively constant is remarkable, as seen in both simulations and experiments.
Additionally, some features should be further investigated. The fact that the GRI 2.11 can be more accurate raises doubts regarding the prompt NO x formation, as discussed by Perpignan and Rao (2019). The GRI 2.11 is supposedly inferior to its successor, GRI 3.0, but the apparent better representation of the prompt NO x pathway might be related to different reaction rates in vitiated environments. The chemistry of NO x formation was shown to be different in previous literature, not only because of thermal NO x abatement (Nicolle and Dagaut, 2006;Fortunato et al., 2018).
In order to understand what causes the non-monotonic behavior of NO x , the rates of formation by each pathway were estimated. The rates of the reactions responsible for NO formation at the end of each pathway were taken into account. It is important to highlight that the pathways are not independent and interact with each other. Therefore, isolating their respective contributions is not entirely possible. Nonetheless, the analysis herein performed serves as an indication. The largest source of uncertainty is concerning the thermal and prompt pathways. The N atoms that are products along the prompt pathway tend to react via the second reaction of the Zeldovich pathway (Lefebvre and Ballal, 2010). Therefore, the NO formation coming from the second and third reactions of the Zeldovich pathway was neglected. The thermal NO contribution is considered to be that of the first equation (O + N 2 ↔ NO + N), which is the rate limiting step of the pathway (Glarborg et al., 2018). As far as the prompt pathway is concerned, the GRI 2.11 mechanism considers the route to follow the HCN → CN → NCO → NO pathway, which was for long believed to be pathway through which the reactions progressed (Lefebvre and Ballal, 2010). Currently, it is known that the NCN route is active instead (Glarborg et al., 2018). In the present analysis, the prompt pathway is considered as it was in the development of GRI 2.11. Therefore, the NO forming reactions shown in Equations (9)-(13) were taken into account to calculate the NO formation rate.
In Figure 14, a comparison between the NO formation rates for conditions a, b, and e is displayed for the case with GRI 2.11 and Y O2 included as a clustering criterion. The figure displays the results along an axial line at the radial position r = 15 mm. This position was found to have the highest rates of formation and, therefore, it was selected. Several conclusions can be drawn from this analysis. The first important fact is that thermal NO decreases as φ is reduced (going from condition a-f ). In fact, the rates of thermal NO reduces as much as 3 orders of magnitude from condition a to f, showing that the contribution of thermal NO cannot explain the overall trend. Secondly, the role of the prompt pathway is shown to be important. The rates of NO production by prompt peak for condition b, which explains why this condition has the highest total value of NO x . Possibly, this peak was predicted for condition b in the simulations while it occurs for condition e according to the experimental data. The prompt pathway was previously shown to be responsible for the peak in NO at lean values for systems operating under FC or with high recirculation (Perpignan et al., 2018b;Perpignan and Rao, 2019).
Thirdly, the role of the NNH pathway is also prominent. Its peak NO formation value is higher at condition b than at condition a. It should be noted that reactions tend to occur further downstream for condition b if compared to the other two conditions. This is an indication of how important the predictions of jet development, recirculation and entrainment coming from the CFD solutions are, since NO formation via pathways other than thermal is dependent on slight variations in composition.

CONCLUSIONS AND RECOMMENDATIONS
The current work presented CFD and CRN simulations of a combustor developed to investigate FC. The FGM model and the EDM were employed in RANS simulations for simulating several operating conditions of the combustor. The CFD results were post-processed and CRNs were built via clustering of the mesh cells and solved with AGNES.
The following conclusions can be drawn from the study: • The RANS CFD simulations performed with FGM are able to replicate experimental data to a good extent. The main discrepancies were found near the burner head region. • CFD simulations performed slightly better for the condition in which a conventional combustion regime was attained (d) if compared to the condition with FC (b). • The CO emission predictions of all simulations performed with AGNES had a fairly good agreement. The proposed approach can achieve good predictions of CO even with computationally cheap and robust CFD modeling as the EDM. • The use of AGNES significantly improved the NO x predictions with respect to CFD, and a reasonable agreement with experiments was achieved. • Simulations showed that, remarkably, the NO x is formed in a relatively narrow region (from z/D = 1.5-2.0, approximately). • The variables chosen as clustering criteria proved to have a significant influence on the obtained results. For a given case, a certain set of criteria could prove to be necessary or optimal with respect to the number of reactors. To the best of our knowledge, this finding has not been previously presented in the literature. • NO formation in the combustor seems to be dictated by the prompt and NNH pathways. The variation of prompt NO x is responsible for the non-monotonic behavior of NO x with φ.
Further investigations should be carried out to improve or clarify the following issues: • The proposed method should be explored further in order to optimize clustering criteria.
• The reproduction of outlet emissions for the employed case may require the use of a CFD model able to perform both in the FC regime as well as in conventional combustion, as different operating conditions result in different combustion regimes. • The validation and subsequent use of chemical reaction mechanisms for highly vitiated conditions would probably aid the performance of AGNES. • The assessment on the effect of including turbulent fluctuations on the CRN is recommendable. This can be done by clustering the computational cells based on fluctuation values and/or by employing a Partially Stirred Reactor approach. • The prompt pathway should be investigated under FC and highly vitiated environments due its key role.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article.

AUTHOR CONTRIBUTIONS
AP conceived the idea behind the paper, performed CFD simulations, and took the lead in the writing of the paper. RS wrote the code for AGNES, performed CRN simulations, and actively contributed with the writing. AG supervised the performed worked, was fundamental in discussions about the methods and results, and actively contributed with the writing.

FUNDING
AP had his PhD research funded by CNPq (National Counsel of Technological and Scientific Development-Brazil) via the Science without Borders program.