Modelling Study of Cycle-To-Cycle Variations (CCV) in Spark Ignition (SI)-Controlled Auto-Ignition (CAI) Hybrid Combustion Engine by Using Reynolds-Averaged Navier–Stokes (RANS) and Large Eddy Simulation (LES)

: The spark ignition (SI)-controlled auto-ignition (CAI) hybrid combustion is characterized by early ﬂame propagation combustion and subsequent auto-ignition combustion. The application of combined SI–CAI hybrid combustion can be used to effectively extend the operating range of CAI combustion and achieve smooth transitions between SI and CAI combustion modes. However, SI–CAI hybrid combustion can produce signiﬁcant cycle-to-cycle variations (CCV). In order to better understand the sources of CCV and minimize its occurrence, the large eddy simulation (LES) and Reynolds-averaged Navier–Stokes (RANS) approaches were employed in this study to model and understand the cyclic phenomenon of SI–CAI hybrid combustion. Both the multi-cycle LES and RANS simulations were analyzed against the experimental measurements in a single cylinder engine at 1500 rpm and a 5.43 bar average indicated the mean effective pressure (IMEP). The detailed analysis of the in-cylinder pressure traces, IMEP, in-cylinder peak pressure (PP), peak pressure rise rate (PPRR) and the crank angles with fuel mass burned fraction at 10%, 50%, 90% and mode transition was performed. The results indicate that overall, the adopted LES simulations could effectively predict the cyclic variations in the hybrid combustion observed in the experiments, while the RANS simulations failed to reproduce the cyclic characteristics at the chosen engine operating conditions. Based on the LES results, the correlation and visualization studies indicate that the cyclic variations in the local velocity around the spark plug lead to the variations in the early ﬂame propagation, which in turn produce temperature ﬂuctuations among the cycles and result in greater variations in the subsequent auto-ignition combustion events.


Introduction
The spark ignition (SI) can be introduced to a gasoline engine that operates with compression ignition combustion operation, which is referred to as controlled auto-ignition (CAI), homogeneous charge compression ignition (HCCI) or gasoline compression ignition (GCI), to assist the control of auto-ignition [1,2]. This SI-CAI hybrid combustion, or spark assisted compression ignition (SACI), has higher thermal efficiency and lower NOx emissions compared to conventional SI combustion. Compared to CAI combustion, the SI-CAI hybrid combustion has a lower heat release rate and broader load range. In addition, this combustion concept enables smooth transitions between SI and CAI mode [3][4][5]. Therefore, in recent years, significant work has been carried out on SI-CAI hybrid combustion. Zhou et al. [6] applied Miller cycle strategies to adjust the ratio of fuel consumed by early spark ignition flame propagation and subsequent auto-ignition combustion. They found that the Miller cycle is an effective approach to optimize the SACI combustion process while high CCV of SI-CAI hybrid combustion. The single-cycle CFD simulations used in such studies were able to identify some key factors that would affect the steady state hybrid combustion, but not the cyclic variation involved. In order to identify the origins of CCV of the hybrid combustion, a preliminary modelling study using large eddy simulations (LES) was performed for 15 cycles and it was found the LES could reflect the cyclic variations in SI-CAI hybrid combustion with limited simulation cycles [28]. In this study, the multi-cycle RANS (12 cycles) and LES (32 cycles) simulations were performed to model the SI-CAI hybrid combustion with large CCV. The detailed analysis and comparisons were carried out by the LES and RANS simulations on the in-cylinder pressure traces, IMEP, peak pressure (PP), peak pressure rise rate (PPRR) and the crank angles with fuel mass fraction burned at 10%, 50%, 90% and mode transition (CA10, CA50, CA90 and CAT). Based on the LES results, the correlation study was performed to unveil the key factors leading to the strong cyclic combustion characteristics of the SI-CAI hybrid combustion process. To the authors' knowledge, it is the first time that such a study has been carried out to investigate the CCV of SI-CAI hybrid combustion with combined multi-cycle RANS and LES engine simulations.

Engine Experiments
The experiment data were obtained on a single cylinder gasoline engine with dual independent variable valve lifts and timing devices on the intake and exhaust camshafts [21,28]. The engine specification is given in Table 1. In order to achieve SI-CAI hybrid combustion with an acceptable peak pressure rise rate, both external and internal EGR were applied in this study. The negative valve overlapping (NVO) was used to enable internal EGR. A piezo electric transducer (Kistler 6125A) via a charge amplifier (Kistler 5018) was used to measure the in-cylinder pressure. An engine operating point at 1500 rpm and 5.43 bar IMEP (in average of 100 cycles) was adopted to study the cycle-to-cycle variations in SI-CAI hybrid combustion, as this operating point is within the transitioning zone between pure CAI and SI modes with a relatively high coefficient of variation (COV) of IMEP at 11.81%. The valve parameters, including the intake/exhaust valve opening timings (IVO/EVO), closing timings (IVC/EVC) and lifts (IL/EL), and other operating conditions are listed in Table 2. The residual gas fraction (RGF) after the intake valve closing (IVC) was calculated to be 25%. The crank angles used in this paper refer to the top dead center (TDC) of the combustion stroke.

