A Comparison of EGR Correction Factor Models Based on SI Engine Data

The article compares the accuracy of different exhaust gas recirculation (EGR) correction factor models under engine conditions. The effect of EGR on the laminar burning velocity of a EURO VI E10 specification gasoline (10% Ethanol content by volume) has been back calculated from engine pressure trace data, using the Leeds University Spark Ignition Engine Data Analysis (LUSIEDA) reverse thermodynamic code. The engine pressure data ranges from 5% to 25% EGR (by mass) with the running conditions, such as spark advance and pressure at intake valve closure, changed to maintain a constant engine load of 0.79 MPa gross mean effective pressure (GMEP). Based on the experimental data, a correlation is suggested on how the laminar burning velocity reduces with increasing EGR mass fraction. This correlation, together with existing models, was then implemented into the quasi-dimensional Leeds University Spark Ignition Engine (LUSIE) predictive engine code and resulting predictions are compared against measurements. It was found that the new correlation is in good agreement with experimental data for a diluent range of 5%-25%, providing the best fit for both engine loads investigated, whereas existing models tend to overpredict the reduction of burning velocity due to EGR.


Introduction
G overnment legislation on controlling emissions from combustion engine powered vehicles is becoming ever more stringent. Therefore, car manufacturers have to keep improving their engines to minimise fuel consumption and emissions. The latest European regulations (EURO VI) came into effect in 2014 and focus primarily on reducing nitrogen oxide (NO x ). An effective strategy to reduce NO x emissions is to use Exhaust Gas Recirculation (EGR), that is, to recirculate cooled exhaust gas into the cylinder. High levels of EGR have been found to reduce not only NO x emissions but carbon monoxide (CO), particulate mass (PM) and particulate number (PN) [1,2,3]. EGR has also been found to reduce the chances of auto-ignition due to lower end gas temperatures [4,5,6,7,8]. The reduction in end gas temperature also allows the compression ratio to be increased, leading to improved engine efficiency [9].
The fact that cooled EGR leads to a decrease in laminar burning velocity is well known. Therefore, in computer simulations of combustion engines, the effect of EGR is typically modelled through a correction factor that reduces laminar burning velocity depending on the applied level of EGR. A variety of different models have been suggested, but none of them are based on engine-like conditions. A summary of the conditions at which published correlation were determined is presented in Table 1.
The correction factor proposed by Metghalchi and Keck [10] was found experimentally, using a constant volume combustion chamber at an initial pressure of 0.1 MPa, for a stoichiometric mixture of isooctane and air for an unburned gas temperature range of 340-440K. The composition of diluent to simulate combustion products was 15% carbon dioxide and 85% nitrogen by volume.
Their model was later modified by Rhodes and Keck [11] for a blended fuel similar to gasoline, known as indolene. Because of the change in mixture, the composition of the diluent was altered to simulate combustion products, with the diluent now containing 20% carbon dioxide and 80% nitrogen by volume. The correction factor was determined for an equivalence ratio of 0.7-1.2 and for initial pressures of 0.1 and 0.2 MPa.
A correlation purely based on simulations was determined by Fu et al. [12] using the CHEMKIN-PRO software and the Frassoldati et al. model [13], which contains 249 species and 7966 reactions. The laminar burning velocity simulations were carried out at pressures and temperatures of up to 0.5 MPa and 500 K, respectively, for a stoichiometric mixture. However, the running conditions at which the EGR correction factor was determined are unclear. The diluent composition is user defined and includes carbon dioxide, nitrogen and water.
The most recently suggested correction factor is based on a combination of simulations and experiments. Bhattacharya et al. [14] used a commercial gasoline (Shell V-Power as available in Germany) that, given the publication date of the paper, should be compliant with EURO VI regulations. For the experimental data, Bhattacharya et al. used a heat flux burner to determine the stretch free laminar burning velocities with the diluent comprising of carbon dioxide and nitrogen. Simulations were performed using the CHEMKIN-PRO package and a chemical mechanism constituting 877 reactions and 171 species [15]. This mechanism has been derived from a more complex mechanism containing 3796 reactions and 874 species [16,17]. The chemical model was validated against the burner experimental data, and a correlation for the EGR correction factor was proposed. Both experimental and numerical data was taken at a pressure of 0.1 MPa and unburned gas temperature of 432 K.
Another recent study on the effects of exhaust gas recirculation of a fuel for a advanced combustion engine (FACE-C) has found that the laminar burning velocity decreases quasilinearly as the levels of EGR increases [18] but does not specify a correction factor.
Currently, to the best of the authors' knowledge, no work exists that compares these models in terms of their accuracy under engine conditions. Middleton et al. [19] simulated the effect of EGR on the laminar burning velocity of isooctane air mixtures at high pressures (0.1-25 MPa), high temperatures (400-1000 K) and high EGR levels (0%-60% by mass) associated with engine studies. Although operating at engine-like conditions, the study did not directly compare different EGR correction factors. Therefore, it is not clear which models are useful to simulate the effect of EGR on combustion in SI engines or whether any of them are predictive at all. This is particularly important since recent works [20,21,22] have started to look at using very high levels of EGR in combustion engines. To enable robust and reliable computer modelling of the effect of EGR, this article provides a comparison of how well different correction factor models match measurements obtained from an engine. It also provides the first derivation of a correction factor directly from engine data. Using the Leeds University Spark Ignition Engine Data Analysis (LUSIEDA) software, we back calculate a correction factor for EGR levels of up to 25%. While this back-calculating approach comes with challenges regarding the proper choice of model parameters, by relying on engine data it captures the effects of the full composition of exhaust gas. The obtained correlation is then implemented into a quasi-dimensional predictive model and compared to the models mentioned above to see which one results in the best fit to the measured pressure trace. We show that while all models provide reasonably good agreement with experimental data at low levels of EGR (5%-10%), only the new found correction factor can match the measured pressure traces and mass fraction burned at higher levels of EGR (20%-25%) as existing models overestimate the reduction in laminar burning velocity due to the effects of EGR to varying degree.

