Numerical Simulation of Laminar and Turbulent Methane/Air Flames Based on a DRG-Derived Skeletal Mechanism

Simulation of turbulent flames using detailed chemical mechanisms is still a challenge in numerical combustion due to the large number of species and the stiffness of the system of governing equations. In this sense, strategies to reduce the size of the detailed model are necessary and one of such models is the wellknown directed relation graph (DRG) method. In the present work, a DRGderived skeletal mechanism developed using only one application for methane/ air simulations is presented and validated for auto-ignition times, laminar flame speed and counterflow flames. The skeletal mechanism is tested for varying the equivalence ratio (φ = 0.4, to 3) and pressure (p = 1 to 150 atm). The temperature spans the range from T = 1000 K to T = 2000 K. The relative error, compared with the detailed mechanism, of our proposed model for ignition delay times and flame speed are less than 10% for most of the parameters. The skeletal mechanism is also used to simulate the piloted turbulent jet Sandia Flame D. Results show that this skeletal mechanism can reproduce the main features of laminar and turbulent methane/air flames. Article info Received: 17 November 2019 Received in revised form: 24 January 2020 Accepted: 2 March 2020


Introduction
The multidisciplinary character of combustion, which combines the interaction of thermodynamics, chemical kinetics, molecular transport and fluid dynamics, requires numerical simulations with a high computational cost. Reacting flows are modelled by a stiff system of partial differential equations [1,2], consisting of the conservation of mass, momentum and energy [3,4]. Mass diffusion and the effects of chemical reactions in the production/ consumption of each species must also be accounted through transport equations, in which the source terms have a high non-linear character and sensitive to any change of parameters, enhancing numerical instability. As the detailed kinetic mechanism may have hundreds or even thousands of species, the full set of equations for a detailed integration consists of hundreds/thousands of equations [5].
Due to the fluctuations and gas expansion, combustion generates instabilities in the flow field, leading to transition to turbulence [6] and thus, in most practical devices, combustion occurs under turbulent conditions [7]. In this case, the challenge is to describe the complicated coupling between chemistry and turbulence. In particular, the chemical source term is influenced by turbulent mixing [2].
Not only turbulence, but the large time scales, covering a range from 10 -10 s to 1 s [8] occurring in the chemical reactions turn simulations with detailed reaction mechanisms prohibitive for practical configurations, and frameworks for reducing the size and complexity of the kinetics are necessary. There are two main categories of reduction techniques for kinetic mechanisms: time scale analysis and the generation of skeletal mechanisms [9]. The latter consists in identifying the important and necessary species and in generating the mechanism only with those. Some examples are the directed relation graph (DRG) [10] and directed relation graph with error propagation (DRGEP) [11], sensitivity analysis based on Jacobian analysis [5] and also artificial neural networks (ANN) [12]. Time scale analysis is used primarily to identify a gap in time scales, so that the system dynamics can be described using only the slow time scales. Examples are the Quasi-Steady-State Assumption (QSSA) [13] and Partial Equilibrium Assumption (PEA) [14], Intrinsic Low-Dimensional Manifolds (ILDM) [15], Flame Prolongation of ILDM (FPI) [16] Reaction Diffusion Manifolds (REDIM) [17] [18], Computational Singular Perturbation (CSP) [19,20], Flamelet approach [21] and its developments, like Flamelet with Progress Variable (FPV) [22] and Flamelet Generated Manifolds (FGM) [23].
Among the skeletal reduction techniques, the directed relation graph (DRG) [10,24] consists in evaluating the error produced when one species is withdraw of the full mechanism. It has a simple theory and it is easy to implement, generating good results. Skeletal mechanism with DRG were developed for several fuels, such as n-heptane and iso-octane [24,25], ethylene [10] and n-decane [25]. Very large detailed mechanism were also reduced using DRG, such as those of biodiesel and biodiesel surrogates [26,27,1,28,29,30,31].
Yang et al. [32] developed a framework of dynamic adaptive chemistry (DAC) to efficiently simulate a turbulent flame, in which the DRG method is invoked for each cell of the computational domain to obtain a small skeletal mechanism that is valid for that local thermochemical condition. Sankaran et al. [33] conducted a direct numerical simulation of a 3D turbulent slot-burner Bunsen flame using a combination of reduction techniques: firstly, DRG is applied, following of sensitivity analysis and computational singular perturbation (CSP) in order to find the final version of the reduced mechanism. However, it is not very common to find applications of skeletal mechanisms resulting from DRG directly in the simulation of turbulent flames.
One of the main characteristics of DRG is its ability to be applied for a specified combustion application. This allows an optimized strategy, that can change the size and accuracy of the reduced model obtained. It has been pointed out by Lu and Law [10] that skeletal mechanisms generated for one type of application can also be used for others.
The present work has two main purposes: to carry on a numerical simulation of a turbulent flame using a DRG-derived skeletal mechanism and to show that this skeletal mechanism developed con-sidering one application, namely, auto-ignition, can also be used to model a turbulent flame. The test case is the well-known methane/air piloted turbulent jet Sandia Flame D [34,35]. Furthermore, validation is also performed for auto-ignition and laminar premixed/diffusion flames.

