A Reduced Chemical Kinetic Mechanism for Toluene Reference Fuels Based On DRGEP and QSSA Methods

: As a gasoline surrogate fuel, the physical and chemical properties of toluene reference fuel (TRF) are relatively simple and stable, and the TRF chemical kinetic mechanism may be used in simulating combustion processes of gasoline. However, simulations using detailed or semi - detailed mechanisms have been limited due to the computational complexity and long computational time. For the construction of the reduced mechanism, the directed relation graph with error propagation (DRGEP) method is used to wipe out insignificant components efficiently, followed by the use of the quasi - steady state assumption (QSSA) method to separat e quasi - steady - state (QSS) species from the kinetic ODEs. In addition, some elementary reactions involving the formation and destruction of H and phenyl methyl radicals are subjected to sensitivity analysis and some kinetic parameters of the relevant elementary reactions are revised. As a result, a reduced mechanism involving 234 reactions and 60 species is developed. Comparing the experimental records with the analog data by utilizing the reduced mechanism, good agreement can be obtained when ignition delay time ( τ ), laminar flame speed ( S L ) and molar fraction of vital species are measured. Moreover, the mechanism may predict S L more accurately under lean mixture ( equivalence ratio φ < 1.0) conditions. The reduced mechanism is small and reliable in performance, which can commendably reproduce the combustion characteristics of gasoline surrogate.