RANS Modelling
The simulations were performed in commercial STAR-CD software. In RANS simulations, the flow and turbulence were modeled with RNG k-ε model [29]. The turbulent mixing effects and chemical kinetics in the SI-CAI hybrid combustion were predicted with a set of premixed flame propagation and auto-ignition models. The three-zones extended coherent flame model (ECFM3Z) [30] was used as the framework of the combustion model as it considers both premixed and diffusion flame propagation as well as auto-ignition combustion.
As detailed in [30], there are three zones (a pure fuel zone, a pure air and residual gas zone and a mixed zone) in ECFM3Z. The mixed zone, where the combustion happens, is determined by the turbulent and molecular mixing between gases in the other two zones. The flame surface density equation is adopted in order to describe the flame propagation process and predict its reaction rate. In order to describe the reaction intensity of flame propagation combustion, average flame surface density is used in this study and defined as the local area of flame surface per unit of volume (m −1 ). The auto-ignition of the mixture is predicted with tabulated chemistry approach [31]. The tabulated database of the auto-ignition delay time was constructed by performing chemical kinetic calculations under different thermodynamic and dilution conditions (pressure: 1-60 bar, temperature: 480-1520 K, equivalence ratio: 0.2-3 and residual gas fraction: 0-90%) with a reduced gasoline surrogate mechanism [32]. The coupled CFD simulations with tabulated database can significantly reduce the central processing unit (CPU) time requirement compared with the CFD simulations coupled with the chemical kinetic mechanism. By using the tabulated chemistry approach, the auto-ignition tendency is defined to describe the close degree of mixture from auto-ignition in each cell, and it ranges from 0 (no tendency to auto-ignition) to 1 (occurrence of auto-ignition). The reaction rate of the auto-ignition combustion can then be determined by the combustion characteristic time once auto-ignition occurs.
The calculation of average flame surface density and the auto-ignition tendency determine the reaction regime in each cell. When the average flame surface density is greater than 0 in a cell, the available fuel/air mixture will be consumed by the flame propagation combustion by following the flame surface density equation. By contrast, when the autoignition tendency in a cell achieves 1, the available fuel/air mixture will be consumed by auto-ignition combustion by following the tabulated chemistry approach. The reaction rates of the flame propagation and auto-ignition also determine the species concentrations during the combustion. The application of these models enables the prediction of the SI-CAI hybrid combustion. The details of the RANS modelling of SI-CAI hybrid combustion can be found in a previous work [28].
In this work, the spark ignition model was upgraded with the consideration of the development of spark kernel [33]. The time delay between the spark and the appearance of the flame surface is modelled by using an indicator function I that starts from 0 and eventually reaches 1.
where ρ * is the ratio between the current gas density and the air density at 1 bar and 300 K. τ f is given by τ f = δ l /U l , where δ l is laminar flame thickness and U l is laminar flame speed. A1 igni , A2 igni are the tuning parameters with default values of 1.0 and 3.0, respectively. A spherical flame surface with radius R k is formed once I reaches 1.
where T b is the burnt gas temperature and T u unburnt gas temperature. F actker and R klimit are the tuning parameters. Then, the flame surface Σ init is distributed in the computational domain with a Gaussian distribution with a mean at X µ = R k and a standard deviation of X σ = 1 × 10 −4 + u'(t − t spark ). The distribution is also constrained by the following equation: where V is the total volume.

LES Modelling
The sub-grid k model [33,34] was used to model the flows and turbulence in the LES simulations. The model has been successfully applied in LES modelling of flow [35], spray [36] and combustion [37] in combustion engines. As recommended by Speziale [34], a transport equation is solved for the sub-grid scale (SGS) kinetic energy to derive an appropriate turbulence velocity scale, which is as follows: where the model constants C 1 and C 2 are derived from turbulence theory and are set equal to 1 and 0.05, respectively. τ SGS,ij is sub-grid scale stress, s ij the rate of strain tensor, k SGS the SGS turbulent kinetic energy. ∆ is the grid filter size and is defined as V cell 1/3 . The LES version of the three-zones extended coherent flame model (ECFM3Z) [37,38] was applied for the LES modelling of SI-CAI hybrid combustion, with the same modelling concept adopted for the RANS simulations of the SI-CAI hybrid, combustion as detailed in previous section. In the LES version of the ECFM-3Z model, a new equation for modelling the SGS flame surface density Σ (Σ = ρσ) is considered as the following [38]: where σ c is a constant with a typical value of 1, D is the molecular diffusivity, v sgs is the SGS viscosity, S C t is the turbulent Prandtl number, S d is the propagation speed relative to the mean flow field of the iso-c surface, u is the Favre filtered velocity vector, n is normal to the iso-surface of the filtered progress variable c, S l is the theoretical laminar flame speed, δ l is the laminar flame thickness; Γ is the ITNFS (net flame stretch) function [39]. β and c * are model constants equal to 4/3 and 0.5, respectively. Σ lam is the laminar part of flame surface density.
As mentioned earlier, the equations for the filtered species, tracers and unmixed species remain unchanged with respect to their RANS counterparts, except that the turbulence timescale k/ε is now the SGS turbulence timescale and not the integral (RANS) value.
In addition, a standard LES spark ignition model [33,38] was applied. The first step is to calculate the flame kernel radius via the same standard ignition model used in RANS formulations. This is usually of mm order of magnitude and can be used in LES as a first approximation. The second step is the calculation of a suitable profile for the burnt gases. The ignition profile is given by where c 0 is a constant. F is a model parameter that can be used to shape the ignition profile. Higher values of F will result in stronger ignition. ∆ is the average combustion filter length. The distribution of the flame surface density is calculated as Therefore, the following condition is satisfied.
When any point in the domain reaches c = 1, the control is transferred to the transport equation and normal flame propagation begins. The S t is the total flame area wrinkled by the turbulence [38,40]. Table 3 shows initial and boundary conditions in the simulations. The inconsistency of the wall temperature was considered, and the temperatures of the cylinder liner and the piston were assessed with the cylinder head temperature and coolant temperature [19]. The intake mixture at the inlet boundary was set as the homogeneous fuel/air mixture at stoichiometric equivalence ratio and the exhaust gas to ensure the same external EGR rate with experiments. The CFD simulations start from 400 • CA after TDC (just before the intake valve opening (IVO) timing) and continue for the 32 consecutive cycles for LES simulations and 12 cycles for RANS simulations. In order to exclude the impact of the initial conditions, the first two cycles of LES and RANS simulations were excluded for analysis. Therefore, only 30 cycles LES and 10 cycles RANS simulations were analyzed in detail.