Brief outline of DRG methodology
The directed relation graph (DRG), developed by Lu and Law in 2005 [10], is a reduction method based on the construction of skeletal mechanisms. The aim of the method is to efficiently solve the species coupling, so that those who has little or none influence on the important species can be removed. DRG reduction is accomplished in a time that is proportional to the number of reactions in the detailed mechanism [24], is easy to implement and has become popular for very huge mechanisms [36].
The index r AB which accounts for the contribution of species B in the production/consumption rate of species A, can be quantified through where ν A,i is the stoichiometric coefficient of A in reaction i, n r is the number of reactions, is the source term (accounting for production/consumption of species) and (1) The index can be seen as a normalization value and defines an error estimate in predicting species A if B is neglected. Therefore, consumption and production reactions must be considered equally [11]. Defining and threshold value ϵ, if r AB > ϵ, then removing species B can induce error in the production of species A, so that species B must be maintained in the skeletal mechanism. One advantage of DRG is that different skeletal mechanisms with different levels of accuracy can be obtained depending on the user-specified threshold ϵ for the conditions under which it is developed [10]. The skeletal mechanism will converge to the detailed as ϵ is approaching zero, and the number of species can vary abruptly as the threshold is varied [2].
Eurasian Chemico-Technological Journal 22 (2020) 69-80 A very important choice on which results depend on is the set of species A for calculation of (1). These are called the target species and are selected as those who have large impact on certain chemical attributes (e.g. H 2 O 2 for high pressures) [11]. An acceptable choice for target species is a combination of fuel molecules, oxygen, main products of combustion or H-radical [37].
The source terms in index (1) depend on concentration and temperatures provided from combustion application simulations. Thus, depending on the final usage of the skeletal mechanism, those values can vary considerably. To obtain a mechanism over a sufficiently wide range of parameters (e.g. pressure, temperature, equivalence ratio and resident time), several values of temperature and concentrations are selected from typical applications, including perfectly stirred reactor and ignition, laminar flame propagation as well as counterflow flames [10]. Consequently, for each application, a skeletal sub-mechanism can be obtained, and the combination of them consists in the application-specific skeletal mechanism. Finally, the union of those mechanisms generates a desired skeletal mechanism that can be applied for the whole range of interest.

Implementation of DRG
A skeletal mechanism was already obtained using the above mentioned DRG strategy and has already been validated for ignition delay times [38]. The CH 4 detailed mechanism used is the San Diego (UCSD) mechanism [39], which consists of n s = 50 reacting species and n r = 247 elementary reactions. Although there are several detailed mechanisms for methane oxidation, where ignition delay-times differ from each other (see [38]), our choice for UCSD as detailed mechanism is restricted to the accuracy with the reduced model.
The set of target species is formed by reactants (CH 4 and O 2 ) and main products (CO 2 and H 2 O), which are considered usual targets species when applying DRG for several fuels [11,25,40,41,42]. With this set of target species we were able to reproduce necessary species, such as H 2 O 2 , important to accurately predict ignition delay times for higher pressures and C 2 H 6 , since for wider range of pressures and mixtures, the reaction where CH 3 reacts to itself to produce C 2 H 6 enhances the temperature necessary for ignition [43].
The detailed mechanism is used to simulate a zero-dimensional batch reactor to calculate the au-to-ignition time. For each time step, DRG is applied using temperature and concentrations from these calculations. Species with r AB > ϵ are retained and stored. The union of the species set from each time step of the reactor computations provides the final set of species. Only the reactions containing species from this set are selected to be in the skeletal mechanism, since our aim is to show that turbulent flames can be reproduced with a skeletal mechanism developed only with auto-ignition application. The readers are referred to [38] for a complete explanation about implementation and validation of the DRG algorithm applied in this work.
The resulting skeletal mechanism consists of 30 species and 131 reactions, obtained with ϵ = 0.11. This value of the threshold was chosen since the relative errors of ignition delay times compared with detailed model are less than 11% for the whole considered range. Also, the value ϵ = 0.11 is consistent with other literature work for methane mechanism reduction with DRG [44]. Nonetheless, others skeletal mechanisms can be constructed. A 19-species amongst 71 reactions was found as the minimal skeletal mechanism in which is possible to obtain reasonable results for low temperature auto-ignition calculations [38]. However, this model fails to reproduce high values of temperature and pressures, due to the absence of necessary radicals (see [38]), which would compromise its application on turbulent flames, the main goal of this paper.
The ignition delay time is defined as the point where there is the steepest gradient of OH concentration, which illustrates the exponential growth of radicals within the reactions. The simulations are calculated for varying the equivalence ratio (ϕ = 0.4, 1, 2, and 3) and pressure (p = 1, 50 and 150 atm). The temperature spans the range from T = 1000 K to T = 2000 K. The same configuration is applied to the detailed and skeletal mechanism and results are presented in Fig. 1. Figure 1 shows an excellent agreement between detailed and reduced model, which can be justified by the good choice for the initial conditions when applying the DRG algorithm, e.g., the set of target species and the user-specified threshold. Some deviations are observed for high pressures (p = 150 atm) and low temperature at the stoichiometric conditions (ϕ = 1). The errors in ignition delay times for atmospheric pressure and different equivalence ratio are all under 11%. For some conditions, errors are less than 0.1%. Further results for auto-ignition in the phase-space are presented in Fig. 2,   [39] (symbols) and skeletal mechanism (lines) in a 2D-projection of the phase-space, for a stoichiometric mixture at p = 50 atm and initial temperature T u = 1200 K.
for a stoichiometric mixture with p = 50 atm and initial temperature T u = 1200 K (values are given in specific mole fraction, i.e., ϕ i = w i /W i , where w i is the mass fraction and W i the molar mass of species i). Again, despite some deviations for H 2 O, the results agree very well with the detailed mechanism.