Methodology
Both the forward model LUSIE (Leeds University Spark Ignition Engine) and the backward model LUSIEDA (Leeds University Spark Ignition Engine Data Analysis) are quasidimensional thermodynamic codes separating the combustion chamber into at least two zones. LUSIE is predictive, modelling pressure traces for specific engine conditions. In contrast, LUSIEDA is a reverse thermodynamic code that back calculates combustion parameters, such as the mass fraction burned, turbulent mass burning velocity and burned gas radius, from experimentally obtained in-cylinder engine pressure data.
In the first part of the article, we used LUSIEDA to back calculate burning velocities from pressure data to find the EGR correction factor correlation from experimental pressure traces in the engine for varying levels of EGR (see Experimental set-up section). Then, the new-found correlation as well as the four existing correction factor models in Table 1 were implemented into the predictive combustion code LUSIE to compute predicted pressure traces. Those were then compared against the original pressure data as well as a second set of data obtained under different running conditions to assess their suitability as a predictive model under engine conditions.

Experimental Set-Up
The engine from which experimental data was taken was a Single Cylinder Research Engine (SCRE) version of the latest Jaguar Land Rover gasoline Ingenium engine, housed at Imperial College London. A geometrical specification of the SCRE is shown in Table 2.
A full description of the engine, ancillaries, sub-systems and experimental methodology is available in the published literature [23,24].
As engine pressure trace data is used to derive the EGR correction factor correlation, the accuracy of the pressure data is important. A description of the pressure acquisition and associated errors is given by Smith et al. [23], stating that "High speed, crank angle resolved data was recorded using AVL Indicom v2.6 as part of an AVL Indiset Advanced Gigabit unit utilising a 14-bit analogue to digital convertor (maximum error of ±0.95 kPa, ±0.061 kPa and ±0.122 kPa for the in-cylinder, intake and exhaust pressure channels). A watercooled Kistler 6041B piezo-electric sensor (accurate to <1% of full scale), mounted flush with the combustion chamber surface, in combination with a Kistler 5064 charge amplifier were used to measure in-cylinder pressure. This 'dynamic' pressure was referenced to the intake manifold pressure (measured using a Kistler 4007 type sensor in conjunction with a Kistler 4665 signal conditioner) measured at the crank angle equidistant between the crank angles of maximum valve lift and intake valve closure." EGR was used throughout the duration of the currently reported study. By definition, the EGR rate was calculated as the ratio of the carbon dioxide measured within the intake manifold to the carbon dioxide content measured in the exhaust stream; both intake and exhaust CO 2 was measured using a Horiba MEXA One emissions analyser. EGR rate was calculated by the MEXA with the result being output to the data acquisition system. When required, EGR was sampled via a low-pressure loop system which resulted in EGR being introduced downstream of the intake air filter but upstream of a split in the intake system that was required to "switch" the single cylinder engine from running under naturally aspirated or boosted conditions. The hot EGR was cooled using a Ford DW12 watercooled EGR cooler with the rate being controlled via a Ford DW12 EGR valve.
Under naturally aspirated conditions, EGR was combined with the intake air before passing through a heat exchanger where the intake charge (air and EGR) temperature was controlled to 45°C, accurate to ±1°C, measured at the intake plenum. Under boosted conditions, intake air and EGR were combined at the same point in the intake system before they were passed through an Eaton V240 supercharger. The same heat exchanger described above was then used to "cool" the intake charge mixture to 45°C ± 1°C, measured at the intake plenum. To control the desired EGR rate, a closed-loop EGR rate control system was developed which allowed the user to define an EGR rate input, with the output defined as the measured EGR rate with the control functionality applied to the EGR valve.
Whilst the level of EGR introduced to the engine is typically very easy to quantify (using the methodology described previously), it is more difficult but not impossible to measure the quantity of residual gas left over in the chamber after combustion. In a paper by Peckham et al. [25], in-cylinder residual concentration was measured using in-cylinder Non-Dispersive Infra Red (NDIR) analysers that allowed crank-resolved measurements of CO and CO 2 and thus the calculation of the in-cylinder residual content. However, without this type of specialized equipment, in-cylinder residual concentration is typically quantified using 1D or CFD simulations.
Residual values typically vary from 5% to 15% with the residual concentration being influenced by valve overlap. Advantages of trapped residuals include a decrease in the generation of in-cylinder NO x emissions and a reduction in pumping losses, whereas disadvantages include an increase in the production of particulate matter and increased cycleto-cycle variation. Throughout the duration of the currently reported study, the crank angles of maximum opening for both intake and exhaust valves remained constant ensuring valve overlap and thus the level of trapped in-cylinder residuals remained broadly constant.
The effect of residuals in the forward model is accounted for by the model parameters for the 0% EGR case in the predictive engine simulations.