Engine Mesh and Numerical Methods
The engine moving mesh, as shown in Figure 1, was generated by ES-ICE software. There are around 870,000 grids with average grid size of 0.96 mm at the bottom dead center and 550,000 grids with average grid size of 0.79 mm at TDC, respectively. The mesh around the spark plug was refined with a grid size of 0.6 mm on average to accurately predict spark kernel formation and early flame propagation process.
The pressure-implicit with splitting of operators (PISO) algorithm was adopted to solve the equations. The monotone advection and reconstruction scheme (MARS) was used to discretize equations of momentum, turbulence kinetic energy and turbulence dissipation. The temperature and density equations were discretized with upwind differencing scheme (UD) and central differencing scheme (CD), respectively. The angular time-step in the simulations was fixed at 0.1 • CA for RANS simulations and 0.05 • CA for LES simulations. center and 550,000 grids with average grid size of 0.79 mm at TDC, respectively. The mesh around the spark plug was refined with a grid size of 0.6 mm on average to accurately predict spark kernel formation and early flame propagation process. The pressure-implicit with splitting of operators (PISO) algorithm was adopted to solve the equations. The monotone advection and reconstruction scheme (MARS) was used to discretize equations of momentum, turbulence kinetic energy and turbulence dissipation. The temperature and density equations were discretized with upwind differencing scheme (UD) and central differencing scheme (CD), respectively. The angular timestep in the simulations was fixed at 0.1 °CA for RANS simulations and 0.05 °CA for LES simulations.
The turbulence resolution parameter (M) proposed by Pope [41] was evaluated at different crank angles to assess the applicability of the LES results. M is defined as the ratio of the sub-grid scale (SGS) kinetic energy to the total kinetic energy. M of 0 denotes a direct numerical simulation (DNS) simulation where all the turbulence scales are resolved, while the M of 1 denotes a RANS simulation where no turbulence scales are resolved and all are modelled. It is recommended that M should be less than 0.2 for an adaptive LES simulation [41]. In this study, M is significantly lower than 0.2 in the LES simulations [28]. The spark plug region shows a relatively higher turbulence resolution parameter M, in particular at the near wall region of the ground electrode with M at 0.14, due to complex geometric features. The majority of the region around the spark plug gap has a turbulence resolution parameter well below 0.08. Figure 2 shows the in-cylinder pressure profiles from the experiments, RANS and LES simulations. The grey envelop highlights the large variations in the experimental results (100 cycles). The dark grey dash curve shows the averaged pressure of the experimental results. As shown in Figure 2a, the RANS simulation results show little variations in the pressure traces and could not predict the strong CCV of the in-cylinder pressures represented by the large grey envelop. In comparison, the LES results in Figure 2b can reproduce the strong cycle-to-cycle variations in the in-cylinder pressure trace very well, although there are some differences in the predicted and experimental lower limit (almost misfire). The turbulence resolution parameter (M) proposed by Pope [41] was evaluated at different crank angles to assess the applicability of the LES results. M is defined as the ratio of the sub-grid scale (SGS) kinetic energy to the total kinetic energy. M of 0 denotes a direct numerical simulation (DNS) simulation where all the turbulence scales are resolved, while the M of 1 denotes a RANS simulation where no turbulence scales are resolved and all are modelled. It is recommended that M should be less than 0.2 for an adaptive LES simulation [41]. In this study, M is significantly lower than 0.2 in the LES simulations [28]. The spark plug region shows a relatively higher turbulence resolution parameter M, in particular at the near wall region of the ground electrode with M at 0.14, due to complex geometric features. The majority of the region around the spark plug gap has a turbulence resolution parameter well below 0.08. Figure 2 shows the in-cylinder pressure profiles from the experiments, RANS and LES simulations. The grey envelop highlights the large variations in the experimental results (100 cycles). The dark grey dash curve shows the averaged pressure of the experimental results. As shown in Figure 2a, the RANS simulation results show little variations in the pressure traces and could not predict the strong CCV of the in-cylinder pressures represented by the large grey envelop. In comparison, the LES results in Figure 2b can reproduce the strong cycle-to-cycle variations in the in-cylinder pressure trace very well, although there are some differences in the predicted and experimental lower limit (almost misfire).