INTRODUCTION
We can make full use of chemical kinetic mechanism to simulate the combustion process. However, simulations using detailed or semi-detailed mechanisms have been limited to the relatively simple fuels due to the computational complexity and long computational time. Meanwhile, in view of the limitation of time scale, stiffness caused by long time scale will directly affect the feasibility of chemical reaction mechanism. Therefore, effective reduction and revision of the elementary reaction parameters are necessary while retaining sufficient species and reaction steps to adequately reproduce relevant physical phenomena under target conditions. Gasoline, one of the conventional fuels for internal combustion engines, is a complex mixture that contains nalkane, iso-alkane, naphthene, olefins, aromatics and oxygenates. It is generally believed that the physical and chemical characteristics of gasoline can be characterized by a few single components, that is, gasoline surrogate fuel [1]. For instance, iso-octane, the simplest surrogate fuel, is often used to represent gasoline [2], coupling with multidimensioned computational fluid dynamics software. Dualistic mixture of iso-octane and n-heptane, called primary reference fuels (PRF), can be used to represent gasoline with a variable octane number. Pitz et al. [3] proposed that gasoline surrogate fuel should contain three essential components: n-heptane, iso-octane and toluene, for the reason that toluene is usually the richest aromatic compound in gasoline based on the importance of developing an accurate chemical kinetic mechanism for numerical simulations [4]. Therefore, the chemical kinetic model of three components, namely, iso-octane, n-heptane and toluene, can be considered the chemical mechanism of the gasoline and applied to simulate the combustion of gasoline. Recently, more and more researches have been carried out on the chemical kinetic mechanism of TRF fuels, shown in Tab. 1 and References [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. Common validation experiments for the chemical kinetic mechanism include flow reactor (FR), rapid press machine (RCM), jetstirred reactor (JSR), shock tube (ST), laminar flame (LF), and homogenous charge compression ignition (HCCI) engine.  [7,8] 49 83 CFR engine Zhang et al. [9] 70 169 ST, HCCI Liu et al. [10] 56 168 ST, LF, FT, JSR, HCCI Ra et al. [11] 113 487 ST, HCCI Andrae et al. [12] 137 633 ST, LF, HCCI Mehl et al. [13] 256 -ST, LF Chaos et al. [14] 469 1221 ST, FR, LF Anderlohr et al. [15] 536 3000 ST, RCM, JSR, HCCI Sakai et al. [16] 783 2883 ST, FR Andrae et al. [17] 1083 4635 ST, HCCI Heghes et al. [18] 1087 4639 ST, HCCI Andrae et al. [19] 1161 4961 ST Mehl et al. [20] 1550 6000 ST, RCM, JSR Andrae et al. [21] 137 633 ST, RCM A detailed TRF mechanism was proposed by Andrae et al. [17] and it was found that the action between toluene and n-heptane in this detailed TRF mechanism was not as important as the earlier conclusions published in 2005 [22]. Chaos et al. [14] developed a TRF kinetic model, which was verified by shock tubes and flow reactors. The comparison between the experimental records and the analog data showed that the direct interaction between the alkane and the toluene molecule was small. This conclusion was also the same as that of Andrae et al. [17], and denied the earlier conclusions of Andrae et al. [22]. Heghes et al. [18] put forward a detailed TRF model based on the model of Andrae et al. [17]. Among them, the submechanism of toluene followed the recommendations of Dagaut et al. [23], and the PRF mechanism was derived from Lawrence Livermore National Laboratory. At the same time, Sakai et al. [16] put forward a comprehensive and specific mechanism, including the interaction between alkanes, olefins and toluene. Mehl et al. [20] proposed the most detailed TRF kinetic mechanism. Then Bhattacharya et al. [2] carried out a simulation using the iso-octane mechanism involving 143 species and 643 elementary reactions in 2017.
Due to the diversity of components of real fuels and the complexity of chemical kinetics, it is time-consuming and costly to perform simulations and detailed chemical kinetics exceed the currently available computational abilities. Therefore, it is desirable to reduce the detailed mechanisms while maintaining their good performance. Kim et al. [5], Lee et al. [6] and Machrafi et al. [7,8] proposed the TRF skeleton models based on the PRF [24] and TRF mechanisms. Based on the PRF skeleton model, Zhang et al. [9] developed the TRF skeleton chemical kinetic mechanism by adding the oxidation mechanism of toluene, and it was verified by RCM, ST and HCCI. A methodology of reducing chemical kinetic mechanism for surrogate fuels, called as semi-decoupling, was proposed by Liu et al. [10,25].
In a previous study, Andrae et al. [21] presented a semi-detailed model (633 reactions and 137 species) for the oxidation of TRF. However, the model was not so compact that it was time-consuming and costly to perform simulations with this mechanism. Besides, there was a little mismatch for the calculated laminar flame speeds and experimental records under lean mixture (φ < 1.0) conditions [21]. Therefore, the aim of the paper is to put forward a new TRF chemical mechanism，which cannot only make the oxidation mechanism compact, but also accurately forecast auto-ignition delay time (τ), laminar flame speed (SL) and distribution of the vital substances.

MECHANISM REDUCTION METHODS
In this study, a reduced TRF mechanism based on the aforementioned model is developed by implementing a reduction scheme with the directed relation graph with error propagation (DRGEP) and quasi-steady state assumption (QSSA) methods. The principles for these methods are as follows.

DRGEP Method
The DRGEP method is first developed on the basis of the directed relation graph (DRG) method [26]. They have the same fundamental principle, that is, the key point is the target substances that depend on the route of others. For example, the route dependence of species A on B is evaluated using the normalized contribution of A to the formation rate of B (rAB). However, the two methods differ in the process of deciding whether to remove a substance from the mechanism. Suppose r AB , r AD and r BC are 0.05, 0.04 and 0.05, respectively. If a redundant substance needs to be removed using the DRG method, the removed material will be D. However, the error of removing D (0.04) is much greater than that of eliminating C (0.0025). Obviously, in order to make the mechanism more reasonable, the substance removed should be C. This is also the core idea of the DRGEP method. If the error caused by removing a substance is less than the user-defined threshold, the substance will be considered as a redundant substance and will not appear in the final skeleton mechanism. Pepiot-Desjardins et al. [26] proposed the normalized contribution of A to the formation rate of B, which is calculated by Eq. (1).

( )
In Eqs. (1) to (3), subscripts i and j represent the i th reaction and j-th species, respectively; subscripts f and b show the forward and reverse directions of the chemical reaction, respectively. R i is the algebraic rates of the i th reaction. C represents molarity. v  and v  are the reactant, product stoichiometric coefficients. K is the rate coefficient. n r,i , n p,i and n r represent the numbers of reactants, products and the total, respectively. However, Eq. (1) only considers the path dependence of the reactions including A and B. Then the Jacobian direct interaction coefficient is defined in Eq. (4). Following the round-off errors in reduction method [27,28], Chen et al. [29] improved the mechanism on multiple gasoline-ethanol surrogates using this approach.

( )
where J A,B is the Jacobian matrix element, is the maximum absolute data in the same set of elements. The semi-regularization sensitivity ratio J A, j is described as below.
where Y A and Y j represent the mass fraction of A and the j th species, respectively. W A and A w  represent the molecular weight and molar formation rate, respectively, ρ is density.

QSSA Method
Further mechanism reduction is necessary to be implemented when the DRGEP method achieves threshold.
The target search algorithm is a way: under specific target situations, it is an effective trial-and-error way to test redundant species and reactions by arranging the maximum mole grade continuously. After wiping out redundant species and a series of related reactions, the maximum relative errors between the current reduced mechanism and the initial mechanism are calculated at all conditions. The QSSA can be used for the construction of the reduced model [30], for instance, when the forward reaction rate and the reverse reaction rate are approximately the same, small changes in concentration can be ignored. In this method, the concentration of QSS species and others can be solved by the non-linear and differential equation systems, respectively.

IMPLEMENTATION OF THE MECHANISM REDUCTION
In this research, a reduced model for the base TRF mechanism is established by carrying out a reduction scheme with the DRGEP and QSSA methods. Moreover, some elementary reactions involving the formation and destruction of phenyl methyl and H radicals are subjected to sensitivity analysis and some kinetic parameters of the relevant elementary reactions are also revised. Then some important cross-reactions in the original mechanism are reserved in the final TRF mechanism.

Base TRF Mechanism
The semi-detailed mechanism (633 reactions and 137 species) by Andrae et al. [21] is applied to the reduction scheme. The mechanism was composed of three parts: the toluene subset, the iso-octane subset and the n-heptane subset. The toluene subset was established based on a benzene mechanism (Alzueta. 2000), which was suitable for flow reactor conditions with air fuel ratios ranging from near theoretical air fuel ratio to very high. The initial temperature (T) can be increased from 900 K to 1450 K and the residence time is approximately 150 microseconds. The entire mechanism was constructed by adding the mechanism of n-heptane (Peters, 2002) and iso-octane (Tanaka, 2003) to the toluene mechanism.

Case Settings and Mechanism Reduction
To sample the semi-detailed chemical mechanism, 24 relevant conditions, listed in Tab. 2, are chosen at different initial pressures (Ps), temperatures and equivalence ratios.
In the DRGEP-stage reduction, O2, CO, CO 2 and several fuel components (C 2 H 5 OH, i-C 8 H 18 , n-C 7 H 16 , C 6 H 5 CH 3 , CH 3 CHO) are selected as the starting species. The desired error tolerances for targeted variables are defined as follows: the threshold of mole-fractions for CO, CO 2 , O 2 and H radical are set to 0.1%, the relative tolerance of C2H5OH, i-C8H18, n-C7H16, C6H5CH3, CH3CHO are all set to 5 %. The first-stage reduction is implemented by setting the above threshold and using the DRGEP method three times. After the first-stage reducing, the mechanism includes 62 species and 244 reactions.
For the construction of the reduced model, the QSSA method is also applied in this part. After setting a threshold (7.5%) for the worst-case error, some substances (like IC4H8, C 5 H 3 ) are treated as quasi-steady state species. The construction of the second-stage reduced mechanism is implemented by removing QSS species. As a result, the second-stage reduced model with 60 species and 226 reactions is developed.

Mechanism Revision
An effective reduction process of this kinetic mechanism for the substitute of gasoline has continuously been carried out. In this section, the reduced model is verified by the widely used TRF fuels (Surrogates A and B) for more practical application. The composition and related parameters of the selected gasoline surrogate fuels are shown in Tab. 3. Laminar flame speeds (S L s) of the target gasoline surrogate fuels have been calculated by utilizing the second-stage reduced mechanism. Fig. 1 shows the comparison of S L between the calculated and the experimental results [31] at different equivalence ratios. The variation of S L with φ calculated for Surrogate A approximately equals the change for Surrogate B, whereas the calculated S L s are about 20% lower than the experimental values at low φs.
By analyzing the sensitivity of TRF chemical kinetic mechanism for S L , the elementary reactions, that significantly affect S L , can be identified. According to the flow rate A-factor sensitivity analysis, several reaction rate constants are changed in order to improve the predictions. The sensitivity analysis indicates that regardless of fuelrich or fuel-lean flames, the following two reactions remain highly sensitive.
O2 + H = H + OH (R108) CO + OH = H + CO 2 (R123) Because R108 is a branch reaction, any elementary reactions that produce H atom increase the chain branching rate of R108. In addition, R123 not only increases H atom but also produces a lot of heat during CO 2 production. As a result, these two reactions play important roles in combustion process. This finding has also been applied to the mechanism of a tri-component fuel that consisted of iso-octane, n-heptane and ethanol [32].

Figure 1
The comparison of laminar flame speed between the second stage mechanism data and the experimental results at different equivalence ratios Moreover, a slight mismatch may be observed on S L between the second-stage reduction model and the original one. One of the reasons is that second-stage reduction model cannot precisely predict the reactions associated with H radical. Therefore, the parameters of some important reactions associated with the H need to be revised.
The adjusted parameters refer to the detailed mechanism of iso-octane proposed by Mehl et al. [20]. It contains relatively complete reactions related to the H atom radical and it accurately reproduces the characteristics of the iso-octane flame, so some significant reactions that affect SL are extracted from the model and they are researched by using sensitivity analysis. For the sake of raising the exactitude of the new mechanism, some important cross reactions (R227~R234) in the original mechanism are reserved in the reduced TRF mechanism.
Modified elementary reactions and kinetic parameters after reducing can be checked in Tab. 4, which also shows data from other studies.  By comparison, a discrepancy is observed among the kinetic parameters obtained from different references. Ultimately, a mechanism with 60 species and 234 reactions is established. Tab. 5 shows an overview of all species.
Comparison of the experimental records and simulation results using our developed model can be seen in the following part.

MECHANISM VALIDATION
The effectiveness of chemical kinetic mechanism is based on whether the mechanism accurately predicts τ, SL and mole fractions of the vital species.

Auto-ignition Delay Time
τ is defined as the interval from the initial temperature (T) until the fuel/air mixture is raised to a given temperature. This definition is close to the time for the temperature inflection point, τ = t(dT/dt)max. It is well known that OH is often used as a mark for ignition of paraffin fuels, so that the moment at which the molar concentration of OH is sharply increased can be used as the ignition point. The calculations are performed with a closed homogeneous constant volume adiabatic reactor.
The calculations τ s by using the final reduction mechanism and the experimental results [34] in ST are compared in Fig. 2 and Fig. 3. In these two comparisons, Surrogate A and Surrogate B (Tab. 3) are still used. Fig. 2 and Fig. 3 present the variation of τ with the increase of T at the initial pressure (P) of 3.04 and 5.07 MPa, respectively, for the stoichiometric TRF/air mixture. As shown in these two figures, the predicted τs and the measured values are consistent within a certain temperature range (850-1150 K), although no experimental values are available at low temperature (T < 850 K).  Fig. 4 and Fig. 5 show the comparison of τ between the experimental data by Herzler et al. [35] with ST and numerical simulation data by using this reduction model. A toluene/n-heptane fuel (65% toluene, 35% n-heptane) (Tab. 3) is used in this part. Three different Ps (1 MPa, 3 MPa and 5 MPa) are selected. The reaction (R27) is found to be very sensitive under lean mixture (φ = 0.3) conditions. In order to get accurate predictions, the rate of the forward reaction has been adjusted in advance (Tab. 5). For φ = 1.0, in order to prevent model predictions from deteriorating, some cross-reactions are re-added. Generally, good agreement can be obtained from the two figures.

Figure 4
Comparison of auto-ignition delay times between the experimental data by Herzler et al. [35] in shock tube and simulation data by using this reduced model. Experimental data normalized with τ exp = f(T) p −0.883 Figure 5 Comparison of auto-ignition delay times between the experimental data by Herzler et al. [35] in shock tube and simulation data by using this reduced model. Experimental data normalized with τ exp = f(T) p −1.06

Laminar Flame Speed
S L is a parameter for a given combustible mixture. It represents the basic diffusion, reaction and exothermic properties and it is often used in the verification of chemical kinetic mechanism. Furthermore, S L plays an important role in engine combustion. Thus, our model can be validated by comparing the simulated laminar flame speed with the experimental laminar flame velocity in the real gasoline/air mixture [31]. The calculations are performed in the CHEMKIN flame speed calculator reactor model. Fig. 6 shows the comparison between the calculated SLs and experimental values [36] for toluene at room temperature and 0.1 MPa. As mentioned earlier, the gap may be observed between the calculated and experiment results at low temperatures, but the error is within the tolerable range (the relative error ≤ 20%).
The next step is to compare the calculated S L s with experimental data for TRF fuels.
The experimental S L s that were measured by Zhao Z et al. [31] and the calculated results by using the reduction kinetic model are shown in Fig. 7 and Fig. 8. In addition, in order to compare the performance of the original mechanism and the present mechanism, the calculated values of the original model [21] are also introduced. It is found that the prediction accuracy of the final mechanism is not worse than the original mechanism at the conditions of T being 353 K and 500 K, and 0.5 ≤ φ ≤ 1.4. Besides, the present mechanism may predict SL more accurately under rich mixture (φ > 1.0) conditions.  Comparing Fig. 7 with Fig. 1, it could be found that the final kinetic model can sufficiently predict S L of the three-component fuel after the revision of reaction parameters. This also illustrates that S L is sensitive to the molecular structure of the retained species. Fig. 9 compares the mole ratio of phenol (C 6 H 5 OH) and 2.4-cyclopentadiene-1-one (C 5 H 4 O) between semidetailed mechanism [21] and this reduced mechanism at two different HCCI operating conditions: operating condition 1 (φ = 0.2857, P = 0.1 MPa, T = 393 K, engine speed (n = 900 r/min) and operating condition 2 (φ = 0.25, P = 0.2 MPa, T = 313 K, n = 900 r/min). The toluene/nheptane fuel (65% toluene and 35% n-heptane) (Tab. 3) is used. As can be seen from Fig. 9, our developed model does a relatively good job in predicting the mole ratios of phenol (C6H5OH) and 2.4-cyclopentadiene-1-one (C 5 H 4 O) of toluene/n-heptane fuel at two different HCCI operating conditions.  Figure 9 Comparison of mole fraction for phenol (C6H5OH) and 2.4cyclopentadiene-1-one (C5H4O) between semi-detailed mechanism [21] and this reduced mechanism in HCCI operating condition 1 (φ = 0.2857, P = 0.1 MPa, T = 393 K, n = 900 r/min) and operating condition 2 (φ = 0.25, P = 0.2 MPa, T = 313 K, n = 900 r/min) Above all, it can be found that the optimized mechanism is predictive in terms of auto-ignition phasing by comparing the predicted τ, S L and mole fractions of some vital species in the optimized mechanism with available experimental data at the selected conditions.

CONCLUSIONS
Starting from a semi-detailed mechanism (137 species and 633 reactions) for a gasoline surrogate fuel comprised of iso-octane/n-heptane/toluene, this study has developed a reduced TRF combustion mechanism, composed of 60 species and 234 reactions. The major results and findings can be summarized as follows: (1) The developed scheme includes a three-stage procedure. On the first stage, the DRGEP method is used to efficiently remove redundant species. On the second stage, the species with short timescales (QSS species) can reach the chemical equilibrium quickly so the total production rate is approximately zero, as a result the differential equations can be replaced by algebraic ones. On the third stage, the kinetic parameters of some elementary reactions involving the formation and destruction of Handphenylmethyl radicals are revised.
(2) The mechanism is sensitive to molecular structure, the formation and destruction of H radical can directly affect the characteristics of the mixture flame.
(3) In the HCCI simulation of the isooctane / n-heptane / toluene fuel mixture, Toluene (C6H5CH3) and Phenyl methyl radical (C6H5CH2) showed a strong negative sensitivity to the ignition delay under lean conditions.
(4) Adjustment of parameters for reactions involving the formation and destruction of Phenyl methyl and H radicals and retention of cross reactions are crucial to the accuracy of the simulation results.
(5) Good agreement is observed by comparing the predicted auto-ignition delay of the optimized mechanism with available ST experimental data. SL of real gasoline and toluene could also be predicted accurately by the optimized mechanism. Moreover, the model accurately predicts the change in the molar fraction of phenol (C 6 H 5 OH) and 2.4cyclopentadiene-1-one (C 5 H 4 O).
The simulation results of this model are in good agreement with experimental observations. This study confirms that the optimized mechanism would be very useful in engine simulations.

SL
laminar flame speed / cm/s τ auto-ignition delay time / ms φ equivalence ratio T initial temperature / K P initial pressure / MPa n engine speed / r/min