Backward Model
The EGR correction factor is defined as where f is the mass fraction of EGR, u l (f ) the laminar burning velocity with EGR and u l (f = 0) the laminar burning velocity without EGR. Therefore, to determine the correction factor from data, we need to determine u l for different levels of EGR and the no-EGR baseline. The in-cylinder residuals were calculated using the model suggested by Cho et al. [26] and included the back calculation to ensure that the unburned gas density was accurate. The hot residuals also effect the in-cylinder temperature at IVC for both the backwards and predictive model.
Several models exist in the literature that link the laminar burning velocity with the turbulent burning velocity using a correlation Eq. (2) with x being a model-dependent constant. In predictive simulations u t is computed from u l . For back-calculation, in contrast, we rearrange the model to find the laminar burning velocity from the turbulent burning velocity determined from pressure data. The full process is shown as a flowchart in Figure 1.
LUSIEDA calculates the pressure rise inside the combustion chamber incrementally during combustion: where P i and P i+1 are experimental pressure data at two consecutive crank angles, ΔP m is the change in pressure due to isentropic compression/expansion, ΔP ht is the change in pressure due to heat transfer, ΔP bb is the change in pressure due to blow-by and ΔP comb is the change in pressure due to combustion. The heat transfer and blow-by in this study (reverse and forward modelling) is calculated using the Woschni [27] and a "flow through orifice" model [28], respectively. An iterative method is used to determine the change in mass burned, Δm b , required for ΔP comb to equal the measured value found from Equation 3. The change in mass burned can then be used to calculate the turbulent mass burning velocity, where ρ u is the unburned gas density and A f is the burned gas area. The burned gas area is calculated by assuming a spherically propagating flame that is truncated by the cylinder. Currently, the backwards model assumes idealised pent chamber geometry, while the predictive model has been developed to accommodate for real engine geometry. LUSIEDA iterates to find a burned gas radius and uses this information to determine the area of the burned gas. Note that LUSIEDA uses a two-zone model approach which assumes that entrained fuel is burned instantaneously. This is because experimental pressure trace data does not allow flame front information to be computed, instead requiring additional optical data which is not readily available in an engine. The reverse analysis approach adopted here has been widely used in engine modelling [29,30,31,32].

