Multidimensional Numerical Simulations of Knocking Combustion in a Cooperative Fuel Research Engine Free

Transportation, alone, accounts for around 57% of the total oil produced annually [1]. Reducing emissions from transportation is crucial to handling the challenges of air pollution, climate change, and increasing demand for mobility. Improving vehicle powertrain efficiency is the most feasible way to meet ever-more stringent CO₂ emission regulations and increasing demand for transportation fuels. Over the years, the U.S. light-duty market has been dominated by gasoline spark-ignition (SI) engines [2]. Reducing engine displacement (downsizing) while maintaining the torque output by operating at high loads (boosting) is considered an effective strategy to improve fuel economy of SI engines as a result of increased thermal efficiency and lower CO₂ emissions. However, high loads also bring about more severe in-cylinder thermodynamic conditions, promoting the likelihood of abnormal combustion phenomena such as knock.

In a SI engine, knock occurs as a result of premature auto-ignition of a portion of the end-gas (unburnt) mixture, before it is consumed by the turbulent premixed flame generated by the spark plug [3]. This can cause very rapid heat release thereby leading to sharp pressure rise and high-speed pressure waves in the combustion chamber. These pressure fluctuations manifest themselves as a metallic clanking noise and can potentially lead to the damage of critical engine parts such as liner and piston. Knocking combustion imposes a stringent constraint on the performance and efficiency of downsized, boosted SI engines since it prohibits the use of more advanced spark timings, higher compression ratios (CRs), and higher boost pressures [4].

The occurrence of knock depends on a number of fuel properties (flame speed, ignition delay, heat of vaporization, lower heating value, etc.) as well as on the prevailing engine operating conditions. One of the most important fuel parameters relevant to SI engines is the antiknock quality. Octane number (ON) is the most widely used measure by which a fuel's resistance to auto-ignition (propensity to knock) is characterized [5]. A higher ON indicates a higher knock resistance. ONs are measured under two different engine conditions in a highly calibrated cooperative fuel research (CFR) engine, and are defined as research octane number (RON) and motor octane number (MON). Primary reference fuels (PRFs) provide the scale for both ONs, with iso-octane assigned RON = MON = 100 and n-heptane defining RON = MON = 0. A test fuel's RON (or MON) is determined by inducing “standard knock intensity” using the test fuel by adjusting some operating parameters of the CFR engine, and then finding a PRF mixture that matches the test fuel's knocking behavior. The ON of the test fuel is defined as the volume percentage of iso-octane in that particular PRF blend [6,7].

The RON and MON tests interrogate different reaction conditions (discussed in the Cooperative Fuel Research Engine Specifications and Operating Conditions section). In general, RON of a practical gasoline fuel is greater than its MON, and the difference between the two is denoted as the fuel octane sensitivity (S = RON – MON). On the other hand, by convention, for paraffinic PRF blends, S = 0. S > 0 for real fuels is mainly attributed to significantly less pronounced negative temperature coefficient behavior exhibited by them as compared to PRFs, which makes them less resistant to auto-ignition under MON conditions [2,8,9]. In addition, recent studies have indicated that heat of vaporization of a fuel can contribute to its octane sensitivity as well [10].

Kalghatgi [11,12] introduced octane index, OI = RON – K*S, as a more appropriate metric to quantify the antiknock quality of fuels by accounting for in-cylinder conditions, which vary with engine design and operating conditions. Here, K is an empirical constant that is representative of the pressure-temperature history of the unburnt mixture in the combustion chamber. By definition, K = 0 at the RON condition and K = 1 at the MON condition. Kalghatgi [13–15] has shown that modern downsized boosted SI engines operate at negative K values (i.e., “beyond RON” conditions), so that a fuel with higher S would be more resistant to knock and enable higher efficiency. This has motivated further studies to investigate potential octane boosters that can be blended with gasoline to increase both RON and S [16,17].