Turbulence modelling -RANS
Reynolds averaged Navier-Stokes (RANS) is a technique with acceptable computational cost in which the balance equations are averaged, and closure models are used to deal with the turbulence [3]. We use the RANS method to simulate the turbulent Sandia Flame D [34,35] within the open-source CFD software OpenFOAM [45] with the solver called ReactingFOAM, in which the pressure-velocity coupling is based on the PI-SO-algorithm [46] (pressure implicit with splitting operators). The RANS equations for continuity, momentum, species mass fractions and energy are given by [47,48] where u is the velocity, ρ the density, p the pressure, τ ij the viscous stress tensor, Y k , v k , are the mass fraction, the diffusion velocity and the reaction rate of species k, respectively, h is the sensible enthalpy, the heat released rate and λ the thermal conductivity coefficient. Here, denotes the variable Reynolds average and the Favre average.
We use the Boussinesq assumption for closing the equations, and the κ -ϵ model for the turbulent viscosity μ t [49]. In this model, two partial differential equations are solved (together with the conservation equations) for the flow turbulent kinetic energy and its dissipation , in order to obtain μ t as [49] (9) where C μ is a coefficient to be defined (see Eq. (9) below). The unclosed fluctuations terms for species mass fractions and enthalpy are modelled using the gradient approach as [49] (7) (8) where Pr t is the turbulent Prandtl number. The diffusion velocity v k is related to the diffusion coefficient, which is assumed equal for all species (Lewis number equal to unity). The numerical setup for simulation is consistent with the description in [48], and the discretization is based on the Finite Volume (FV) method.
For the interaction between chemistry and fluid flow, we use the partially stirred reactor (PaSR) approach. In this model, each computational cell is divided into a reacting and a non-reacting zone. The reacting zone is modelled as a perfectly stirred reactor (PSR) [3] and a parameter determining the mixing time scale, C mix , is required, which can be derived from the turbulent Reynolds number and the coefficient C μ as Constant C mix needs to be defined a priori, and for turbulent flows, this value can vary from 0.001-0.3 [50]. In this work, we use C mix = 0.3, the same value used in literature for this type of simulation [48]. The simplified PaSR model was chosen since it is already implemented in Open-FOAM, and therefore no modification to the code is necessary. Besides, the focus of this work is the reduced chemistry, and not the modeling of turbulence-chemistry interaction. In this way, the readers can easily reproduce the results with our skeletal mechanism (see supplementary material).