Turbulent Burning Velocity Model
The original Zimont model describes a flame which transitions from a laminar flame to the thin reaction zone and wrinkled flame of the flamelet region and eventually to a region where the flame brush thickens [33]. The Zimont model reads: where u t0 is the fully developed turbulent burning velocity, A is a user defined constant, u ′ is the turbulent RMS velocity and Da is the Damköhler number. The Zimont model was updated to the Zimont-Lipatnikov model through the inclusion of a flame development factor and an extra laminar burning velocity term, u l , to correct for when the turbulence approached zero [34]. The Zimont-Lipatnikov model reads: where u k ¢ is the flame development factor multiplied by the turbulent RMS velocity, also known as the effective turbulent RMS velocity. In the back calculation of the correction factor, for the sake of simplicity, we use the mean u tr value over all experimental pressure traces instead of running 300 instances of the backward model and averaging those. The Damköhler number, Da, depends on the laminar burning velocity: Eq. (7) where τ t is the turbulence time scale, τ c is the chemical time scale, L is the integral length scale of turbulence and κ is the thermal diffusivity. The parameters L and u ′ for the studied engine were determined using CFD simulations [23] and were parsed into the code from an ASCII file.

Eq. (9)
In predictive simulations using the forward model, the constant A is fixed by tuning the model to measurements. This is not possible for the backward model but we can fix A to match more recently published models for the turbulent burning velocity. Recent research showed that the turbulent burning velocity depends on the stretch rate Markstein number, Ma sr , and the Karlovitz stretch factor [35,36]. Their model, which we refer to as the U/K correlation, reads .

Ma sr
Eq. (12b) Since the Markstein number is dependent upon the mixture composition, temperature and pressure, making direct use of their correlation in the backward model is challenging. No data exists in the literature about this dependence at the temperatures and pressures experienced in an Internal Combustion Engine (ICE). Furthermore, there is also no data regarding a similar fuel to that used in the present study or the effect that increasing levels of EGR diluent has on Ma sr . However, we can fix the constant A in the Zimont-Lipatnikov model so that it closely aligns with the U/K correlation for a given Markstein number. The Damköhler number is inversely proportional to the Karlovitz stretch factor, with an extra factor to account for the difference in turbulence length scales used: where λ is the Taylor microscale. If we ignore the added laminar burning velocity term that accounts for u ′ → 0, the Zimont-Lipatnikov model becomes:

Eq. (17)
A plot of the Zimont-Lipatnikov model with A = 0.33 is shown in Figure 2. The two models give very similar results and a value of A = 0.33 is therefore used for the backward model.
Note that the stretch rate Markstein number for isooctane was found to be around 6 at 0.1 MPa and 358 K [37] with an error of ±1. While our measurements were taken at a higher pressure and higher temperature, it has been shown that increases in pressure will lead to a decrease in Ma sr [37,38], whereas an increase in temperature leads to an increase in Ma sr [37]. Therefore, the value of Ma sr = 5.73 used in this study for gasoline seems reasonable but the lack of data in the literature prevents a more quantitative assessment. Experiments aiming to determine Markstein numbers at high temperatures and pressure would be an important area for future research.
It is worth stressing that the only free parameter in the backwards model is A and that the resulting correction factor is not particularly sensitive to it: changing A by ±10% was found to have a negligible effect on the back-calculated EGR correction factor, with a maximum change of 1.5% over all EGR values.