A better understanding of the interaction between fuel properties and engine operating conditions is necessary to develop knock mitigation strategies and extend the high load limit of advanced SI engines. In this context, multidimensional modeling and computational fluid dynamics (CFD) can complement engine experiments and serve as a valuable tool to provide more insights into the knocking phenomena. Successful prediction of knock requires accurate modeling of turbulent flame propagation as well as end-gas auto-ignition. One of the commonly used knock modeling frameworks employs the G-equation (also known as “level set”) combustion model [18,19] to describe the propagation of the premixed turbulent flame front, and a well-stirred-reactor model [20] with detailed/reduced chemical kinetic mechanisms to capture auto-ignition in the end-gas. Tan and Reitz [19] implemented the G-equation model into the KIVA-3V code. This model was further extended for knock prediction by including chemical kinetics calculation ahead of the flame front [21–23]. More recently, Ali et al. [24] used this hybrid approach to investigate the effect of pre-ignition on the likelihood of super knock in a direct-injection SI engine. In these previous works [19,21–24], the fuel laminar flame speed, an input to the G-equation model, was calculated using empirical correlations [25,26], which are available and valid for simple fuels only, such as iso-octane, thereby precluding the study of more complex multicomponent fuels or fuel blends with this type of modeling approach.

In the present work, a numerical approach was implemented to predict knock in a CFR engine. As mentioned earlier, the CFR engine is practically relevant to studying fuel effects on auto-ignition propensity at typical RON/MON conditions. The model can also be effectively used to probe “beyond RON” conditions relevant to advanced SI engines. In addition, the virtual CFR engine model allows for a detailed investigation of the physical and chemical mechanisms governing knock in a simple engine geometry. For combustion modeling, the hybrid approach combining G-equation and multizone [27] models was employed. Additionally, a novel methodology of accounting for fuel effects on flame propagation based on tabulated laminar flame speed was incorporated, thereby ensuring that the model could be applied to any fuel of interest. The numerical model was used to predict the critical compression ratio (i.e., CR at knock onset) under RON/MON conditions, for different PRF blends. A new skeletal PRF mechanism was developed and used to both capture end-gas auto-ignition and generate the laminar flame speed tables.

The remainder of the paper has the following structure: The Cooperative Fuel Research Engine Specifications and Operating Conditions section describes the CFR engine specifications and the RON and MON test conditions studied, along with available experimental data. Details of the numerical setup, combustion model, and skeletal mechanism are provided in the Modeling Approach section. Subsequently, results from the numerical simulations are discussed. Finally, the paper concludes with a summary of the main findings.

Cooperative fuel research engine is a carbureted single-cylinder engine, which is used to measure a fuel's RON and MON based on standardized protocols [6,7]. The CFR engine specifications and operating conditions for RON/MON tests are detailed in Tables 1 and 2, respectively. RON and MON tests represent two different pressure-temperature trajectories. The MON trajectory has a higher temperature at a given pressure as compared to the RON trajectory [28,29].

In the present work, RON/MON test conditions for two PRF blends, PRF60 and PRF80, were considered for the numerical study. Simulations were performed at different compression ratios for each PRF blend and operating condition. The objective was to accurately capture the transition from no knock to knock with increasing compression ratio and compare the predicted “CR at knock onset” against the experimental data available [28]. The CR sweeps simulated in the present work along with the experimental curves of the critical compression ratio are shown in Fig. 1. A threshold peak pressure fluctuation of 0.5 bar was used as the criterion for knock.

A commercial 3D CFD code, CONVERGE (version 2.3) [30], was used to perform numerical simulations for the closed part of the cycle, from IVC to EVO. Computational mesh for the CFR engine geometry is shown in Fig. 2. The compression ratio was varied by changing the position of the piston at bottom-dead-center while keeping the stroke constant. A uniform mixture and temperature distribution was specified at IVC. The temperature at IVC was specified to be 380 K and 422 K, for RON and MON conditions, respectively [28]. On the other hand, a constant pressure of 1 bar at IVC was prescribed for all cases. An air–fuel ratio of unity was chosen as the maximum knock intensity is typically achieved close to stoichiometric mixture conditions [29]. The piston and cylinder head temperatures were assigned a temperature of 440 K [31–33], whereas 410 K was used as the cylinder wall temperature [34–37], based on previous experimental data reported in the literature.