Numerical Results and Discussion
In this section, we present the numerical results for laminar and turbulent flames using the skeletal mechanism developed within the framework explained above. Flame speed of premixed flat flame and counterflow flame are simulated in order to validate the skeletal mechanism for laminar flames. These calculations are relevant in the context of validation as reactors and auto-ignition provide analysis for several values of temperature and pressure, covering ignition and extinction of the system for low and mean temperatures [51]. These analyses should be performed for different equivalence ratios, so rich, lean and stoichiometric mixtures can be analysed.
Nevertheless, it is important to extend the validation to simulations that involve the diffusive transport among the species that are not included in the reduction step. In this sense, the calculation of flame speed in premixed flames, through the modelling of a free flat flame, and simulation of counterflow diffusion flames must be employed. For the turbulent case, which consists in the goal of the present work, the well-known Sandia Flame D [34,35] is used as test case and comparisons are performed with experimental data.

Laminar premixed flame
The simulations of this section are performed with two different solvers. For the ignition-delay times and the premixed flame, we use the opensource software Cantera [52] and for the counterflow diffusion flame, we use the code INSFLA [53]. However, we have also tested the premixed flames cases in INSFLA and non-premixed cases in Cantera, and the same results can be obtained. The deviations between detailed and reduced mechanism are calculated as where f is the variable of interest. The solution of a freely propagation premixed flat flame with Cantera [52] follows the axisymmetric stagnation-flow equations [8]. The finite-difference method is used to discretize the flow equations and form a system of non-linear algebraic equations. A damped Newton method is them applied to solve the system. If the Newton iteration fails to find the steady-state solution, a pseudo-transient problem with a larger domain is attempted to be solved [52].
Laminar flame speeds for different values of equivalence ratio (ranging from ϕ = 0.5 until ϕ = 2.0) and temperature of unburnt mixture (T u = 298, 400, 500, 600, 700 and 800 K), for both the skeletal and detailed mechanism have been calculated. Figures 3 and 4 display the comparison between detailed and skeletal mechanism, and a good agreement can be observed. In Halter et al. [54], experimental data for flame speed is presented for a mixture of CH 4 /Air. Although in this work the pressure is approximately 10 atm (0.1 MPa), it consists a good validation for the skeletal mechanism developed in this work. Thus, in Fig. 3, we also compare the results with the experimental data [54], where a good agreement is also observed.
Although the maximum error calculated is 16%, for a very lean flame (ϕ = 0.5) and T u = 298 K, as can be seen in Fig. 5, for all other conditions the deviations remain small, all under 10%, which is an acceptable accuracy for the skeletal mechanism. For some conditions, errors are below 0.1%. It is important to emphasize that the skeletal mechanism was not developed using premixed flat flames as applications.

Laminar counterflow diffusion flame
Another simulation that was carried out to validate the skeletal mechanism, although it was not Fig. 3. Laminar flame speed calculated in a freely propagating flat flame at atmospheric pressure and T u = 298 K. Star symbols are results for detailed mechanism [39], cross symbols for experimental data [54] and lines for skeletal mechanism.
Eurasian Chemico-Technological Journal 22 (2020) 69-80 used in the DRG algorithm, is a counterflow laminar non-premixed flame. This is a one-dimensional flame consisting of two opposed flows, where a laminar flow leaves the fuel duct and stagnates against the flow leaving from the oxidizer duct.
One main characteristic of counterflow flames is that the flow velocity along the center line near the stagnation region varies linearly with distance [51]. Thus, the flow can be characterized by a single parameter, namely its velocity gradient a (unity: s -1 ), which constitutes the local strain rate. The strain rate can be used to perturb the counterflow flame. The higher a is chosen, the more the flame structure is perturbed by transport processes and its structure changes [55]. An extinguished flame can be obtained by enhancing the strain rate so that the reaction zone moves towards the wall, losing heat for it and consequently quenching, or by increasing both velocity and strain rate, so the mixing characteristic time scale becomes large than the reaction time scale.
The boundary conditions for the counterflow flame are the same one used for the Sandia Flame D, that is, 25% methane and 75% air at the left boundary, with temperature equals to T = 294 K and pure air at the right boundary, with temperature equals to T = 292 K. The pressure is p = 1 atm for both sides. Figure 6 shows the results of a coun- Fig. 6. Results of a counterflow flame, with strain rate a = 92 s -1 . Symbols are results for detailed mechanism and lines results for skeletal mechanism. terflow flame with strain rate a = 92 s -1 , which is a stable flame far from the extinction limit. It can be seen that the skeletal mechanism reproduces very well the full mechanism, for all the quantities considered. Figure 7 shows the results for a flame with strain rate a = 503 s -1 , a stable flame close to the quenching limit a = 651 s -1 . The skeletal mechanism can also reproduce this flame with great capability, which is an important feature when simulating turbulent flames, due to the range of local extinction and re-ignition of the latter. The main conclusion from these simulations is that the skeletal mechanism can reproduce with accuracy the existing diffusion processes in reactive flow simulations.
To further analyze the errors between detailed and skeletal mechanisms, we show in Fig. 8 the temperature and mass fractions of H 2 O, CO 2 and CO over different strain rates, with quantities calculated at stoichiometric mixture fraction in the counterflow flames. We can observe higher deviations for the mass fraction of CO, which can be associated with the simplification and removal of reactions in the skeletal mechanism, which modifies the reaction path of production and consumption of CO. We will see in subsection 4.3 that this difference can also be observed in the turbulent flame.
Thus, the results presented in this section show that the skeletal mechanism can reproduce the main characteristics of combustion processes, and the possible phenomena that can appear in a turbulent flame. The results of the skeletal mechanism are because it is not the smallest mechanism developed by DRG, but to the one with better accuracy for the selected applications. The consequence was that more species were retained to appear in the final reduced model, as explained in Section 2.2. The next section will show that this skeletal mechanism can also reproduce a turbulent flame with low degree of local extinction.