Forward Model
The predictive quasi-dimensional model LUSIE splits the combustion chamber into three zones with the extra zone, known as the entrainment zone, based on the work by Blizard and Keck [39]. Fresh gas is entrained into the flame at a rate: where m e is the mass of gas entrained into the flame brush, A fe is the flame surface area and u te is the turbulent entrainment burning velocity. The rate of mass burned is related to the mass entrained by: where τ b is the characteristic burn-up time: where C τ b is a constant. LUSIE has been introduced in detail and validated for numerous engines under an array of running conditions in previous publications [23,30,40,41,42,43].
The tunable constants within the forward model are the turbulent burning velocity constant, parameter A from the Zimont-Lipatnikov model, the characteristic burn-up time constant, C τ b , and the flame growth time. The turbulent burning velocity constant and characteristic burn-up time constant were tuned for the zero EGR case only.
The flame growth time is tuned to match the 0%-2% mass fraction burned for the simulations. The flame growth time needed to increase (slower early combustion) as the level of EGR increased to avoid the forward model over-predicting the early rate of combustion. For the engine used in the present study, we found that a linear relation between flame growth time and level of EGR worked well. This is analogous to the flame growth multiplier (FGM) factor used in the study by Robertson et al. [44] where it was also found that, as EGR increases, the time taken from 0% to 2% mass fraction burned increased. A more general study of how to model the impact of EGR on kernel delay would be another interesting direction for future research.

Results and Discussion
Laminar burning velocities were calculated using Equation 8 for a range of EGR values from the experimentally derived turbulent burning velocity values. To avoid spark effects and flame deceleration due to interaction with the cylinder walls, the burning velocity measurements were taken when the burned gas radius was 10-30 mm. The back-calculated laminar burning velocity plotted against in-cylinder pressure is shown in Figure 3. A logistic function was used to fit the data points and then to calculate the EGR correction factor for the varying levels of EGR shown in Figure 4.
The calculated EGR correction changes as the flame develops over time and pressure in the cylinder increases. This raises the question at what pressure in Figure 4 the correction factor should properly be chosen. The flame develops in three stages: initial acceleration, then propagation at an approximately constant speed and finally deceleration due to wall effects [29]. The correction factor should account for the impact of EGR on steady-state flame development, whereas the impact on early flame development is modelled by the time taken to form the flame kernel. Figure 5 shows the turbulent mass burning velocity plotted against the burned gas radius. While the steady-state phase with near-constant u tr is not very pronounced in the case studied here, it occurs around the peak of u tr which corresponds to a pressure of 3.0 MPa for the 0% EGR case. Therefore, we determine the proposed correction factor at that pressure value.
Obtained correction factors are plotted against EGR fraction in Figure 6, together with models from the literature. All models show an approximately linear decrease of the correction factor with EGR. To confirm that the correction factor we derive from data is not sensitive to the number of acquired cycles, we confirmed that when analysing only data from 100 or 200 cycles (in contrast to the full data from 300 cycles used here), the resulting correction factor is very similar to the one shown in Figure 6 (less than 1.5% average difference).
The models by Bhattacharaya et al. and Metghalchi and Keck align very closely. Both the model by Rodes and Keck and by Fu et al. give noticeably smaller correction factors for higher levels of EGR than the first three. The model derived in the present study finds that the EGR correction factor is larger at high levels of EGR compared to the literature. This could be attributed to the composition of the EGR containing species like hydrogen which burn much faster than those included in the simulated EGR within the literature. This is supported by a study by Mannaa et al. [18] who found that synthetic EGR reduced the laminar burning velocity more than real exhaust gas.

Eq. (21)
where f is the mass fraction of EGR, which was calculated experimentally as the ratio of CO 2 in the intake manifold to the ratio of CO 2 in the exhaust stream. The correlation found in the present study is similar to that suggested by