CONVERGE uses a modified, cut-cell Cartesian method for grid generation directly during runtime. In addition, it has the capability to include fixed embedding of cells, i.e., increasing the grid resolution with respect to the base grid size a priori and adaptive mesh refinement to refine areas where the subgrid field is the largest [30]. In this study, a base mesh size of 1 mm was used. One level of fixed boundary embedding was specified near the cylinder head, piston, and wall (cell size of 0.5 mm), while three levels of fixed spherical embedding (radius of 0.5 mm) were employed to resolve the spark (cell size of 0.125 mm). The spark timing for RON and MON conditions was set to 13 deg BTDC and 26 deg BTDC, respectively. Two levels of adaptive mesh refinement were included based on the velocity and temperature subgrid scales of 1 m/s and 2.5 K, respectively. This was essential for resolving the propagating turbulent flame front and resulted in the minimum grid size of 0.25 mm. The maximum cell count per simulation was 5× 10⁶ approximately. It has been shown in previous Reynolds-averaged Navier–Stokes (RANS) [38–40] and large eddy simulation studies [41,42] that a minimum grid size of ∼0.25–0.5 mm is sufficient for simulating normal SI combustion and knock. In-cylinder turbulence was modeled using the RANS based renormalized group k-ε model [43] with wall functions. The model proposed by Han and Reitz [44] was employed to account for wall heat transfer. The combustion model and PRF skeletal mechanism used for the simulations are described in the Combustion Model and Skeletal Mechanism for Primary Reference Fuel sections.

In order to capture the high frequency local pressure oscillations induced by knocking, the maximum value of Mach Courant–Friedrichs–Lewy (Mach CFL) number was set to 1.0 after spark timing. Moreover, several monitor points were set up at multiple locations near the edge of the combustion chamber to obtain pressure fluctuation information, which is analogous to measured signals using a pressure transducer in experiments. The simulated local pressures at these locations were used as the primary parameter to quantify knock intensity (KI). In addition, the in-cylinder temperature and formaldehyde (CH₂O) [23] distributions were used to visualize knock.

where u′ is turbulent velocity fluctuation, s_L is the laminar flame speed, Da is the Damköhler number, and a₄, b₁, and b₃ are standard modeling constants. A more detailed description of the G-equation model can be found in Refs. [18] and [30].

The laminar flame speed in Eq. (3) is generally calculated using empirical correlations [24–26,40]. However, a major limitation of these correlations is that they are valid for only simple fuels, such as iso-octane, ethanol, methanol and their blends. Therefore, in this work, a more general methodology to incorporate fuel effects was developed based on tabulated flame speeds. In particular, a lookup table for s_L as a function of temperature, pressure, and equivalence ratio was generated for each PRF blend investigated, based on a skeletal kinetic mechanism (discussed in the Skeletal Mechanism for Primary Reference Fuel section), using the CONVERGE 1D laminar flame speed solver [30]. Interpolation schemes were implemented in a user-defined function to provide appropriate values of s_L to Eq. (3) from these lookup tables. The ranges of T/P/Ф used for generating the tables are listed in Table 3. It must be noted that the dimensions of the flame speed tables can be readily extended to include exhaust gas recirculation as well, if needed.

On the other hand, both auto-ignition in the end-gas and postflame chemistry in the burnt region were captured using the SAGE detailed chemistry solver [46] along with the multizone approach, with bins of 5 K in temperature and 0.05 in equivalence ratio [27]. Although it does not utilize an explicit turbulent combustion closure [47–49], recent work has shown that it works well if a fine mesh is employed, as this significantly reduces subgrid effects [39,50–53].

In our previous work [54], a skeletal mechanism with 165 species and 839 reactions was developed based on a 1389-species detailed mechanism for gasoline surrogate mixtures with octane numbers in the range 60–100 [55], employing directed relation graph (DRG) [56–58], DRG-aided sensitivity analysis (DRGASA) [59], and isomer lumping [60] reduction techniques. This mechanism was valid for toluene primary reference fuels and toluene primary reference fuel-ethanol blends.