Validations of RANS and LES Results with Experimental Measurements
In order to quantitatively compare the experimental and simulation results, Figure 3 shows the IMEP histograms of the consecutive cycles from experiments (100 cycles), LES simulations (30 cycles) and RANS simulations (10 cycles). The corresponding average value x, standard deviation σ and COV of IMEP are also included in Figure 3. Overall, the distribution pattern of the IMEP among the cycles is similar between the experiments and LES simulations and shows strong CCV of IMEP. Both experiments and LES simulations have more cycles with higher IMEP values at the right part of the distribution. However, it is noted that the current LES modelling could not predict the weak cycles with significantly lower IMPE values spreading over a wider region in the left part of the distribution. Overall, the predicted average IMEP of 6.67 bar by LES is higher than the measured average IMEP of 5.43 bar. In addition, the COV in IMEP is also under-predicted at 5.25% by LES simulations, compared to the measured value of 11.81%. In comparison, the average IMEP is well predicted by the RANS simulations at 5.8 bar, compared to the experimental result at In order to quantitatively compare the experimental and simulation results, Figure 3 shows the IMEP histograms of the consecutive cycles from experiments (100 cycles), LES simulations (30 cycles) and RANS simulations (10 cycles). The corresponding average value x, standard deviation σ and COV of IMEP are also included in Figure 3. Overall, the distribution pattern of the IMEP among the cycles is similar between the experiments and LES simulations and shows strong CCV of IMEP. Both experiments and LES simulations have more cycles with higher IMEP values at the right part of the distribution. However, it is noted that the current LES modelling could not predict the weak cycles with significantly lower IMPE values spreading over a wider region in the left part of the distribution. Overall, the predicted average IMEP of 6.67 bar by LES is higher than the measured average IMEP of 5.43 bar. In addition, the COV in IMEP is also under-predicted at 5.25% by LES simulations, compared to the measured value of 11.81%. In comparison, the average IMEP is well predicted by the RANS simulations at 5.8 bar, compared to the experimental result at 5.43 bar. However, the RANS simulations produced negligible COV in IMEP at 0.08%, and failed to reproduce the cyclic variations, as observed in the experiments and LES.    In order to further evaluate the performance of the LES and RANS simulations, Figures 4 and 5 compare the in-cylinder peak pressure (PP) and the peak pressure rise rate (PPRR) with the corresponding crank angles, respectively. Similarly, the average value x, standard deviation σ and COV are also included in these figures. Figure 4 reveals the nonlinear relationship between PP and the corresponding crank angle of the experimental results and the presence of two distinct combustion types, i.e., weak combustion at the bottom-left and strong combustion at the upper-right. According to the heat release analysis of the experimental data [42], the weak combustion cycles are dominated by the slow flame propagation with less/no subsequent auto-ignition combustion, leading to increased PP, as well as more complete combustion with the crank angle. In comparison, the strong combustion cycles are dominated by auto-ignition heat release, and the peak pressure is produced by the auto-ignition combustion. Therefore, the peak pressure is reduced with the auto-ignition event moving towards to the expansion stroke. This phenomenon was also observed in typical spark ignition engines [43]. Although the LES simulations show promising agreement of the average PP (44.35 bar) and the corresponding COV (22.2%) with experimental observations (42.11 bar and 26.2%, respectively), the detailed comparison of the scatter distribution with the experiments indicates the great challenge of modelling the nonlinear phenomenon of cyclic hybrid combustion. As shown in the figure, the LES simulations could predict PP and the corresponding crank angle very well for the strong cycles. However, the weak combustion cycles were not captured with the current LES models. The reasons could be the overestimation of the flame propagation speed and auto-ignition combustion speed at extreme conditions near the experimental lower limit, as shown in Figure 2.   The two distinct combustion types can also be observed in Figure 5. The weak combustion cycles with peak pressure rise rate (PPRR) below 2 bar/ • CA are located before TDC, while the strong cycles with PPRR ranging between 1 and 20 bar/ • CA are located after TDC. Although there are noticeable differences between the predicted and measured average PPRR and COV, the LES simulations were able to predict the two distinct combustion types. However, the LES produced fewer number of weak cycles and the corresponding crank angles closer to TDC than the experimental observations. In addition, the LES predicted more retarded crank angles of PPRR for some strong cycles, indicating very the late auto-ignition processes for those cycles in LES.  The two distinct combustion types can also be observed in Figure 5. The weak bustion cycles with peak pressure rise rate (PPRR) below 2 bar/°CA are located b TDC, while the strong cycles with PPRR ranging between 1 and 20 bar/°CA are lo after TDC. Although there are noticeable differences between the predicted and mea average PPRR and COV, the LES simulations were able to predict the two distinct bustion types. However, the LES produced fewer number of weak cycles and the sponding crank angles closer to TDC than the experimental observations. In additio LES predicted more retarded crank angles of PPRR for some strong cycles, indicating the late auto-ignition processes for those cycles in LES.
As shown in Figures 4 and 5, the RANS simulations could not capture the two d combustion types of the SI-CAI hybrid combustion. The predicted average PP (39.6 is close to the experimental value (42.11 bar), but the predicted PPRR (1.63 bar) is s cantly lower than the experimental result (4.04 bar/°CA). The predicted COVs of P PPRR by RANS are far lower than the experimental observations. Figure 6a shows the average crank angles with the mass fraction of fuel burn 10%, 50% and 90% (CA10, CA50 and CA90). The predicted CA50 by LES is slightly  Figure 6a shows the average crank angles with the mass fraction of fuel burned at 10%, 50% and 90% (CA10, CA50 and CA90). The predicted CA50 by LES is slightly later than the experimental value and the predicted CA10 and CA90 are earlier than the experimental values, indicating slower early flame propagation but faster subsequent auto-ignition dominated combustion process by the LES. The predicted faster auto-ignition process in the expansion stroke also contributed to higher average IMEP in the LES simulations compared to the experimental results, as shown in Figure 3. The predicted CA50 and CA90 of RANS simulations are close to the experimental values, although the CA10 is earlier than the experimental result. It can be observed in Figure 6b that the measured cyclic variation is gradually enlarged as the combustion process proceeds and the standard deviation (σ) is increased from 2.58 • CA for CA10 to 12.34 • CA for CA90. This trend was fairly well reproduced by the current LES results, with the predicted σ increasing from 2.2 • CA for CA10 to 7.78 • CA for CA90. The slightly lower predicted σ in CA10 and CA90 indicate that both the spark ignition/initial flame development and late auto-ignition combustion models in the current LES need to be improved to capture their sensitivities to changes in thermofluidic properties in the cylinder. Compared to the LES, the RANS simulations could not reproduce the enlarged cyclic changes during the hybrid combustion process.
In order to evaluate the variation in the transition of combustion mode from spark ignition to controlled auto-ignition mode, the crank angle of the transition (CAT pressure ) was calculated based on the pressure trace. As shown in Figure 7a, the CAT pressure is defined as the transition point where the local first order derivative of the pressure trace is at the minimum, while the second order derivative just exceeds 0.1. For the non-transitioning cycles in which no transition point is detected, the CAT is simply denoted as 30 • CA. This was chosen for the purpose of plotting and no CAT was detected after 30 • CA. was fairly well reproduced by the current LES results, with the predicted σ increasing from 2.2 °CA for CA10 to 7.78 °CA for CA90. The slightly lower predicted σ in CA10 and CA90 indicate that both the spark ignition/initial flame development and late auto-ignition combustion models in the current LES need to be improved to capture their sensitivities to changes in thermofluidic properties in the cylinder. Compared to the LES, the RANS simulations could not reproduce the enlarged cyclic changes during the hybrid combustion process. In order to evaluate the variation in the transition of combustion mode from spark ignition to controlled auto-ignition mode, the crank angle of the transition (CATpressure) was calculated based on the pressure trace. As shown in Figure 7a, the CATpressure is defined as the transition point where the local first order derivative of the pressure trace is at the minimum, while the second order derivative just exceeds 0.1. For the non-transitioning cycles in which no transition point is detected, the CAT is simply denoted as 30 °CA. This was chosen for the purpose of plotting and no CAT was detected after 30 °CA.
In order to validate the applicability of the calculation of the transition point based on cylinder pressure trace, the LES simulation results were analyzed and the crank angle with 1% cells auto-ignited is defined as the transition crank angle (CATignition) based on the ignition cells. The results are compared with the CATpressure, as shown in Figure 7b. The    In order to validate the applicability of the calculation of the transition point based on cylinder pressure trace, the LES simulation results were analyzed and the crank angle with 1% cells auto-ignited is defined as the transition crank angle (CAT ignition ) based on the ignition cells. The results are compared with the CAT pressure , as shown in Figure 7b. The results indicate that the cylinder pressure-based approach shows almost the same results of CAT as the autoignition cells-based approach. Therefore, the pressure-based approach was applied to analyze the CAT in both experiments and simulations for consistency. Figure 8 shows

Analysis of the Cyclic Phenomenon of Hybrid Combustion
The above critical validation study shows that the adopted LES simulations could capture the basic combustion characteristics of the SI-CAI hybrid combustion observed in the experiments. More specifically, two distinct combustion types (strong/weak cycles), the variation in combustion mode transitioning, and the enlarged cyclic phenomenon throughout the combustion process can be well captured in the current LES modeling. Thus, the LES simulation results will be further explored in this section to investigate and understand the observed cyclic phenomenon of the SI-CAI hybrid combustion.
In order to analyze the cyclic variation in the initial flame development, a spherical zone with 20 mm diameter around the spark plug gap is defined as the spark zone to capture the initial conditions. The calculation of the COVs of local velocity magnitude (VM), temperature, residual gas fraction (RGF) and pressure was performed for the whole combustion chamber, as well as the spark zone. It is found that the local temperature, RGF and pressure showed little variations (COVs < 0.8%) but the VM shows significant variations with COV of 8.66% in the whole chamber and 22.75% in the spark zone.
Furthermore, a correlation study was performed for the in-cylinder velocity magnitude, temperature and RGF before spark ignition timing against CA10, CA50 and CA90 to identify the most influential factors affecting the CCVs of the hybrid combustion process. Strong correlations are found between the VM in spark zone at −50 °CA (just before the spark ignition timing) and CA10, CA50 and CA90, as shown in Figure 9. In particular, CA10 is strongly dependent on the local VM in the spark zone and the Pearson correlation coefficient (r) was around 0.8. It is also noted that the slops of the fitting curves gradually increase from CA10 to CA90, indicating the enhanced impact of local VM in the spark

Analysis of the Cyclic Phenomenon of Hybrid Combustion
The above critical validation study shows that the adopted LES simulations could capture the basic combustion characteristics of the SI-CAI hybrid combustion observed in the experiments. More specifically, two distinct combustion types (strong/weak cycles), the variation in combustion mode transitioning, and the enlarged cyclic phenomenon throughout the combustion process can be well captured in the current LES modeling. Thus, the LES simulation results will be further explored in this section to investigate and understand the observed cyclic phenomenon of the SI-CAI hybrid combustion.
In order to analyze the cyclic variation in the initial flame development, a spherical zone with 20 mm diameter around the spark plug gap is defined as the spark zone to capture the initial conditions. The calculation of the COVs of local velocity magnitude (VM), temperature, residual gas fraction (RGF) and pressure was performed for the whole combustion chamber, as well as the spark zone. It is found that the local temperature, RGF and pressure showed little variations (COVs < 0.8%) but the VM shows significant variations with COV of 8.66% in the whole chamber and 22.75% in the spark zone.
Furthermore, a correlation study was performed for the in-cylinder velocity magnitude, temperature and RGF before spark ignition timing against CA10, CA50 and CA90 to identify the most influential factors affecting the CCVs of the hybrid combustion process. Strong correlations are found between the VM in spark zone at −50 • CA (just before the spark ignition timing) and CA10, CA50 and CA90, as shown in Figure 9. In particular, CA10 is strongly dependent on the local VM in the spark zone and the Pearson correlation coefficient (r) was around 0.8. It is also noted that the slops of the fitting curves gradually increase from CA10 to CA90, indicating the enhanced impact of local VM in the spark zone on the combustion process. It should be noted that the correlation of the average VM in the whole cylinder is relatively weak, with r below 0.7. All the other parameters (temperature and RGF) before spark ignition timing show very weak correlation (Pearson correlation coefficient r below 0.3) with combustion parameters. Although the average in-cylinder temperature at 670 °CA has very weak correlations with CA10, CA50 and CA90, Figure 10 shows strong correlations of the average temperature at TDC (T TDC) with CA50, CA90 and CAT with a Pearson coefficient over 0.94, which indicates the dominated impact of the induced thermal condition by the early flame propagation on the auto-ignition combustion at a later stage of hybrid combustion. It should be noted that CA10 is earlier than TDC; therefore, the correlation of T TDC with CA10 is not shown here.
The correlation study reveals that the early phase flame propagation (CA10) of the hybrid combustion is controlled by the VM around the spark plug. The flame propagation induces thermal conditions, then dominates the subsequent combustion process (CA50 and CA90) as well as mode transitioning (CAT), as the auto-ignition process is more sensitive to temperature.  Although the average in-cylinder temperature at 670 • CA has very weak correlations with CA10, CA50 and CA90, Figure 10 shows strong correlations of the average temperature at TDC (T TDC ) with CA50, CA90 and CAT with a Pearson coefficient over 0.94, which indicates the dominated impact of the induced thermal condition by the early flame propagation on the auto-ignition combustion at a later stage of hybrid combustion. It should be noted that CA10 is earlier than TDC; therefore, the correlation of T TDC with CA10 is not shown here. Although the average in-cylinder temperature at 670 °CA has very weak correlations with CA10, CA50 and CA90, Figure 10 shows strong correlations of the average temperature at TDC (T TDC) with CA50, CA90 and CAT with a Pearson coefficient over 0.94, which indicates the dominated impact of the induced thermal condition by the early flame propagation on the auto-ignition combustion at a later stage of hybrid combustion. It should be noted that CA10 is earlier than TDC; therefore, the correlation of T TDC with CA10 is not shown here.
The correlation study reveals that the early phase flame propagation (CA10) of the hybrid combustion is controlled by the VM around the spark plug. The flame propagation induces thermal conditions, then dominates the subsequent combustion process (CA50 and CA90) as well as mode transitioning (CAT), as the auto-ignition process is more sensitive to temperature. In order to understand the combustion process that is influenced by in-cylinder flow motion and temperature, five typical combustion cycles with both strong and weak combustion processes were selected for further analysis. The corresponding mass fraction burned (MFB) traces are shown in the Figure 11. As shown in the figure, the combustion  The correlation study reveals that the early phase flame propagation (CA10) of the hybrid combustion is controlled by the VM around the spark plug. The flame propagation induces thermal conditions, then dominates the subsequent combustion process (CA50 and CA90) as well as mode transitioning (CAT), as the auto-ignition process is more sensitive to temperature.
In order to understand the combustion process that is influenced by in-cylinder flow motion and temperature, five typical combustion cycles with both strong and weak combustion processes were selected for further analysis. The corresponding mass fraction burned (MFB) traces are shown in the Figure 11. As shown in the figure, the combustion process is faster for Cycles 28, 6 and then 13. Cycles 14 and 25 are the slowest with much later combustion phasing (CA10 to CA90), as well as CAT.  The distributions of flow magnitudes and vectors at 670 °CA (just before spark ignition timing) are shown in Figure 12. The sections used for the distribution plots are also marked in the figure. It is visible that the flow velocity magnitude in the area around the spark plug is lower (in blue color) for Cycles 14 and 25. The flow velocity vectors captured in-cylinder flow motions and they are very different at the spark plug region among the cycles, which implies the potential impact on the flame kernel formation and initial development of the flame propagation process. In Cycle 12, there is a strong flow across the spark plug region towards the upper part of the horizontal plots. An anti-clockwise tumble towards the spark plug gap is also observed in the vertical plot. Similarly, strong flow motions are observed across the spark plug region towards the upper part of the cylinder in the horizontal plot of Cycle 28, while a clockwise tumble motion towards the spark plug gap is formed in the vertical plot. In comparison, a relatively weaker flow motion pointing to the spark plug is observed in Cycle 13, while the clockwise tumble in the vertical plot is slightly away from the spark plug gap. However, there is no coordinated flow motions across the spark plug region in Cycles 14 and 15. The distributions of flow magnitudes and vectors at 670 • CA (just before spark ignition timing) are shown in Figure 12. The sections used for the distribution plots are also marked in the figure. It is visible that the flow velocity magnitude in the area around the spark plug is lower (in blue color) for Cycles 14 and 25. The flow velocity vectors captured in-cylinder flow motions and they are very different at the spark plug region among the cycles, which implies the potential impact on the flame kernel formation and initial development of the flame propagation process. In Cycle 12, there is a strong flow across the spark plug region towards the upper part of the horizontal plots. An anti-clockwise tumble towards the spark plug gap is also observed in the vertical plot. Similarly, strong flow motions are observed across the spark plug region towards the upper part of the cylinder in the horizontal plot of Cycle 28, while a clockwise tumble motion towards the spark plug gap is formed in the vertical plot. In comparison, a relatively weaker flow motion pointing to the spark plug is observed in Cycle 13, while the clockwise tumble in the vertical plot is slightly away from the spark plug gap. However, there is no coordinated flow motions across the spark plug region in Cycles 14 and 15.
As a result of the in-cylinder flow motions, the combustion process varies significantly, as captured in Figure 13, which shows the iso-surface of the flame front and auto-ignition sites at different crank angles. The iso-surface of the flame surface density (100 m −1 ) is plotted in dark blue to capture the flame front, while the iso-surface of the auto-ignition tendency (0.98) is plotted in brown to indicate the potential auto-ignition sites in the next time step. In this way, the interactions of flame propagation on the auto-ignition at specific crank angles can be analyzed. It should be noted that these auto-ignition sites are not necessary the first auto-ignition sites in a specific cycle. As a result of the in-cylinder flow motions, the combustion process varies significantly, as captured in Figure 13, which shows the iso-surface of the flame front and autoignition sites at different crank angles. The iso-surface of the flame surface density (100 m −1 ) is plotted in dark blue to capture the flame front, while the iso-surface of the autoignition tendency (0.98) is plotted in brown to indicate the potential auto-ignition sites in the next time step. In this way, the interactions of flame propagation on the auto-ignition at specific crank angles can be analyzed. It should be noted that these auto-ignition sites are not necessary the first auto-ignition sites in a specific cycle.
The combined information in Figures 12 and 13 provide valuable information on the interactions of in-cylinder flow motions and combustion processes. The results indicate that the initial flame front is directly affected by the flow motions. For example, in Cycle 6, the early flame kernel and subsequent flame front are driven by the flow motions towards the upper part of the cylinder in the horizontal plot (as shown in Figure 12). At TDC, the flame front already spreads to the cylinder wall of the upper part of the cylinder, while the flame front is in the middle of the lower part of the combustion chamber, as shown in Figure 13. The further development of the flame propagation process triggers the auto-ignitions (in brown color) at the rim of the flame front in the upper part at 10 °CA.
A similar flame development process can be observed in Cycle 13, where the initial flame kernel at −40 to −30 °CA is blown to the upper part of the chamber and the majority of the flame front is kept in the upper part until 10 °CA, when the flame propagates into the lower part of the chamber. However, it is evident that the flame propagation process in Cycle 13 is slower compared to Cycle 6. This is caused by the relatively weaker flow motion in the central region for Cycle 13 compared to Cycle 6, as shown in Figure 12  motion across the spark plug. It is also obvious that the flame front of Cycle 25 (e.g., at TDC and 10 °CA) has less winkles and is close to a spherical shape due to less disturbance from the flow fields. However, this leads to the slowest combustion process among the five cycles, as shown in Figure 11. The strong flow motion in Cycle 28 crosses the spark plug region and moves towards the upper right corner, as shown in Figure 12. As a result, the flame front is driven towards the upper right corner from the very beginning of flame kernel formation. It is particularly obvious at −10 °CA when the flame front is biased towards to the upper right corner of the chamber. Figure 13 also shows important information on the mode transitioning from flame propagation mode to auto-ignition. Interactions of the flame and auto-ignition sites are well captured in Cycle 6, Cycle 14 and Cycle 25, as the selected crank angles for plotting are very close to the CAT. More specifically, it can be observed that the auto-ignition sites are located very close to the flame front in Cycle 6. In Cycle 14, the auto-ignition sites at 30 °CA are very close to the cylinder walls (in particular, the cylinder head and cylinder liner), while the flame fronts are much more weakened. In comparison, the auto-ignition sites at 30 °CA in Cycle 25 are closer to the piston top, while the flame front spread more to the cylinder wall, indicating a slower combustion speed and later transitioning to the auto-ignition mode, as shown in Figure 11.
It should be noted that in Cycles 13 and 28, the CAT is away from the selected crank angles for plotting. Therefore, the flame is much more weakened, while the auto-ignition combustion also approaches the end and only the mixtures close to the wall are left to be auto-ignited at 20 °CA for Cycle 13 and 10 °CA for Cycle 28. The auto-ignition process is affected by local thermal conditions. Figure 14 shows the temperature distributions at −50 °CA and TDC, and the same horizontal section as Figure  12 is used for the plotting. Overall, the temperature distribution at −50 °CA shows very limited differences among the cases, as the temperature range of the plots is 520 K to 580 K. In addition, it was found that a higher temperature in the center does not necessarily produce a fast flame propagation. For example, Cycle 14 has a relatively higher temperature in the center compared to Cycle 6, while its flame propagation is slower than Cycle 6, as shown in Figures 11 and 13. This also indicates that the early flame propagation stage is less affected by the cylinder thermal condition. Instead, it is more influenced by the incylinder flow motions, as analyzed above.
However, the development of the early flame propagation stage inevitably influences the in-cylinder thermal conditions. As shown in the plots of temperature distributions at TDC, the area of the high temperature region of the two slow cycles (Cycles 14 and 25) is smaller compared to the faster cycles. As the result, the temperature of the unburnt mixture in the outer region of the slow cycles is also lower. For example, the temperature of the near-wall region of Cycles 14 and 25 is around 600 K (in purple), while the near-wall temperature is higher for Cycles 6, 13 and 28. The cylinder temperature distributions, in turn, affect the formation of auto-ignition sites. Therefore, Cycles 14 and 25 have relatively A similar flame development process can be observed in Cycle 13, where the initial flame kernel at −40 to −30 • CA is blown to the upper part of the chamber and the majority of the flame front is kept in the upper part until 10 • CA, when the flame propagates into the lower part of the chamber. However, it is evident that the flame propagation process in Cycle 13 is slower compared to Cycle 6. This is caused by the relatively weaker flow motion in the central region for Cycle 13 compared to Cycle 6, as shown in Figure 12. In Cycles 14 and 25, the flow velocity is low in the central region around the spark plug and no clear flow motion is formed across the spark plug region. As a result, the flame kernel formation and flame development are weaker compared to Cycle 6 and Cycle 13. It is visible that Cycles 14 and 25 produce a smaller flame area, in particular from TDC, compared to Cycles 6 and 13. In addition, the flame fronts in Cycles 14 and 25 are positioned in the center of the chamber around the spark plug, due to less impact by the bulk flow motion across the spark plug. It is also obvious that the flame front of Cycle 25 (e.g., at TDC and 10 • CA) has less winkles and is close to a spherical shape due to less disturbance from the flow fields. However, this leads to the slowest combustion process among the five cycles, as shown in Figure 11. The strong flow motion in Cycle 28 crosses the spark plug region and moves towards the upper right corner, as shown in Figure 12. As a result, the flame front is driven towards the upper right corner from the very beginning of flame kernel formation. It is particularly obvious at −10 • CA when the flame front is biased towards to the upper right corner of the chamber. Figure 13 also shows important information on the mode transitioning from flame propagation mode to auto-ignition. Interactions of the flame and auto-ignition sites are well captured in Cycle 6, Cycle 14 and Cycle 25, as the selected crank angles for plotting are very close to the CAT. More specifically, it can be observed that the auto-ignition sites are located very close to the flame front in Cycle 6. In Cycle 14, the auto-ignition sites at 30 • CA are very close to the cylinder walls (in particular, the cylinder head and cylinder liner), while the flame fronts are much more weakened. In comparison, the auto-ignition sites at 30 • CA in Cycle 25 are closer to the piston top, while the flame front spread more to the cylinder wall, indicating a slower combustion speed and later transitioning to the auto-ignition mode, as shown in Figure 11.
It should be noted that in Cycles 13 and 28, the CAT is away from the selected crank angles for plotting. Therefore, the flame is much more weakened, while the auto-ignition combustion also approaches the end and only the mixtures close to the wall are left to be auto-ignited at 20 • CA for Cycle 13 and 10 • CA for Cycle 28.
The auto-ignition process is affected by local thermal conditions. Figure 14 shows the temperature distributions at −50 • CA and TDC, and the same horizontal section as Figure 12 is used for the plotting. Overall, the temperature distribution at −50 • CA shows very limited differences among the cases, as the temperature range of the plots is 520 K to 580 K. In addition, it was found that a higher temperature in the center does not necessarily produce a fast flame propagation. For example, Cycle 14 has a relatively higher temperature in the center compared to Cycle 6, while its flame propagation is slower than Cycle 6, as shown in Figures 11 and 13. This also indicates that the early flame propagation stage is less affected by the cylinder thermal condition. Instead, it is more influenced by the in-cylinder flow motions, as analyzed above. later auto-ignition timing, as indicated in both Figures 11 and 13. In addition, more mixture is burned at CAT in order to trigger the auto-ignition, due to the delayed combustion event towards the expansion stroke.

Conclusions
In this study, LES and RANS simulations were employed to model the cyclic phenomenon of SI-CAI hybrid combustion. A detailed analysis was carried out on the incylinder pressure traces, IMEP, in-cylinder peak pressure (PP), peak pressure rise rate (PPRR), the crank angles with fuel mass fraction burned at 10%, 50%, 90% and mode transition to evaluate the capability of LES and RANS simulations to reproduce the cyclic variations in hybrid combustion. LES simulations were then further used to investigate and understand the observed cyclic phenomenon of SI-CAI hybrid combustion. The main conclusions can be summarized as follows.
(1) The adopted LES simulation is capable of predicting the cyclic variations in hybrid combustion observed in the experiments, with strong CCVs of IMEP and very similar However, the development of the early flame propagation stage inevitably influences the in-cylinder thermal conditions. As shown in the plots of temperature distributions at TDC, the area of the high temperature region of the two slow cycles (Cycles 14 and 25) is smaller compared to the faster cycles. As the result, the temperature of the unburnt mixture in the outer region of the slow cycles is also lower. For example, the temperature of the near-wall region of Cycles 14 and 25 is around 600 K (in purple), while the near-wall temperature is higher for Cycles 6, 13 and 28. The cylinder temperature distributions, in turn, affect the formation of auto-ignition sites. Therefore, Cycles 14 and 25 have relatively later auto-ignition timing, as indicated in both Figures 11 and 13. In addition, more mixture is burned at CAT in order to trigger the auto-ignition, due to the delayed combustion event towards the expansion stroke.

Conclusions
In this study, LES and RANS simulations were employed to model the cyclic phenomenon of SI-CAI hybrid combustion. A detailed analysis was carried out on the incylinder pressure traces, IMEP, in-cylinder peak pressure (PP), peak pressure rise rate (PPRR), the crank angles with fuel mass fraction burned at 10%, 50%, 90% and mode transition to evaluate the capability of LES and RANS simulations to reproduce the cyclic variations in hybrid combustion. LES simulations were then further used to investigate and understand the observed cyclic phenomenon of SI-CAI hybrid combustion. The main conclusions can be summarized as follows.
(1) The adopted LES simulation is capable of predicting the cyclic variations in hybrid combustion observed in the experiments, with strong CCVs of IMEP and very similar histogram distribution patterns. However, RANS simulations fail to reproduce the cyclic characteristics of the hybrid combustion at the chosen engine operating conditions. (2) The analysis shows the nonlinear relationship of PP/PPRR against the corresponding crank angles in the experimental results and the presence of two distinct combustion types (strong/weak cycles). The measured cyclic variation is gradually enlarged as the combustion process proceeds from CA10 to CA90. (3) The LES simulations are able to produce an excellent agreement between the predicted average PP (44.35 bar) and the COV (22.2%) and experimental observations (42.11 bar and 26.2%, respectively). The LES simulations are also able to predict the two distinct combustion types and CCVs in the PPRR-crank angle plot, and the enlarged cyclic variations from CA10 to CA90. (4) The CAT pressure, defined as the transition point where the local first order derivative of the pressure trace is at the minimum, while the second order derivative just exceeds 0.1, shows very similar results of the crank angle with 1% cells auto-ignited (CAT ignition) in LES based on autoignition cells. The distribution of CATs of the current cycle and next cycle in LES shows a very good agreement with experimental results, indicating excellent capture of the combustion mode transitioning with the LES. (5) However, the detailed comparison between the LES and experimental data also highlights the challenges of modelling the nonlinear phenomenon of SI-CAI hybrid combustion. For example, the weak cycles are not well captured with the existing LES modelling, which produce late CAT, higher average IMEP and lower COV compared to the experiments. The reasons could be the overestimation of the flame propagation speed and auto-ignition combustion speed at extreme conditions near the experimental lower limit. (6) The correlation study reveals that the local velocity magnitude (VM) around the spark plug before spark ignition timing has a large effect on the early flame propagation and CA10, while the induced thermal conditions by early flame propagation dominates the subsequent auto-ignition combustion in the late stage, as evidenced by the strong correlation between the temperatures at TDC and CA50, CA90 and CAT.