Eq. (23)
where f mole is the mole fraction of EGR has the same functional form but, due to different parameters, leads to significant differences in the obtained correction factors (compare for Figure 6). We experimented using a linear fit to our data instead of the quasi-linear fit but this led to worse agreement with the measured pressure traces, in particular for low levels of EGR. Eq. (25) where X i is the mole fraction of a single component in total diluents, n is the total number of diluents, μ 1 − μ 4 and ϕ m are correlation coefficients. For our case (ϕ = 1), the correlation is given by: where the mole fraction of each diluent was calculated using the isooctane balance equation This data was then divided into fast, middle and slow combustion cycles while removing data in-between these regions. This approach is commonly used when dealing with cyclic variability [40,42,45] and allows for a more detailed comparison of the simulated pressure curves. Fast, middle and slow cycle are determined by the peak pressure values, with the fast combustion cycles defined as P P max m ax ³ +s , where P max is the mean peak pressure value and σ is the standard deviation. The middle cycles are defined as P P P  Figure 7. The mean experimental pressure trace was calculated by finding the average pressure at each time step (0.1 crank angle degrees) for the 300 cycles. The simulation is in good agreement with experimental data and represents a middle combustion cycle.
The pressure traces and mass fraction burned profiles were then simulated for increasing levels of EGR, ranging from 5% to 25% (by mass) in increments of 5%. For the correlations derived using mole fractions of diluent (Rhodes and Keck and Fu et al.), the mass fraction user input was converted to the mole fraction within the LUSIE code, ensuring the correction factor used the appropriate value of EGR. Figures 8-12 show the simulated pressure traces and mass fraction burned profiles using the different models with  Figure 13 shows the root mean square error between simulated and experimental pressure values plotted against level of EGR for all models.
It can be seen that at the lower levels of EGR (5%-10%), all of the correlations provide a reasonable fit to the experimental data although the Rhodes and Keck model predicts a slower combustion cycle. This is somewhat unsurprising as all models give fairly similar correction factor values at 5% EGR as shown in Figure 6.
Substantial differences between models arise as the level of EGR increases to values between 15% and 25%. The very small correction factors in the Rhodes and Keck and Fu et al. models result in an increasingly poor match with experimental data, with errors increasing roughly linearly with EGR level.
The new correlation, the Bhattacharya et al. and Metghalchi and Keck model gives reasonable approximations throughout with similar errors up to around 15% EGR. Bhattacharya et al. and Metghalchi and Keck produce closely aligned pressure traces, tending slightly toward slower cycles for high levels of EGR. The new correlation from the present study improves on these two models for both the 20% and 25% EGR cases predicting pressure trace and mass fraction burned profiles much closer to the mean than any of the models from the literature. The difference at higher levels of EGR is possibly attributed to the diluent composition itself, with the current study recirculating actual exhaust gas as opposed to using a simulated EGR consisting of nitrogen (N 2 ) for the numerical model and a mixture of nitrogen and carbon dioxide for the heat flux burner and constant volume combustion bomb experiments.   The forward modelling was repeated for part load engine conditions to demonstrate that the new correlation is predictive at engine conditions that are different to those from which it was derived. The zero EGR case was tuned for the part load condition to match a mean cycle. Just as for the full load case, no further tuning was done for the simulations with EGR. The resulting RMS error for the pressure traces are presented in Figure 14. The results show again that the new correlation derived from engine conditions produces the lowest RMS errors across the range of EGR used, matching the results for the full load conditions shown in Figure 13.
From the four existing models, the Metghalchi

Conclusions
The article introduces a new method to determine correction factors for exhaust gas recirculation (EGR) from engine pressure trace data by using a reverse thermodynamic model. Using this approach, a new correlation for the effects of exhaust gas diluent on the laminar burning velocity of a EURO VI specification gasoline is derived and compared against existing models by Metghalchi and Keck [10], Rhodes and Keck [11], Fu et al. [12] and Bhattacharya et al. [14]. It is found that existing models tend to overestimate the impact of EGR on laminar burning velocity, probably because their derivation did not consider all chemical species contained in exhaust gas. This is in agreement with the study by Mannaa et al. [18] who found that synthetic EGR led to a greater reduction in laminar burning velocity when compared to real exhaust gas.
The new correlation and models from the literature were then implemented into a predictive combustion code to compare their predictive modelling capabilities. For the new correlation, simulated pressure traces and mass fraction burned profiles show good agreement with experimental data over the full range of 5%-25% EGR (by mass.) under full load and throttled engine conditions. In contrast, the models by Rhodes and Keck and Fu et al. show reasonable agreement up to 5% EGR but suffer from rapidly rising errors for higher EGR levels. The models by Metghalchi and Keck and Bhattacharaya et al. agree reasonably well even at higher levels of EGR but are outperformed by the new correlation. In line with the too high correction factors, we find that, in the forward model, the four existing correction factors reduce the burning velocity too much, thus tending to model slower combustion cycles instead of mean cycles as EGR increases. For the Rhodes  Definitions/Abbreviations aTDCgx -after top dead centre gas exchange bTDCgx -before top dead centre gas exchange