In the present work, DRGASA was reperformed to remove toluene and ethanol related species and reactions to further reduce the size of the skeletal mechanism since only PRF mixtures were studied here. The target parameters in DRGASA included ignition delay time of constant-pressure autoignition and the extinction residence time in perfectly stirred reactors for PRF mixtures with ON between 60 and 100. The reduction was performed within the same parameter range as earlier [54]: pressure from 1 to 100 atm, equivalence ratio from 0.5 to 1.3, and initial temperature of 750–1800 K for auto-ignition, thereby incorporating the low temperature chemistry. As a result, a PRF skeletal mechanism with 121 species and 538 reactions was obtained.

The new PRF skeletal mechanism was validated for both homogeneous ignition delay and laminar flame speed against experimental data reported in the literature. A comparison of ignition delays predicted by the detailed and skeletal mechanisms with shock tube experimental data [61] is shown in Figs. 3–5, for PRF60/air, PRF80/air, and PRF90/air mixtures, respectively. It can be observed that the skeletal mechanism predicts experimental data very well and is in good agreement with the detailed mechanism.

Figures 6 and 7 present experimental laminar flame speed data for iso-octane/air [62–67] and PRF90/air [64] mixtures, respectively, for different equivalence ratios at initial pressure of 1 atm and initial temperature of 298 K, along with the corresponding computed results. Overall, the numerical flame speeds are also in close agreement with the experimental values. These validation studies indicate that the skeletal mechanism is suitable for capturing end-gas auto-ignition as well as flame propagation characteristics for PRFs with reasonably good fidelity.

Numerical simulations were performed for the compression ratio sweeps denoted in Fig. 1, using the numerical setup described earlier. The local pressure traces and in-cylinder visualizations were used to characterize the combustion behavior in each case. The capability of the virtual CFR engine model was assessed by investigating whether the compression ratio at knock onset was predicted accurately for the chosen PRF blends.

In order to differentiate between knocking and nonknocking cases, the criterion of a threshold KI was employed. Based on previous studies [68,69], a value of 0.5 bar was used in the present work so that KI > 0.5 bar corresponded to a knocking condition. KI was quantified using the maximum amplitude of pressure oscillations (MAPO) [70]. To determine MAPO, the local pressure evolution at the monitor points was first analyzed and the one with the highest peak pressure fluctuation was chosen. It was observed consistently for all knocking cases that the peak pressure fluctuation was highest for the monitor point placed near the edge of the combustion chamber diametrically opposite to the spark plug. This was mainly attributed to that location being farthest from the propagating flame front thereby allowing for most severe end-gas auto-ignition. Therefore, the local pressure trace at this point was considered for knock characterization. Fast Fourier transform of the local pressure trace was calculated and higher frequencies in the range 4–20 kHz [71,72] (relevant to knock) were filtered out using a band-pass filter. Subsequently, the knocking pressure was obtained by computing the inverse fast Fourier transform of the filtered pressure signal and its maximum absolute value was defined as the MAPO.

Figure 8 depicts the temporal evolution of local in-cylinder pressure for PRF60/air mixture over the chosen CR sweep under RON operating condition. It can be seen that a clear transition from no knock to knock occurred with an increase in CR. The results of the corresponding MAPO analysis in Fig. 9 show that KI = 0.16 bar for CR = 5.5 and KI = 1.38 bar for CR = 5.6. This indicates that the critical CR would lie somewhere in between these two CRs. On the other hand, the experimental value of critical compression ratio is 5.55, which also falls in the same CR interval. This indicates that the numerical results are in good agreement with experiments within the uncertainty of 0.1 CR units.