Turbulent case
The methane/air piloted turbulent jet Sandia Flame D [34,35] is used as experimental test-case to simulate the turbulent flame. This flame shows a small degree of local extinction, which allows a useful comparison for the reduced model presented in this work. The flame consists of a fuel jet with D = 7.2 mm diameter, 25% methane and 75% air. The coaxial pilot is a mixture of C 2 H 2 , air, CO 2 and N 2 , and is operated at lean condition (ϕ = 0.77), with the same nominal enthalpy and equilibrium composition as methane/air at this equivalence ratio. The Reynolds number of the jet is 22400, which represents a jet velocity of 49.6 m/s. In the present work, mesh generation and boundary and initial conditions are consistent to [48]. For the PSR model we also use C mix = 0.3 as in [48]. Be-Eurasian Chemico-Technological Journal 22 (2020) 69-80 sides, for simplicity, we do not consider radiation heat transfer. Figures 9 and 10 display the simulation results using both detailed and skeletal mechanisms for species profiles and temperature along the center line of the jet, compared with the experimental measurements. One can observe in Fig. 9 that the temperature is over predicted for x/D > 45, a difference justified by the radiation heat loss. Other methods, such as probability density functions (PDF) can also be used to improve this result [56].
Except for CO, Figs. 9 and 10 show that all species mass fractions have a good agreement between experimental, detailed and skeletal mechanisms, particularly the simulations results. The mass fraction of CO is under-predicted compared to detailed mechanism, although both simulations over-predict the measurements. However, this behavior is expected, since it is similar to the laminar case of counterflow diffusion flame presented in Fig. 8. As already mentioned, the difference is due to the elimination of reactions in the CO oxidation path.
We also compare in Figs. 9 and 10 the results of the skeletal mechanism against a numerical simulation using the chemistry modeling technique REDIM, that was used to simulate the Sandia Flame D using a RANS/transported PDF framework [56]. Unlike the DRG, the REDIM uses the thermodynamic and transport information of the detailed system to obtain a low-dimensional manifold imbedded in the state space that can be used to describe the full dynamics of the system. A full explanation of implementation can be found in [56].
One can observe in Figs. 9 and 10 that the RED-IM can describe the experiment with very good accuracy for all profiles for x/D < 60. In comparison with DRG, REDIM has a better global accuracy and in the regions where this is lost, the reason is due to the mixing model of the turbulence modeling. Results with the skeletal mechanism can be improved if a better method is used for the turbulence simulation.
The differences between the two methods are not only restricted to kinetics. Since DRG maintain the timescales of the original mechanism, the stiffness of the ODE system to be integrated is similar. This is not the case for REDIM, or other slow manifold-based methods [38], with which the system stiffness is indeed reduced. However, a different REDIM needs to be obtained for different simulations, while the skeletal mechanism has a more global character, with the advantage that it can be used for a whole range of combustion applications. As shown here, a skeletal mechanism developed using a not complicated homogeneous reactor application can be used to model the main characteristics of more complex systems, such as laminar and turbulent flames.

Conclusions
The main outcome of the results presented in the last sections is that the skeletal mechanism presented in this paper is universally applicable, i.e., can be used to reproduce the main features of combustion processes. Despite not using information about diffusion or extinction processes in its generation, the skeletal mechanism can be used to simulate a turbulent flame with low to moderate degree of local extinction.
The great advantage of using skeletal mechanism is the reduction of species mass transport equations to be integrated, which allows the achievement of reliable simulations in a time that is proportional to the number of species. Consequently, there is a decrease in the computational effort necessary for the simulations. Moreover, the derived skeletal mechanism can be used for different approaches and in other combustion applications.