Methane/Air Auto-Ignition Based on Global Quasi-Linearization (GQL) and Directed Relation Graph (DRG): Implementation and Comparison

ABSTRACT In this work, two different reduction methods for simplification of chemical kinetics, the Global Quasi-linearization (GQL) and the Directed Relation Graph (DRG), are implemented and compared for ignition of CH /air homogeneous reacting system. The GQL approach is an automatic and scaling invariant method for a global analysis of the combustion system time-scale hierarchy that takes into account all the original species and elementary reactions, while the DRG is a reduction method for construction of skeletal mechanism that aims at removing unimportant species and thus elementary reactions from the original detailed kinetics. The auto-ignition process is used for validation in a wide range of system parameters and initial conditions which are important for gas turbines, dual-fuel diesel engines, etc. The present study shows that for methane combustion system studied in this work DRG requires dimensions at least two times higher than GQL to ensure the same level of accuracy for the whole range of considered system parameters. And, the advantages and disadvantages of both methods are discussed.


Introduction
Auto-ignition is a spontaneous ignition of a combustible mixture.This is an important topic in combustion science and engineering, because it is very sensitive to the chemical kinetics (Turányi and Tomlin 2016).Mathematically, auto-ignition is described by a system of ordinary differential equations (ODEs) using chemical reaction mechanisms (Warnatz et al. 1996).However, such detailed formulation results in two difficulties with numerical implementation.One is that typical chemical reaction mechanisms involve a large number of different elementary reactions between many species (Turányi and Tomlin 2016), leading to a high dimension of the ODE system.For instance, even for methane oxidation, up to a hundred species with of order of thousand elementary reactions (see e.g.Princeton mechanism; Shen et al. 2015;Zhao et al. 2017) should be accounted for a reliable description.
The other drawback is the high stiffness of the governing system of ODEs, which stems from different characteristic time-scales due to non-linearity of elementary reaction rates (Goussis and Maas 2011).Therefore, both difficulties cause an extreme complexity in numerical calculation and typically require large CPU times.For a general 2D or 3D reacting flows with complex geometry, such detailed chemical mechanisms remain computationally prohibitive (Goussis and Maas 2011;Tomlin, Turányi, Pilling 1997;Turányi and Tomlin 2016).Thus, model reduction should aim at both reducing the stiffness and number of the differential equations in the governing system of ODEs to be integrated into the numerical simulation.
The history of model reduction of chemical kinetics goes back to three conventional reduction methods, namely rate-determining step (RDS) (Campbell 1994;Dumesic 1999), quasi-steady state assumption (QSSA) (Bodenstein, 1913) and partial equilibrium assumption (PEA) (Chapman and Underhill 1913).In order to implement the assumptions applied in these methods, one must decide which species can be assumed being at quasi-steady state or which elementary reactions can be assumed equilibrated.Verification of these assumptions before the numerical experiments are conducted is difficult.Thus, only post-processing can verify the validity and accuracy, requiring a significant amount of time validation for human resources.The task becomes even more challenging when the assumptions are valid only for limited ranges of initial conditions or in a certain stage of system evolution (e.g. during an induction period of the ignition).It was shown in (Boivin 2011), for example, that even for the hydrogen/air system, one needs two different QSSA strategies to describe the auto-ignition processes for low and high initial temperatures.
Recently more advanced techniques have been developed.They can be mainly subdivided into two groups.One group deals with reducing the number of species and elementary reactions to obtain corresponding skeletal mechanisms.A very transparent representative of this group is the Directed Relation Graph (DRG) (Lu and Law 2005).The other group uses timescales and the concept of low-dimensional invariant manifolds, allowing to keep information about all species while reducing the dimension of chemical kinetics.Among them are the Intrinsic Low Dimensional Manifold (ILDM) method (Maas and Pope 1992b,a), the Method of Invariant Manifolds (MIM) approach (Gorban and Karlin 2003), the Global Quasi-linearization (GQL) (Bykov, Gol'dshtein, Maas 2008), the Minimal Entropy Production approach (Lebiedz 2004), the Computational Singular Perturbation (CSP) method (Lam and Goussis 1989, 1991, 1994), Rate-Controlled Constrained Equilibrium (RCCE) (Jones and Rigopoulos 2005) etc.Of course, there are lots of other methods which do not fall into these two groups.For example, Analytically Reduced Chemistries (ARC) (Felden 2017) and Optimized single-step chemistry models (Er-Raiy et al. 2018) are two methods aiming at constructing global reaction steps.
In this work, two model reduction approaches that belong to these two different groups, namely the DRG and the GQL, are investigated, compared and discussed for a homogeneous reacting system.The first approach focuses on finding out those species which have little influence on the important species and can be neglected (Lu and Law 2005).The GQL method is a manifold-based reduction method based on characteristic time-scale analysis, allowing a global identification of fast and slow invariant manifolds (Bykov, Gol'dshtein, Maas 2008).Both methods are implemented, validated and compared for CH 4 /air auto-ignition system.Note that although there are comparisons of both methods, we should emphasize that since both methods are developed and based on different concepts, the main target of this work is to show their own advantages and disadvantages.Furthermore, we only investigate the applicability of both methods in pure homogeneous system.The simulation of both methods for non-homogeneous reacting flows is beyond the scope of the current paper.
This work is also accompanied by a supplementary material, containing reduced chemistry we validate and compare in this work.

Mathematical model
In this work, we focus on auto-ignition processes under adiabatic and isobaric condition described by a pure homogeneous system of ODEs.In such case, the resulting mathematical model can be represented by the following system of ODEs (Turns et al. 1996): it can also be written in vector notation in the form: In the equations above, h denotes the enthalpy of the system, p the pressure, _ ω i the molar rate for the formation of a chemical species i due to chemical reaction in [mol/ (s Á m 3 )], where ϕ i its specific mole number defined by ϕ i ¼ w i =M i in [mol/g], w i is the mass fraction and M i is the mole mass of species i in [g/mol].This mathematical model describes the evolution of n s chemical species (n ¼ n s þ 2, with two additional variables describing thermodynamic properties such as h and p here) participating in n r elementary reactions.
As mentioned above, due to different characteristic time-scales which result in strong non-linearity of elementary reactions, the governing ODE systems Eq.(1) or Eq.( 2) is highly stiff.The characteristic time-scales can be obtained by determining the eigenvalues of the Jacobian matrix JðψÞ (Maas and Pope 1992b;Turányi and Tomlin 2016).Suppose the eigenvalues are sorted in the descending order as: where n 0 means the number of eigenvalues equal to zero, corresponding to the conserved quantities.Then, the characteristic time-scales at considered local point ψlocal can be obtained by calculating the inverse of each non-zero eigenvalues: τ ¼ 1=jλj.Usually in the reacting system, the time-scales differ largely, and the integration algorithm requires that the integration time step at this local point must be at most equal to the smallest characteristic time-scale 1=λ n .Therefore, quantitatively we can define the stiffness of ODE system Eq.( 2) at one local point as (Turányi and Tomlin 2016): The numerical simulation is performed using the in-house code HOMREA (Maas and Warnatz 1988), which is based on the semi-implicit extrapolation integration method (LIMEX) for the solution of a system of the ordinary differential equations (ODEs) and differential-algebraic equations (DAEs) (Deuflhard 1985;Deuflhard, Hairer, Zugck 1987).

DRG approach
The directed relation graph (DRG), developed by Lu and Law back in 2005 (Lu and Law 2005), is a reduction method to construct skeletal mechanisms.The aim of this method is to efficiently resolve the species coupling, such that those species having a little or none influence on a pre-defined list of target species can be removed.
The DRG is applied to identify the unimportant species.By using a numerical criterion the unimportant species are identified and elementary reactions associated with those species are removed.As it was originally suggested in (Lu and Law 2005), the contribution of species B in the production/consumption rate of species A can be quantified through the normalized index r AB , given by: where ν A;i is the stoichiometric coefficient of A in reaction i, _ ω i is the reaction rate and δ B;i is: if the ith elementary reaction involves species B; 0; otherwise: The terms in the denominator of Eq. ( 5) are the contribution of reactions to the production/consumption of species A, and the terms in the numerator are those from the denominator that involve species B (Lu and Law 2006a).The dependence of species A and B is compared with a small threshold value (0 < ( 1) (Lu and Law 2005), and if index r AB is bigger compared to this threshold, then removing species B can induce an error in the production of species A, so that species B must be maintained in the skeletal mechanism.Usually, species A are chosen as those who have some desirable chemical attributes that the reduced mechanism should reproduce (Lu and Law 2005).The resulting skeletal mechanism obtained has errors according to the user-specified threshold for the conditions under which it is developed (Pepiot-Desjardins and Pitsch 2008).Therefore, mechanisms with different levels of accuracy can be obtained by setting different values for .From a practical point of view, index r AB defines a direct error estimate in predicting species A if B is neglected.
As observed from Eq.( 5), species A in the index r AB should be specified.These A species are called target species and such choice is quite arbitrary.A common guideline for the choice of target species is to select those species who have large impact on certain chemical attributes (e.g.H 2 O 2 for high pressures) (Lu and Law 2005).In general, a suitable choice for target species is a combination of fuel molecules, oxygen and main products of combustion (Tosatto, Bennett, Smooke 2013), i.e. species whose accurate description is of interest and should also be provided.H-radical can also be chosen when dealing with fuels with large group of isomers which can be found in other literature (Luo et al. 2010).
To obtain a mechanism over a sufficiently wide range of parameters, such as pressures, temperature, equivalence ratios and resident times, a group of points in the parametric space are sampled for typical applications, including perfectly stirred reactor and ignition (Lu and Law 2005).Consequently, for each application, a skeletal sub-mechanism can be obtained for each point considered, and the union of these provides the applicationspecific skeletal mechanism.Finally, the union of those mechanisms generates a desired skeletal mechanism that can be applied for the whole range of interest.
Different formulations for the index r AB can also be defined.One can use the maximum over the summation in Eq. ( 5), since that the original formula proved to be inaccurate for large isomers groups in error estimation (Løv°as 2012;Luo et al. 2010) when used with, e.g.biodiesel mechanisms.According to Lu and Law (Lu and Law 2006b), there are also two more alternatives for r AB : the first considering the net rate of production for species A in the denominator of ( 5) and the second dealing with the forward and backward coefficients of each reversible reaction as two different reactions.While the first leads to only small differences in the results, the latter can produce completely invalid-reduced mechanisms.
Since in DRG the same level of importance is assumed for all species, without considering kinetic parameters when setting the coupling, some extensions were developed to overcome this issue, such as the DRG with error propagation (DRGEP) (Pepiot-Desjardins and Pitsch 2008), DRG with expert knowledge (DRGX) (Lu et al. 2011), DRG with sensitivity analysis (DRGASA) (Zheng, Lu, Law 2007) and DRGEP with sensitivity analysis (DRGEPSA) (Niemeyer, Sung, Raju 2010).These methods combine the framework of DRG with additional tools to improve the skeletal mechanism.For instance, DRGEP is also a graph-search-based method, however, its implementation differs from DRG.In DRGEP, a modified version of index ( 5) is used so that production and consumption are no longer considered equal (Pepiot-Desjardins and Pitsch 2008).Thus, one species depend on another based on the contribution to overall production or consumption rate (Niemeyer, Sung, Raju 2010), and the strategy allows a better selection of the chemical paths for a specific set of target species (Pepiot-Desjardins and Pitsch 2008).The DRGEPSA combines the framework of DRGEP and sensitivity analysis.The first phase is based on DRGEP, and the second phase uses a brute-force sensitivity analysis (chosen due to computational effort (Niemeyer, Sung, Raju 2010)) to those species that the index (5) are smaller but close to the specified .These are called limbo species, and they are first removed one-by-one from the mechanism to find the induced error in the simulation without that species.The limbo species are sorted in ascending order and removed from the mechanism, while the global error is evaluated after each removal (Niemeyer and Sung 2011).The advantage of this strategy was that DRGEPSA generates a smaller skeletal mechanism compared with DRG for n-heptane and iso-octane, with the same accuracy (Niemeyer, Sung, Raju 2010).In this work, as our main goal is to compare the standard DRG technique with detailed mechanism with moderate size, we use only the standard DRG.

Global quasi-linearization (GQL) approach
In contrast to the DRG approach that aims at reducing the number of species and constructing a skeletal mechanism, GQL keeps all information of species and elementary reactions in their original mechanism and aims at reducing the dimension through constructing an approximation of slow invariant manifold.It follows the main idea of the ILDM method about timescale separation that can be found by applying eigenvalue decomposition of the Jacobian matrix of the source term (see ( 1c)) (Maas andPope 1992b, 1992a).Thus, the decomposition into fast and slow subsystems and the fast and slow variables themselves is no longer empirical as in conventional approaches such as QSSA or PEA, but strictly reflect the dynamical properties of the considered system.From this point of view, singular perturbation theory (Skinner 2011) can be applied to identify these fast and slow subsystems (Bykov, Goldfarb, Gol'dshtein 2006).Therefore, the GQL approach and its global analysis of the system timescale hierarchy solve the problem with the selection of species and the approximation of the fast evolution of the system.
The GQL procedure is still based on several assumptions which need to be verified (Bykov, Gol'dshtein, Maas 2008;Bykov and Maas 2009), namely: • the system dynamics can be decomposed, i.e. there is a small parameter controlling difference in characteristic timescales, which determines the disparity of timescales.
• The vector field FðψÞ of the system (3) possesses a linear transformation of the original coordinate system, and it can be decomposed into sub-fields describing processes having different characteristic timescales.Thus, according to the singular perturbation theory (Bykov, Goldfarb, Gol'dshtein 2006;Skinner 2011), this linear transformation can be used to decompose the system into fast and slow sub-systems explicitly.
• The decomposition is valid everywhere inside our domain of interest.This domain of interest can be described as the domain with physically relevant values of system state variables for considered application range (e.g.positive concentrations, temperatures) (Bykov and Maas 2009).
As it is often observed for a typical system of combustion, due to the time-scale analysis and singular perturbation theory (Skinner 2011) these three assumptions are generally applicable for chemical reacting systems, and the following procedure allows us to find a global linear approximation of the vector field FðψÞ resulting in the dynamical decomposition and simultaneously verify these assumptions.If the second assumption is valid, then there is a matrix T GQL as a linear approximation in the domain of the vector field (1c): which has a simple geometrical interpretation mapping state vectors ψ to vectors FðψÞ and describes the linear transformation needed.Note that a linear and an accurate approximation of the source term are not possible, because the source term is strongly non-linear.And the main idea of the GQL approach is not to use the linearization form of source term to solve the ODE system, but to use this linearization approximation to define the subspaces suitable for decomposition and, thus, describing the slow manifold accurately enough (see details later).Then, under the assumption that the dynamical system has time scales with different orders of magnitude, this global linear transformation matrix T GQL can be used to identify and validate the assumptions about constant decomposition and to approximate relevant subspaces of fast and slow motions.The algorithm is the same as ILDM employing the Jacobian matrix of the system (Maas andPope 1992b, 1992a) with only one remarkable difference that the GQL approach is implemented to a constant matrix T GQL .By using, e.g. the Schur decomposition (Golub 1992) we obtain two groups of eigenvalues corresponding to different time scales: where Z s (n s Â m s ) are the invariant slow subspaces corresponding to the group of relatively "small" eigenvalues N s (m s Â m s ), and Z f (n s Â m f ) are the invariant fast subspaces corresponding to the group of relatively "large" eigenvalues containing eigenvalues with their absolute real parts ordered from the smallest to the largest.Thus, GQL and ILDM have the same algorithm using the timescale analysis and employing the Jacobian matrix to define the fast and slow subspace.However, for ILDM such algorithm must be applied at each state and the local Jacobian matrix is applied, while for GQL one needs to apply Schur decomposition for one time, and the resulted Jacobian matrix is globally valid.Note that although a local basis of eigenvectors for strongly non-linear system such as combustion changes drastically, e.g.before and after ignition, the subspaces still remain approximately the same with a high accuracy, which will be demonstrated later.
The singular perturbation hidden parameter can be estimated by the gap between the largest eigenvalue in group Λ s and the smallest eigenvalue in group Λ f (Bykov and Goldshtein 2013): Thus, we can transfer the original ODE systems into a singular-perturbated system by introducing new coordinates for state vector: here Ũ ¼ ðU 1 ; :::; U m s Þ T is the vector consisting of new variables describing slow motions and Ṽ ¼ ðV 1 ; :::; V m f Þ T the fast motions.Note that here the fast variables Ṽ defining actually also the quasi-steady state variables which are formulated in a linear combination of species.The fast subsystem is then given as: The zeroth-order slow subsystem reads: Equations ( 11) and ( 12) can be transformed back into its original coordinate system in the form of DAE systems, which are also implemented in HOMREA.The fast subsystem is described in original coordinate as follows (Yu, Bykov, Maas 2018): where And the slow subsystem is described in original coordinate in the form of DAE system as follows (Yu, Bykov, Maas 2018): where To summarize, the GQL approach provides a global decomposition that is valid for the whole domain of interest.Therefore, the applied subspaces Z s and Z f are constant, and the resulted Q f in Eq.( 13) and Q s in Eq.( 14) are constant matrix which can be used for all different conditions.
One important feature for the GQL approach, also for other manifold-based simplified chemistries, is that the fast process is decoupled from the dynamic system, leading to reduce the stiffness of the differential equations.In the GQL approach, fast processes are calculated by using the algebraic equations, and the slow processes are numerically integrated.Quantitatively we can also determine the eigenvalues of the integrated system: Then, the stiffness can be calculated as S ¼ jλ m s =λ 1 j, similar to Eq.( 4).
In the following, a simple example which is similar to (Lu and Law 2008) will be discussed.As shown in (Lu and Law 2008), non-quasi steady state (non-QSS) species can be misclassified for such system and advanced techniques must be used.Here we will show that the GQL approach can describe the global decomposition present in this system very accurately.The system consists of two steps involving three species A, B, and C: The time interval of interest here is 0 t 1.The rate constants are given as: where is a small parameter.The use of nonconstant k f and k r is to give a nonlinear dependence such as temperature dependence to the example.As explained in (Lu and Law 2008), although k f ) k due to the sufficiently small parameter , the species B cannot be classified as QSS species because its short species lifetime is induced by the second reaction step (B Ð C) which is actually a partial equilibrium reaction.However, the GQL can describe the correct fast variables and avoid such problem.Later in section 5.3, we will provide a more detailed explanation about the implementation of T GQL , here we only give the resulting GQL matrix.The system and the GQL invariant manifolds are implemented in a Matlab code which is also provided in the supplementary material.Following the algorithm in (Yu, Bykov, Maas 2018) we find the T GQL by evaluation the local Jacobian matrix of the corresponding source term at point ð½A; ½B; ½CÞ T ¼ ð0:3525; 1:9222; 0:9864Þ T .Then, the T GQL is given by Based on this GQL matrix, one can use the eigenvalue/eigenvector decomposition mentioned above to determine the fast and slow invariant manifold.In the left figure of Figure 1 we show the case for ¼ 0:1 ( 1 in the [A]-[B] projection.The red-dashed line represents different detailed solutions.The green line represents the fast motion predicted by GQL determined by Zs Á ψ ¼ Zs Á ψ0 (c.f.Eq.11), and the blue line the GQL slow manifold determined by Zf Á FðψÞ ¼ 0 (c.f.Eq.12).We notice that the GQL fast and slow manifolds can represent such system dynamics quite accurately.In the first phase (fast motion), the state moves along the fast manifold (green line) and is projected onto the slow manifold (intersection of fast and slow manifold).Then, in the second phase (slow motion), the state moves entlang the slow manifold (blue line) towards its final state (e.g.equilibrium point).
In the right figures of Figure 1 we show the time histories of species concentrations [B] and [C] for the case that at initial time t ¼ 0 the species concentrations are given as ½A 0 ¼ ½B 0 ¼ ½C 0 ¼ 1.We notice that during the early stage there is a deviation between detailed solution and reduced solution, and then the deviation becomes negligible.The deviation at the beginning stage occurs because the initial state is projected on the GQL slow manifold due to the fast motion, which can be inaccurate if is not small enough.For a sufficiently small this deviation also becomes smaller, and the GQL reduced solution has a much better agreement with detailed solution.Through this small example, we show how the GQL method performs in the case where the standard approach such as QSSA cannot be used.The suggested approach generically identifies an asymptotic limit and provides most appropriate basis for the system decomposition in this limit.

Results and discussion
There exist various chemical kinetic models for the CH 4 oxidation.In Figure 2, we compared the ignition delay times for Leeds (Hughes et al. 2001)  times can differ for about 20% at high temperatures (e.g. for Φ ¼ 0:4 and Φ ¼ 1:0) and for even more than 90% at low temperatures for all equivalence ratios.In this work, we restrict the accuracy of the reduced chemistry, such that the relative errors of ignition delay times over the whole range of interest is less than 6%, which is much smaller than the relative deviation between different detailed models of chemical kinetics and can be achieved by both DRG and GQL methods.

Validation
The CH 4 detailed mechanisms used in this work are the San Diego mechanism (UCSD 2014), which consists of n s ¼ 48 reacting species and n r ¼ 247 elementary reactions (without He and Ar).Therefore, the total dimension of Eq.( 2) is n ¼ n s þ 2 ¼ 50, consisting of n c ¼ 6 conserved quantities (two conserved quantities related to the enthalpy h and the pressure p, and four related to element conservation for atoms of C, O, H, and N), and the dimension of reacting subspace is n À n c ¼ 50 À 6 ¼ 44.The reactor simulation is performed varying equivalence ratios from 0.4 to 3.0, and initial pressures and temperatures are varied from 1 to 150 bar and 1000 to 2000 K, respectively.This range of interest corresponds to the studied conditions in (Li and Williams 2000), which is important in gas turbines applications.

Implementation and validation of DRG
The DRG algorithm is performed as follows: first, the target species, which will play the role as A species in the index r AB , are defined.The detailed mechanism is used to simulate an ideal gas reactor to catch the ignition time.For each time step, the DRG is applied using concentrations and temperature from the reactor calculations, and the species with the index greater than the threshold is retained and stored.The final set of species consists in the union of the species set from each time step of the reactor computations.The reactions containing only those species then are selected to be in the skeletal mechanism.This framework follows precisely the one explained in the original work of Lu and Law (2005).After the data from reactor simulations are obtained, the total computational time needed to generate the skeletal mechanism for the whole range of parameters is less than 1 min for using 2000 sample points.These values are based on a regular computer (intel i7 with 16GB RAM), within a serial framework.The set of target species is formed by fuel CH 4 and the oxygen O 2 , as well as the main products CO 2 and H 2 O.These are the usual targets species that are chosen when applying DRG for several fuels in other studies such as (Liang, Stevens, Farrell 2009;Niemeyer and Sung 2011;Niemeyer, Sung, Raju 2010;Pepiot-Desjardins and Pitsch 2008;Shi et al. 2010).It should be noted here that of course other species can also be considered as target species.This, however, leads to a larger skeletal mechanism and makes no obvious enhancement within the studied range.
Two skeletal mechanisms were obtained using the aforementioned strategy.One consists of 19 species and 71 reactions (DRG-19 chemistry) and the other 30 species and 131 reactions (DRG-30 chemistry).The only difference in the implementation is the threshold considered for the generation of each skeletal mechanism.Figure 3 shows the number of species in the skeletal mechanism as a function of .As the threshold approaches unity, the target species remain, but only those that appear together in the reactions are chosen.For values of > 0:5, species are strongly coupled with the target species and are retained in the final mechanism, since its elimination would lead to unrealistic results.Usually, a choice above > 0:2 can cause elimination of important species (Lu and Law 2006b), and this is where the strongly coupled groups of species are withdrawn together from the mechanism.Since the groups of species should be either kept or eliminated together, the choice of the threshold should consider these jumps (Lu and Law 2005).Thus, a conservative choice for the threshold is between 0:1 and 0:2.
In order to validate our implementation, we compare the results obtained for the reduction curve with the work of (Lu and Law 2008) in Figure 3.Although in that work the skeletal mechanism for methane is obtained using GRI 3.0 (without NOx chemistry) with 35 species as detailed mechanism, and the range of parameters covers different values (equivalence ratios from 0.5 to 1.5, pressure and temperature from 1 to 30 bar and 1000 to 1600 K, respectively), it can be used as an acceptable comparison.Both species labeled in Figure 3 (C 3 H 8 and CH 3 CHO) yield similar values of for elimination: for C 3 H 8 , ¼ 0:03 in this work and ¼ 0:02 in Lu and Law (2008), and for CH 3 CHO ¼ 0:14 in this work and ¼ 0:1 in Lu and Law (2008).Besides, for > 0:1, where the number of species become similar, there is the same pattern of decrease in both reductions.In Lu and Law (2008), a 30-species skeletal mechanism was obtained for ¼ 0:13, while in the present work, the same number of species were obtained with ¼ 0:11.
Thus, for the 19-species skeletal mechanism, we chose ¼ 0:21.Further increase of leads to an elimination of other species, which reduces the validity of the reduced model, producing errors above acceptable-specified accuracy.The 19-species skeletal mechanism is the minimal in which is possible to obtain reasonable results, e.g., low temperature in auto-ignition calculations and lean to stoichiometric mixtures (see Figure 4).However, this skeletal mechanism has some large deviations for high values of temperature and pressures and, therefore, we chose to expand to a large skeletal mechanism for better accuracy.
Decreasing the threshold in the implementation of DRG, for ¼ 0:11, a 30-species skeletal mechanism (DRG-30 chemistry) is generated to obtain more accurate results.It can be observed that at this stage a large number of species (e.g.CH 3 CHO, C 2 H 5 , CH 3 CH 2 O, etc.) was added to improve the mechanism.The choice for ¼ 0:11 is because C 2 H 6 is eliminated at ¼ 0:12 (see Figure 3), and C 2 H 6 is an important intermediate for rich mixtures.The reaction where CH 3 reacts with itself to produce C 2 H 6 enhances the temperature necessary for ignition (Reid, Robinson, Smith 1985).Therefore, the accuracy of results, specially for higher temperatures and pressures, is lost without this species.Furthermore, this value is consistent with others found in the literature concerning methane mechanism reduction (Lu and Law 2008;Luo, Lu, Liu 2011) In Figure 4, the results of ignition delay times based on DRG-19 chemistry (dashed lines) and DRG-30 chemistry (solid lines) are compared.Indeed, we observe clearly that DRG-30 chemistry performs better at higher temperatures.Figure 5 shows further comparison of the auto-ignition results from DRG-19 chemistry and DRG-30 chemistry at T 0 ¼ 1800 K, p 0 ¼ 50 bar and stoichiometric condition in the phase space.It can be observed that the DRG-30 chemistry has again a better agreement with the detailed mechanism than DRG-19 chemistry.

Implementation and validation of GQL
The detailed procedure for the implementation of GQL can be found (e.g. in Bykov and Maas 2009;Yu, Bykov, Maas 2018), here only main steps are briefly outlined for completeness.In general, the implementation of GQL can be performed in six steps: • Building up the T GQL ; • Calculation of the eigenvalues and eigenvectors of T GQL ; • Explicit decomposition of the system by using invariant subspaces; • Integration of the system based on GQL reduced chemistry using Eq.( 14); • Comparison of the solutions based on GQL reduced models to the detailed system solutions for varying initial conditions for validation purpose.
Among these steps, the most important and challenging step is to build up the T GQL .Several algorithms have been suggested, developed and validated.One possibility (Bykov, Gol'dshtein, Maas 2008;Bykov and Maas 2009) is to select n linearly independent vectors ½ψ 1 ; :::; ψn in such a way that the corresponding vector field ½ Fðψ 1 Þ; :::; Fðψ n Þ is also linear independent.
Then T GQL can be determined by using T GQL ¼ ½ Fðψ 1 Þ; :::; Fðψ n Þ Á ½ψ 1 ; :::; ψ1 À1 .However, this algorithm can cause degeneration and thus may lead to an artificial decomposition ( Bykov, Gol'dshtein, Maas 2008).To overcome this problem, in (Yu, Bykov, Maas 2018) another algorithm was proposed based on the mean value theorem.In this algorithm, the linear approximation of the original vector field is described by one local Jacobian matrix, and thus the local Jacobian matrix can represent the global behaviors of the system dynamics.In this case, T GQL ¼ Jðψ local Þ, where ψlocal is the state where one determines the local Jacobian matrix.In this work, the T GQL is determined by using the second approach, namely Furthermore, in the 4th step, the system is only integrated using Eq.( 14), while the equation for fast sub-system Eq.( 13) is not performed.This is because with the optimal dimension of GQL reduced chemistry, numerical simulations showed that the time for the initial states to relax towards the GQL slow manifolds lasts extremely short (in order of magnitude Oð10 À8 sÞ), and the initial states are already very close to the GQL slow manifolds.Therefore, integration of fast sub-system using Eq.( 13) is not necessary, and integration of the slow sub-system using Eq.( 14) is sufficient (see (Yu, Bykov, Maas 2018) for details).
The total CPU (also based on intel i7 and 16GB RAM) to for construct this 14D-GQL reduced chemistry takes 12 min.It takes 8.4 s on average to build up the T GQL and the consecutive Q s .The integration step and validation for different dimensions take 6.7 s on average for one initial condition and mixture composition.And this Q s has been also validated furthermore for 12 different initial conditions (from the lowest to the highest temperatures and pressures) and mixture compositions (from the leanest to the richest).
The above-mentioned procedure has been performed yielding a reduced model dimen- 6 shows typical convergence curves of relative errors of ignition delay times against dimension of GQL reduced chemistry for several initial conditions.It can be clearly observed that although for some cases (e.g.T 0 ¼ 2000K, p 0 ¼ 50bar and Φ ¼ 1:0, diamond-dashed line) a 13D GQL reduced chemistry is good enough to predict the ignition delay times with relative errors less than 6%, a 14D GQL reduced chemistry ensures the relative errors less than 6% for all cases.Hence, this 14-dimensional GQL chemistry has been applied, which is the minimal dimension one can obtain to ensure the accuracy (less than 6%) for the whole-studied parameter range.Figure 7 shows the comparison of ignition delay times between detailed (symbols) and 14D GQL reduced chemistry (lines).It can be observed that for all cases, GQL reduced chemistry can predict the ignition delay times very well.Although only several fuel/air equivalence ratios and pressures are shown in Figure 7, other conditions have also been tested, and all relative errors are less than 6%.

Comparison of DRG versus GQL
In Table 1, the reduced and remaining dimensions based on GQL and DRG are listed.
To compare the performance between DRG and GQL reduction method, we select DRG skeletal with 30 species (DRG-30 chemistry) because it guarantees the relative errors of ignition delay times to be less than 6% for the whole-considered range.Note that for DRG-30, although the skeletal mechanism consists of 30 reacting species, the number of remaining dimensions is 26 because n c ¼ 4 conserved quantities related to element conservation for atoms of C, O, H, N must be excluded.Besides, the dimension of skeletal mechanism is consistent with other found in the literature (Lu and Law 2008), although there the DRG skeletal mechanism is constructed based on GRI 3.0 (without NOx chemistry).It must be emphasized here that the GQL and DRG reduced chemistries derived here can guarantee the validity for the studied range (1000 K T 0 2000 K, 1 bar p 0 50 bar, 0.4 Φ 3.0).For both methods, one can also derive reduced chemistry with even lower dimension for certain specific application ranges.And the CPU cost to construct the reduced chemistry can be even smaller.
In Figure 8 we compare the species concentration (in specific mole number) in phase space for H 2 O, CO 2 , H 2 O 2 , and C 2 H 6 .We observe that for all these four species, results  other species can be better predicted by using GQL reduced chemistry.The fact that the GQL reduced model performs very well also demonstrates that the slow and fast subspace decomposition for optimal dimension does not change before and after ignition, although a local basis of eigenvectors can change drastically before and after ignition.As already mentioned, the high stiffness of the governing system of ODEs causes an extreme complexity in the numerical solution.Although this is still achievable for homogeneous system, it is still not realistic to use fully implicit solvers to integrate the system for 2D or 3D turbulent reacting flows.Therefore, a model reduction should also reduce the stiffness of the ODEs to be integrated into the numerical simulation.In Figure 10 we show one example for the stiffness of the integrated system of the ODEs along the ignition process for T 0 ¼ 1700 K, p 0 ¼ 1 bar and stoichiometric condition.Although only one example is shown here, stiffness under other conditions was checked and the same conclusion can be drawn.It can be seen that the stiffness of the system by using detailed chemistry (symbols) increases from order of magnitude Oð10 7 Þ dramatically to Oð10 9 Þ right after the ignition delay time.The GQL chemistry (solid line) causes the stiffness varying from Oð10 2 Þ to Oð10 6 Þ, which is at least a factor of 1000 lower than for the detailed chemistry.However, it is interesting to see that the stiffness by using the DRG-30 chemistry (dashed line) has a similar stiffness compared to the stiffness by using detailed chemistry.This can be attributed to the fact that the DRG approach only considers the importance of species and does not identify long chemical paths involving fast processes (Lu and Law 2006a).Therefore, the resulting skeletal mechanism preserves the characteristic time-scales.The reduction of dimension and stiffness of ODE system leads to a reduction of the computational CPU cost.Since DRG reduced chemistry reduces the number of ODEs (c.f.Table1) on one hand and preserves the stiffness on the other hand, its CPU cost is around 80% of the CPU cost by using detailed chemistry.For GQL reduced chemistry, it reduces both the dimension (c.f.Table1) and the stiffness (c.f. Figure 10) of the governing system of ODEs, its CPU cost is around 45% of the CPU cost by using detailed chemistry.
To summarize, both DRG and GQL have their own advantages and drawbacks.For DRG reduced chemistry, one obtains skeletal mechanism, which can be easily implemented in commercial software such as Chemkin or OpenFOAM.However, reduced dimensions and selection of target species are a crucial point which largely effect the accuracy.Moreover, because DRG is based on user-specified criterium for , those species that are important for the reaction systems but have small values of r AB may also be removed, which also results in inaccuracy.For GQL reduced chemistry, all species and elementary reactions are resolved based on the construction of an invariant low-dimensional manifolds.Thus, as already shown in this section, both the time evolution (e.g.ignition delay times) and species concentrations can be predicted with good accuracy.Moreover, the stiffness of the integrated system of ODEs using GQL approach is largely reduced.In practical application, the original thermokinetic states can be transformed into slow and fast variables.This allows to use special integration techniques such as inner iteration algorithms to speed up the simulations (Arvidsson, Løv°as, Mauss 2007;Chen and Tham 2008;Wang and Rogg 1993).

Conclusions
Two different approaches, the Direct Relation Graph (DRG) and the Global Quasi-linearization (GQL) methods, were outlined and discussed, and their performance was compared for the autoignition problem of CH 4 -air system.The results of the comparison were validated in a very wide range of initial conditions (1000 K T 0 2000K, 1 bar p 0 150 bar and 0.4 Φ 3.0), corresponding to relevant working load of gas turbine.
The DRG is a technique to generate skeletal mechanism, with an algorithm that is of easy implementation and can generate results very rapidly and can be applied to reproduce a single application combustion system or in a more global situation, covering several combustion processes.Thus, different sizes of skeletal mechanism can be developed.The main advantage of DRG is that the skeletal mechanisms constructed using homogeneous reactors as applications can be used to reproduce more complex system, specially involving diffusive transport processes.It can be observed that, for a methane/air ignition process, a small skeletal mechanism, consisting of only 19 species, is not enough to reproduce with prescribed accuracy the ignition process for high pressures and temperatures.A larger skeletal mechanism (here with 30 species) yields much better results.
The GQL method, on the other hand, considers all species and elementary reactions.The reduction has been achieved by finding out slow invariant manifolds given in an implicit form.For methane/air ignition process, the relative errors by using GQL reduced chemistry with 14 dimensions were less than 6% for the whole condisered range.
DRG and GQL are based on different concepts and were developed for different purposes.So we focus on a comparison of them only for two aspects.One is the stiffness of the ODE system to be integrated.It showed that reduced chemistry based on DRG has a stiffness similar to the stiffness based on detailed chemistry because the DRG skeletal mechanism still preserves the characteristic time-scales.The stiffness of the ODE system based on GQL reduced chemistry, however, is reduced (for the example considered here) at least by a factor of 1000 compared to the detailed chemistry.The reduction of stiffness (for GQL reduced model) together with the reduction of dimension (for both methods) causes for both methods a faster numerical calculation.The DRG skeletal mechanism (with 30 species remaining) needed around 80% of the CPU time when using detailed chemistry.The GQL reduced chemistry (with a dimension of 14) uses around 45% of the CPU cost compared to detailed chemistry.The other comparison focuses on the species in phase space.Since GQL can retain all species, it represents the species concentration in phase space better, compared with the DRG skeletal mechanism.Nevertheless, the DRG has the important advantages that it can be easily applied in other commercial computational software such as Chemkin or OpenFOAM.

Figure 3 .
Figure 3. Number of species in the skeletal mechanism against the threshold value for DRG.The gray lines indicate the species that are retained for each skeletal mechanism.

Figure 4 .
Figure 4. Auto-ignition calculation for the DRG-19 chemistry (dashed lines) and DRG-30 chemistry (solid lines) compared with detailed mechanism (symbols), for different pressures and equivalence ratio.

Figure 5 .
Figure5.System state space in 2D projection.Comparison of the auto-ignition results from DRG-19 chemistry (dashed lines), DRG-30 chemistry (solid lines) and detailed chemistry (symbols) at T 0 ¼ 1800 K, p 0 ¼ 50 bar and stoichiometric condition in the phase space.

Figure 6 .
Figure 6.Relative errors in % versus dimension of GQL reduced chemistry for several initial conditions.

Figure 7 .
Figure 7.Comparison of ignition delay times between detailed (symbols) and GQL reduced chemistry (lines).

Figure 8 .
Figure8.System state space in 2D projection.Comparison of the auto-ignition results from GQL chemistry (solid lines), DRG-30 chemistry (dashed lines) and detailed chemistry (symbols) at T 0 ¼ 1800 K, p 0 ¼ 50 bar and stoichiometric condition in the phase space.

Figure 9 .Figure 10 .
Figure9.System state space in 2D projection.Comparison of the auto-ignition results from GQL chemistry (solid lines), DRG-30 chemistry (dashed lines) and detailed chemistry (symbols) at T 0 ¼ 1800 K, p 0 ¼ 50 bar and stoichiometric condition in the phase space.

Table 1 .
Chemical reduction methods and their dimensions.