Figure 10 shows the in-cylinder temperature and CH₂O distributions on a vertical cut plane (passing through the spark plug) for CR = 5.65, at different instants. Figures 10(a)–10(c) clearly show the buildup of CH₂O radical in the end gas, attributed to low temperature chemistry. This increases the reactivity of the end-gas mixture. Subsequently, the CH₂O radicals are consumed rapidly and the end-gas auto-ignites ahead of the flame front (Figs. 10(d)–10(f)), leading to the knocking behavior.

The local pressure traces corresponding to the parametric cases corresponding to PRF80/air mixture at RON condition are presented in Fig. 11. In this case, the critical CR is bracketed by CR = 6.1 (KI = 0.44 bar) and CR = 6.15 (KI = 0.81 bar) based on the MAPO analysis (as shown in Fig. 12), whereas the corresponding experimental value is 6.11. Hence, the transition to knock is again predicted very well, within the uncertainty of < 0.05 CR units. Although not shown here (for the sake of brevity), similar trends were observed for the two PRF blends under MON conditions as well, with the numerical values of critical CR corresponding well to the experiments.

The present numerical study demonstrates that the virtual CFR engine model is able to capture the sensitivity of knock propensity to fuel chemistry, for different fuels (PRFs in this case) and under different operating conditions (RON and MON). It must be noted that the initial/boundary conditions for the numerical simulations were based on previous literature and only closed-cycle simulations were performed. Therefore, as part of ongoing work, a more detailed validation of the computational model is being carried out against CFR engine experiments [73] using multicycle simulations. In addition, the numerical model will also be employed in future studies, to assess the capability of turbulent auto-ignition regime theory [74–80] in predicting super-knock and pre-ignition phenomena.

In this work, a numerical approach based on multidimensional CFD was developed to capture knocking combustion in a CFR engine. For combustion modeling, a combination of two different approaches was employed. The turbulent flame propagation was modeled using the G-equation approach, whereas a multizone model was used to capture end-gas auto-ignition. A novel framework based on tabulated laminar flame speed was developed to account for fuel chemistry effects on flame propagation more accurately as compared to empirical correlations. Numerical experiments were performed for two PRF blends, PRF60 and PRF80, at different CRs under RON and MON conditions. A compact 121-species PRF skeletal mechanism was developed from a more detailed gasoline surrogate mechanism using the DRGASA reduction technique and was extensively validated against literature data on homogeneous ignition delay and laminar flame speed. Subsequently, the skeletal mechanism was used to model end-gas auto-ignition and generate laminar flame speed tables for the PRF blends. The CFD simulations were able to predict the transition from no knock to knock as the compression ratio was increased. The critical compression ratio (compression ratio at knock onset) predicted by the numerical simulations was in very good agreement with available experimental data for both PRF blends. This study, therefore, demonstrated the capability of the computational model to capture the sensitivity of knock propensity to fuel chemistry for the RON and MON conditions explored.

The authors would like to thank Dr. Benjamin Wolk of Sandia National Laboratories and Professor J. Y. Chen of University of California, Berkeley for providing the surface mesh of the CFR engine geometry. The authors are also thankful to Dr. Marco Mehl of Lawrence Livermore National Laboratory (LLNL) for providing the experimental data relevant to the RON and MON tests. In addition, useful discussions with Dr. Chris Kolodziej and Dr. Thomas Wallner of Argonne National Laboratory are gratefully acknowledged. The authors would like to acknowledge the HPC center at the National Renewable Energy Laboratory for computing time on the Peregrine cluster as well as the Laboratory Computing Resource Center at Argonne National Laboratory for computing time on the Blues cluster that were used in this research.

The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (Argonne). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DEAC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. This research was partially funded by DOE's Office of Vehicle Technologies and Office of Energy Efficiency and Renewable Energy under Contract No. DE-AC02-06CH11357. The authors wish to thank Gurpreet Singh, Kevin Stork, and Leo Breton, program managers at DOE, for their support. This research was conducted as part of the Co-Optimization of Fuels & Engines (Co-Optima) project sponsored by the U.S. Department of Energy (DOE) Office of Energy Efficiency and Renewable Energy (EERE), Bioenergy Technologies and Vehicle Technologies Offices.