A Combined Numerical and Experimental Investigation of Cycle-to-Cycle Variations in an Optically Accessible Spark-Ignition Engine

A combined numerical and experimental investigation is carried out to analyze the cycle-to-cycle variations (CCV) in an optically accessible spark-ignition engine with port fuel injection. A stable and an unstable operating condition is considered. Well-established turbulence, combustion, and ignition models are employed in the large-eddy simulations (LES). High-speed measurements of the velocity field via particle image velocimetry and flame imaging in the tumble plane are conducted in the experiments. A detailed comparison between LES and experiments is carried out, including the in-cylinder pressure, the flow fields, the spatial flame distribution, and the fields conditioned on fast and slow cycles. Good agreement is achieved for the variables considering all cycles; yet, some discrepancies are observed for the conditionally averaged quantities. A systematic quantitative correlation analysis between the selected influencing variables and the CCV is presented, in which the influencing variables are extracted from different length scales (r = 3 mm, 12 mm, and 43 mm) and the CCV are distinguished between the early flame kernel development and later flame propagation. Even though the most relevant influencing parameters are different for the two operating conditions, the location of the coherent vortex structure is found to be important for the CCV of both cases.


Introduction
Environmental pollution and global warming have prompted the need for the use of highefficiency and low-emission propulsion technologies. Besides electric powertrains in vehicles, internal combustion engines (ICE) will still be used in the near future due to the advanced state of ICE technologies including their reliability, durability, and the existing mature supply chain, manufacturing infrastructure, and recycling facilities (Leach et al. 2020). Indeed, ICE operated with non-fossil fuels (e.g. e-fuels, biofuels, or hydrogen) will play an important role in the future transport sector, especially in combination with electrification for heavy-duty vehicles (Onorati et al. 2022). To reduce greenhouse gas emissions, thermal efficiencies are essential for ICE. Concepts such as exhaust gas recirculation (EGR) (Fontana and Galloni 2010;Kargul et al. 2019) and ultra lean-burn (Luszcz et al. 2018;Ye et al. 2021) are promising to improve the efficiency of spark-ignition (SI) engines. However, cycle-to-cycle variations (CCV), which are characterized by the nonrepeatability of combustion events for consecutive engine cycles under the same nominal operating conditions (Heywood 2018), represent a limit to the achievable dilution level of stable operations, and thus limit the improvements of engine efficiency and emissions. The predictive description of CCV and fundamental understanding of its causes are still limited but are essential for the development of high-efficiency ICE.
The significant enhancement of engine diagnostics enables detailed experimental investigations of the in-cylinder process and CCV (Buschbeck et al. 2012;Zeng et al. 2015;Bode et al. 2017;Zeng et al. 2019). Important sources for CCV were identified, which include the kinetic energy of the flow field at stoichiometric operating conditions (Buschbeck et al. 2012), the flow-spray interaction for direct-injection spark-ignition (DISI) engines (Zeng et al. 2015), and flow features in specific regions (Bode et al. 2017;Zeng et al. 2019). Experimental studies typically focus on the global quantities, such as the incylinder pressure or heat release rate, and two-dimensional quantities, such as the mixture composition, velocity, and flame location in selected planes, due to the limited access to the three-dimensional fields. However, the availability of in-cylinder optical measurements allows a detailed comparison with numerical studies, where the three-dimensional fields are available. Such combined numerical and experimental investigations are promising to gain further understanding of CCV.
The advance of large-eddy simulations (LES) as an inherently unsteady simulation technique reflecting the stochastic nature of turbulent flows, enables numerical studies on CCV, such as qualitative analysis considering selected fast and slow cycles and systematic quantitative correlation analysis. Granet et al. (2012) carried out LES of stable and unstable engine operation conditions and demonstrated that LES can be used to distinguish a stable operation condition from an unstable one. Zembi et al. (2019) performed LES of three engine operation conditions with different values of relative air-fuel ratio. It was shown that LES has the capability to predict the transition from stable to unstable lean operation conditions. Enaux et al. (2011) reproduced the range of the variations of the in-cylinder pressure observed in experiments using LES. By comparing the selected typical fast and slow cycles, it was found that the velocity fluctuations at the spark plug play a significant role for early flame kernel development and the overall combustion. Using LES of a stable engine operation condition, Zembi et al. (2022) also showed that valuable indicators of the combustion rate can be derived based on the local flow fields around the spark plug. Under the investigated conditions, greater velocities lead to faster combustion processes. Zhao et al. (2018) conducted correlation studies using LES and showed that the velocity field close to the spark plug determines differences in flame propagation from cycle to cycle. Truffin et al. (2015) performed correlation studies to systematically investigate the influence of the global (referring to the whole cylinder) and local (referring to only the region close to the spark plug) variables at spark timing, such as the averaged local and global temperature, pressure, and flow features, on the indicated mean effective pressure (IMEP). Significant correlations were identified and it was found that the most relevant influencing variables are notably different for varying operating conditions. Using LES of two relevantly different engines, d 'Adamo et al. (2019) demonstrated significant correlations between the global as well as the local variables and the duration of the first 1% of mass fraction burned (MFB1). They also found that the most influencing factors for determining CCV are case specific. In the literature, reactive LES results for CCV investigations are commonly validated with experiments using global variables, such as the in-cylinder pressure, or mean velocity fields as well as flame images. However, detailed comparisons between LES and experiments, especially including both the velocity fields and flame images conditioned to the fast and slow cycles, are rare. The correlation analyses of CCV usually focus on a single combustion parameter, such as IMEP, MFB1, or the maximal in-cylinder pressure. In addition, the large-scale flow features are typically represented by global quantities, such as the tumble and swirl ratio based on the center of in-cylinder mass considering gas motion inside the whole cylinder, which are not suitable for some engine geometries or operating conditions. A systematic analysis of multiple aspects is expected to further increase the understanding of CCV, such as the separate evaluation of combustion phases in combination with their main different governing effects and proper description of the coherent vortex structure.
The present study aims to systematically assess the influence of variables extracted from different length scales on the CCV of different combustion phases. For this purpose, a combined numerical and experimental investigation of an optically accessible SI engine with a detailed comparison between LES and experiments has been carried out. The comparison includes the flow fields, the spatial flame distribution, and especially the fields conditioned on the fast and slow cycles. The method presented by Buhl et al. (2017) has been employed to detect the three-dimensional vortex structure in the complex in-cylinder flows. Stable and unstable operating conditions have been considered.
The remainder of the paper is organized in the following manner. First, the experimental methods and operating conditions are given in Sect. 2, followed by the numerical methods in Sect. 3. Then, a detailed comparison between LES and experiments, including the incylinder pressure, the velocity fields, and the flame propagation, is presented in Sect. 4. Afterwards, CCV of the early flame kernel development and the flame propagation are analyzed along with their sources in Sect. 5. Finally, the paper finishes with a summary and conclusions.

Experimental Methods and Operating Conditions
The engine considered in this study is the well-characterized single-cylinder four-valve optically accessible engine at TU Darmstadt (Baum et al. 2014). The research engine employs a Bowditch piston extension, a flat quartz glass piston surface, and quartz glass cylinder liner to allow optical access from the bottom and the sides of the cylinder. The engine's spray-guided cylinder head configuration with a pent-roof and a bore and stroke each of 86 mm results in a compression ratio of 8.7:1. Using two large plenums upstream 1 3 of the intake pipe, vibrations and acoustic oscillations are controlled to provide precise boundary conditions. Two part-load operating conditions with an intake pressure of 0.4 bar and an engine speed of 1500 rpm were selected to represent a stable (stab) and an unstable (unst) condition for fired engine experiments of stoichiometric iso-octane fuel introduced via port fuel injection. The stab case is characterized by the continuous firing of a stoichiometric fuel-air mixture. For the unst case, nitrogen and carbon-dioxide are mixed with the intake air to simulate 12.9% external EGR. To ensure a constant amount of total EGR, the unst case employed skip-firing in which only every 7th cycle was ignited in order to flush out residual exhaust gases due to internal EGR before the next ignition cycle. Furthermore, the unst case has an earlier ignition timing. The ignition timing of both conditions was optimized such that the average cycle's 50% mass fraction burned (MFB50) occurred at 8 crank-angle-degrees (CAD). In this study, 0 CAD stands for compression top dead center. The resulting coefficient of variation (COV) of the maximum in-cylinder pressure for the unst case is 12.8% compared to 3.51% for the stab case. Although the stab case has a higher amount of residual gases due to the internally induced EGR by the valve overlap, the condition is nevertheless more stable than the carefully controlled external EGR case due to the different spark timings. The later spark timing of the stab case results in a stable operation and an average exhaust temperature of more than 250 • C greater than the unst case. Engine specifications and operating conditions are summarized in Table 1.
High-speed planar Mie scattering of oil droplets was used to measure the velocity field via particle image velocimetry (PIV). The burned gas areas where oil droplets were evaporated represent the position of the flame. The velocity and burned gas measurements of the stab case were conducted using the setup (and same vector processing) as described in detail by Hanuschkin et al. (2021) and Dreher et al. (2021), who measured image pairs from − 360 to 0 CAD in increments of 5 CAD (1.8 kHz). For the unst case, oil droplets were illuminated by two laser sheets from pulsed Nd:YVO 4 cavities and were captured by a Phantom v2640 ultra high-speed camera equipped with a Sigma lens (APO Macro 180 mm F2.8 EX DG OS HSM, set to f/5.6), 532 nm bandpass filter, and a correction lens of focal length f = + 2000 (to counteract the astigmatism caused by the curved cylinder glass) at a sampling rate of 4.5 kHz ( − 180 to − 4 CAD , increments of 2 CAD). Burned gas regions were captured using the second frame of each image pair and indicate the position of the flame; therefore, the rest of this work will refer to extracted burned gas regions as the flame. The different laser and optical systems used coupled with the improved painting of the back of the cylinder glass resulted in fewer reflections in the PIV images of the unst case. Therefore, the resulting flame images of the unst case require less masking before analysis.
Vector processing of the unst case was conducted using the software DaVis 10.1.2 and employed a multi-pass cross-correlation technique (first 2 passes: 64 × 64 , 50% overlap; last 2 passes: 32 × 32 , 75% overlap; after each pass: peak ratio criterion of 1.3, universal outlier detection of 7 × 7 , and a vector group removal criterion of size 5). Further processing of velocity fields such as obtaining instantaneous line profiles or the mean vector field was conducted using MATLAB. Likewise, extraction of the flame for both cases from burned gas regions was conducted using an in-house MATLAB script whose procedure is as follows: After masking and dewarping (3rd order polynomial), raw images were grouped by their respective cycle then the intensity counts of each image were normalized by the maximum intensity achieved in the respective cycle. Afterwards, a sliding entropy filter (5 px × 5 px) was applied, followed by the division of each image by the first image without a flame. Next, the pixels whose intensities were less than 1 were set to 0, an erosion (disk radius of 4) removed small structures around all edges, a minimum area criterion of 200 px was applied, all holes were filled within each flame, and finally a dilation of the same size of the erosion was applied to restore the flame size.

Numerical Methods
In this study, LES of the optical engine was performed using a commercial 3D-CFD code, CONVERGE (version 3.0) (2009), which is based on a finite-volume approach. The computational domain is discretized with a Cartesian mesh with local mesh refinement as shown in Fig. 1. The largest cell size in the cylinder is 0.5 mm. Local refinement is applied in the region close to the cylinder walls, near the spark plug, and in the crevice, with a cell size of 0.25 mm. Similar resolutions were employed in previous engine LES studies (d'Adamo et al. 2019;Truffin et al. 2015;Zembi et al. 2019) and proposed by Zembi et al. (2019) as the best compromise to resolve more than 80% of turbulent kinetic energy with a computational cost not overly expensive. A coarser mesh is used in the intake and exhaust ports with the largest cell size of 4 mm. The computational domain is composed of about 2.5 million cells at top dead center (TDC) and about 7.5 million cells at bottom dead center (BDC). Simulation of a full engine cycle of 720 CAD took about 60 h on 480 cores.
The intake and exhaust boundaries of the computational domain are specified at the positions where experimental boundary conditions are available. In particular, measured 1 3 instantaneous static pressures and average temperatures are imposed at the boundaries. Since the variations of those boundary conditions over different cycles are negligible ( < 0.1 bar for pressure and < 1 K for the intake temperature), they are kept constant for multiple LES cycles. Blow-by is expected to be small in the investigated engine (Baum et al. 2014) and is not considered in the LES. The temperature on the walls of the combustion chamber and on the walls of the manifolds are specified as 60 • C and 36 • C in the LES, respectively. Dirichlet boundary conditions are applied for the crevice and cylinder walls, while the velocity wall model by Werner and Wengle (1993) and the heat transfer model by Amsden et al. (1989) are used for the walls of the manifolds where mesh resolution is coarser.
The compressible Favre-filtered governing equations of mass, momentum, energy, and species are solved using the finite-volume method and discretized by a blending second-order scheme in space. The Pressure Implicit with Splitting of Operator (PISO) method (Issa 1986) Fig. 1 Computational domain and mesh of LES is used for the pressure-velocity coupling. A second-order Crank-Nicolson (Crank and Nicolson 1947) time integration is employed. Details about the discretization schemes and the iterative algorithms can be found in the CONVERGE manual (2020). The one-equation eddy viscosity model (Yoshizawa and Horiuti 1985;Menon and Calhoon Jr 1996) is used for the LES sub-grid scale (SGS) closure, where a transport equation for the SGS kinetic energy, k sgs , is solved. The G-equation model (Pitsch 2002) is applied to describe the turbulent flame propagation, where the flame front is tracked using a levelset approach by solving a transport equation of the variable G (Peters 2009). The flame front is defined at G = 0 . � G < 0 and � G > 0 indicate the regions of unburned and burned mixtures, respectively. Chemical equilibrium (Pope 2003(Pope , 2004 is assumed for the burned mixtures. The transport equation for G reads where s t is the turbulent flame speed. − D t̃ + u s t can be understood as the propagation speed of the filtered flame front. The filtered flame front curvature ̃ and turbulent diffusivity D t are given by and where sgs is the SGS turbulent dissipation modeled by C (k sgs ) 3∕2 ∕Δ . Δ is the filter width. The SGS Schmidt number Sc and the modeling constants C and C are set as 0.78, 0.0845, and 0.05, respectively. The turbulent flame speed s t is given by Pitsch (2002) where s l is the laminar flame speed, is the dynamic viscosity, t is the turbulent viscosity, u ′ is the SGS velocity, and b 1 and b 3 are modeling constants. Eq. (1) is solved using the levelset approach. Outside the flame surface, the scalar is required to satisfy | ∇G |= 1 . The G -field is reinitialized as unburned in every cycle before ignition. A look-up table for the laminar flame speed s l has been generated with FlameMaster (Pitsch 1998) using a detailed iso-octane mechanism (Cai et al. 2019). The modeling constants b 1 and b 3 in the turbulent flame speed model are calibrated to match the mean in-cylinder pressure in the experiments. The necessity of such calibration might be attributed to the poor prediction of the SGS velocity and the underresolved flame wrinkling in the performed LES. Due to the different thermodynamic conditions of the stab and unst cases, b 1 and b 3 have different values for both cases. The ignition is realized by adding an energy source of 100 mJ within a sphere with a radius of 0.5 mm. Positive G is initialized in the region with a temperature higher than a threshold of 4000 K. The same ignition is used for both cases. Because of this simplified ignition model, which resulted in faster early kernel development, the ignition timing is calibrated to match the experimental mean pressure trace, especially for the first 10-20 CAD after ignition, as proposed by Esposito et al. (2020). An increase of b 1 or b 3 or earlier ignition timing leads to higher in-cylinder pressure. The calibrated model parameters are summarized in Table 2 for both cases. For the stab case, consecutive cycles are simulated. For the unst case, the skip-firing mode is realized in the LES by first simulating consecutive non-reactive cycles without ignition and then using these results as initial conditions for the reactive simulations to reduce the computational cost.

Comparison of LES and Experiments
This section presents a detailed comparison of the LES and the experiments. More than thirty cycles were simulated for both cases. The first two simulation cycles are not considered in the analysis, since they are influenced by the assumed initial condition (Truffin et al. 2015). In the experiments, about 200 and 1000 cycles are measured for cases stab and unst, respectively.

In-Cylinder Pressure
One of the most important aspects for the validation of engine simulations is the agreement in terms of in-cylinder pressure with experiments, which is reported for both stab and unst cases in Fig. 2. As displayed in the figure, the two cases have different spark timings. Note that the model parameters, b 1 , b 2 , and the ignition timing are calibrated to match the mean pressure trace in the experiments. For case stab, the standard deviation of the in-cylinder pressure is well reproduced, even though a slight deviation in the first half of the combustion stroke in terms of pressure slope is observable. This might be attributed to the limitation of the combustion model, which does not accurately account for the initial flame kernel development and the transition from laminar to turbulent flame propagation. However, further development of the model to account for this effect was not within the scope of this study. For case unst, larger differences between the LES and the experiments exist, especially for the variations of the maximal in-cylinder pressure, p max , which might be attributed to the different proportions of the misfire cycles. Only 1% of cycles in the experiment recorded misfire, which is defined as a cycle whose gross IMEP ( IMEP g ) was negative. In the simulations, one cycle out of thirty-two, about 3%, had misfire. Such difference might be attributed to the limited number of the simulated cycles. Since the pressure of the misfire cycle is significantly lower than the other cycles, the mean and the standard deviation of pressure can be strongly influenced by the proportion of the misfire cycles. If the single misfire cycle in the LES is assigned a weight of 1/3 so that the proportion of misfire cycles in the LES equals that of the experiments, a better agreement is observed in the cylinder pressure, as shown in the supplementary material ( Fig. S-1). However, the main trend does not change. Larger pressure variations are obtained in the LES compared to experiments. In addition, some differences are observed in the expansion stroke, which might be due to the poor prediction of incomplete combustion and heat loss. However, this will not affect the analysis of the CCV in this study, as the important range for the analysis of the CCV is the time before p max is reached. Therefore, it can be confirmed that the simulations are able to distinguish the stable and unstable engine operating conditions, which was also reported in Granet et al. (2012), Truffin et al. (2015), and d' Adamo et al. (2019). A quantitative comparison of the mean, the standard deviation, and the timing, CAD p max , of p max are provided in Table 3, where LES* stands for statistics evaluated by assigning a weight of 1/3 to the single misfire cycle in LES. The variability is assessed with the coefficient of variation (COV), which is defined as the ratio of the standard deviation to the mean value. As expectable from Fig. 2, while p max is well predicted for both cases, CAD p max is retarded for case unst. Figure 3 compares the distribution of p max for both cases. It is seen that the distribution is not symmetric. In the experimental data, a long tail in the low pressure side in the distribution is observed for both cases, especially for case unst. The same trend is also captured in the simulation. However, the extreme low and high pressure cycles are not predicted in the simulation, which might be attributed to the limited combustion model for the early flame kernel development and the limited number of the simulated cycles.

Flow Fields
In this section, the velocity fields during the intake and compression strokes are evaluated. Figure 4 compares the 2D velocities projected onto the tumble plane passing through the cylinder axis and their standard deviations at − 270 CAD in the intake stroke for case stab. Good agreement can be observed for the overall flow structure and magnitude of the velocity fields. The high intake velocity and the stagnation regions are well predicted as can be observed in Fig. 4a. For further comparison, the two velocity components along the white dashed horizontal lines at different vertical locations are compared in Fig. 4b and c. Very good agreement is achieved for the averaged velocities (designated by the thick lines in Fig. 4b). In Fig. 4b, the range of the vertical axis is chosen to show the velocity fluctuations of individual cycles. Comparisons of the averaged velocities with a smaller range of the vertical axis are shown in the supplementary material (Figs. S-3 and S-4). The fluctuations of the velocity components in the simulation also agree fairly well with the experiments in terms of absolute values and trends, even though stronger fluctuations are observed for the LES. Furthermore, the standard deviation of the velocity components for the experiment is relatively smooth along the x-axis, while it fluctuates for the LES. These are probably due to the large difference in the sample sizes between the sets in LES and experiment. Similar results are obtained for case unst, which are shown in the supplementary material ( Fig. S-2). Figure 5 compares the velocity in the compression stroke close to spark timing for both cases. For case stab, the averaged velocities ( Fig. 5a and b) and the standard deviations (Fig. 5c) are fairly well predicted, albeit the magnitudes of the mean values are slightly higher in the LES. In addition, a vortex is observed in the simulation at the top right corner near the exhaust valves as shown in Fig. 5a. However, in the experiment, this vortex is outside the field of view. Such a difference in the location of the vortex explains the discrepancies in w for x > 20 mm in Fig. 5b. Similar results are observed for the velocity fields of the case unst. The global flow structure is captured by the LES. However, also in this case,  the magnitudes of the averaged velocities are slightly higher for the LES and the location of the vortex center slightly differs (Fig. 5d), namely it is closer to the spark plug and thus a large difference in w for x > 10 mm is observed (Fig. 5e).

Spatial Flame Distribution
After the comparison of the flow fields, this section analyzes the flame propagation for both cases. A quantitative comparison of the flame propagation can be conducted based on the flame probability, P f , which describes the local probability of the occurrence of burned gas at a given CAD, defined as  Figure 6 compares the flame probability of LES and experiments for both cases considering all cycles of each case.
(5) Since the two cases have different spark timings, combustion duration, and thus frequencies of PIV measurement, flame images are shown for different CADs for the two cases. It should be noted that the experimental flame images for the stab case had to first be heavily masked due to strong laser reflections on the spark plug; for the unst case, the experiment was conducted with a different laser and camera imaging setup, as described in Sect. 2, and reflections were better mitigated. Good agreement between numerical and experimental flame propagation is observed for both stab and unst cases. For case stab, the spatial flame distribution is concentrated with high peak P f values and the flames tend to always first propagate to the left of the spark plug, which is observed in both the LES and experiments. Flames in LES tend to occur further to the top left compared to the experiments, which is consistent with the slightly overpredicted velocities from the LES before ignition as shown in Fig. 5. For case unst, the flame is more distributed throughout the field of view in this plane with lower overall peak in P f values. Compared to the stab case, the flame on the left side of the spark plug is more concentrated and the flame on the right of the spark plug is more distributed due to the flame interactions with the flow, which has a different state near the spark plug at ignition (Fig. 5). Furthermore, similar spatial flame distributions are seen in both LES and experiments. The typical flame propagation in a 3D projection obtained in LES for a single cycle of both cases is illustrated in Figs. 7 and 8. For both cases, the flames show asymmetry with respect to the tumble plane until − 10 CAD. This is visible especially for the unst case, where the flame is strongly stretched by the vortex on the exhaust side resulting in multiple pockets of burned gases.

Conditionally Averaged Flow Fields and Spatial Flame Distributions
In addition to the comparison of statistics from all cycles, it is interesting to compare statistics of the flow fields and the spatial flame distributions conditioned to the highest and lowest pressure cycles. Figure 9 compares the velocity fields before ignition Fig. 6 Flame probability P f in the tumble plane (y = 0) considering all cycles. The empty region in the experiments is due to reflection conditioned on the cycles with the 20% highest and lowest p max for both cases. For case stab, the same trend as in the experiments is predicted in the LES: larger velocities in the tumble plane correlate to high p max . For case unst, the correlation of p max and the location of the vortex, which is characterized by the low velocity magnitude at the top right corner, is different between LES and experiment. While the vortex center is larger and thus closer to the spark plug for high pressure cycles in the experiment, it is closer to the spark plug for low pressure cycles in LES. In addition, for the unst case, the experimental flows in the horizontal direction across the spark plug gap are greater, which accounts for the downward shift of the vortex center when comparing high pressure cycles to low pressure cycles. This difference in conditionally averaged flows between the experiments and LES is mainly caused by discrepancies in the combustion process, which influences the selection of high and low pressure cycles. It can be attributed to the simplified ignition model in the simulation, which is an energy source at the spark plug. While the spark channel in reality is convected by the flow and thus the flame kernels are initiated at different locations, the process of flame initiation occurs at the same spot in the LES and possible convection of the flame by the flow can happen only at a later stage. A more sophisticated ignition model (Dahms et al. 2009) will be employed in future studies. Figure 10 shows the spatial flame distribution conditioned on the high and low pressure cycles. For case stab, similar results are observed between the simulation and experiment. High pressure cycles exhibit larger regions of burned gas and the flame is mainly on the left side of the spark plug for both high and low pressure cycles.   Fig. 9 Phase-averaged 2D velocity fields at − 30 CAD for case stab (a) and − 50 CAD for case unst (b) conditioned on the 20% highest and lowest p max cycles. The empty region in the experiments is due to reflection and the evaporation of the PIV particles However, for the unst case, different behavior is observed in the simulations compared to the experiments. In the experiments, flames for high pressure cycles are more likely to be located on the right side of the spark plug until − 18 CAD, which is consistent with the greater horizontal flow across the spark plug and the closer location of the vortex center to the spark plug as shown in Fig. 9. The flames for low pressure cycles tend to spread to the left side of the spark plug. In the simulations, flames tend to spread to the left of the spark plug for both high and low pressure cycles. This might be attributed to the simplified ignition and flame-wall-interaction models and the discrepancy of the predicted flow fields in the vicinity of the spark plug electrodes in the simulation.

Fig. 10
Flame probability P f in the tumble plane (y = 0) for case stab (a) and unst (b) conditioned on the 20% highest and lowest p max cycles. The empty region in the experiments is due to reflection

Analysis of CCV
After the detailed comparison of LES and experiments shown in the previous section, the sources of the CCV in the simulations are analyzed here. Despite the differences between the LES and experiments, the analysis of the CCV in the LES may shed new light on the causal chain of engine CCV with the available three-dimensional data.

Influencing Variables
The influencing variables are introduced in this section, which are extracted at the spark timing on different scales, including the quantities evaluated in the spherical region ( r = 3 mm) centered between the spark plug electrodes, the parameters for the threedimensional coherent vortex structure evaluated locally considering a spherical neighborhood ( r = 12 mm), and the large-scale parameters considering the whole cylinder ( r = 43 mm).

Small Region Close to Spark Plug
Ignition and the very early combustion phase are crucial for CCV (Schiffmann et al. 2018). Therefore, variations of the thermodynamic and fluid-dynamic conditions in the region close to the spark plug are expected to be one of the major sources of CCV. In this study, the local fields close to the spark plug in a spherical region with a radius of 3 mm centered between the spark plug electrodes (Fig. 11) are considered. The influence of the following variables will be assessed in Sect. 5.2: spatially averaged temperature T sp , SGS kinetic energy k sgs sp , pressure p sp , resolved velocity components u sp , v sp , and w sp , resolved velocity magnitude |u| sp , and fuel mass fraction Y f,sp in the region close to spark plug.

Coherent Vortex Structure
Another group of influencing variables can be defined based on the coherent vortex structure. In this study, the vortex structure is detected by the function, 3p , defined as (Buhl et al. 2017) where ̂ o is given by ⟨ ⟩ V ( ) refers to the average velocity in the volume V around point . In this study, V is defined as a sphere with a radius of r = 12 mm centered at . p and ̂ p are the position and velocity vector projected onto the plane normal to the axis of rotation given by and (6)

3
The position vector is defined as = o − . The direction of the axis of rotation is given by where the function 3 ( ) reads (Gohlke et al. 2008) | 3p | = 0 means there is no rotation, while the center of an axisymmetric and uncurved vortex structure has unity | 3p | . Examples of the vortex structure of a single cycle close to spark timing are shown in Fig. 12 for both cases. Different vortex structures can be identified. For the unst case, since the ignition time is earlier, the tumble motion can be identified by the coherent vortex core. Conversely, since the ignition for the stab case occurs later, the organized tumble motion breaks into a vortex pair rotating in the x-direction.
To quantify such vortex structures, the locations of the vortex centers on the planes across the spark plug are analyzed. For the stab case, the two vortex centers on the plane, x = x sp , are considered since the vortex pair is rotating in the x-direction. The vortex center on the plane, y = y sp , is analyzed for the unst case because the rotational motion is roughly around the y-direction. The vortex centers are determined by the maximal | 3p | on these planes. The distances between spark plug and the vortex centers, d y<0 , d y>0 , and d, illustrated in Fig. 12, are considered. The local tumble ratio, LR, which describes the strength of the local rotational flow around the vortex center, is also evaluated as where c is the vortex center and ref is the angular speed of the engine.

Large-Scale Variables
The large-scale variables include the averaged thermodynamic parameters of the in-cylinder mixture and the large-scale flow motion considering the whole cylinder ( r = 43 mm). Fig. 11 Illustration of the spherical region close to spark plug with a radius of 3 mm For both cases, variations of the in-cylinder pressure and the averaged temperature at spark timing are negligible (about 0.001 MPa and 2 K), and thus not analyzed in this study. Since the amount of EGR is well controlled and homogeneous for the unst case using skip-fire operation, the influence of the mass fraction of EGR, Y EGR , is only considered for the stab case. To quantify the large-scale flow motion, the tumble ratio (TR), swirl ratio (SR), and cross-tumble ratio (CR) are considered, which are defined by and where V cyl. is the in-cylinder volume. xz , xy , and yz are position vectors referring to the center of mass of the in-cylinder gases projected onto the xz, xy, and yz planes, respectively. xz , xy , and yz are the velocity projected onto the xz, xy, and yz planes, respectively. The positive and negative values of these large-scale flow motion parameters may cancel each other when the flow field exhibits symmetry. As shown in Fig. 12, the flow fields are roughly symmetric about the plane, y = 0 . Therefore, the positive and negative contributions are considered separately. Note that these values not only quantify the strength of the large-scale rotational motion but are also influenced by the location of the coherent vortex center relative to the center of mass. Fig. 12 Vortex core and streamlines around the vortex core of a single cycle at the spark timing for both cases 1 3

CCV Correlations
To systematically analyze the CCV correlations in LES, CCV are distinguished between the early flame kernel development and flame propagation. The former is quantified by the burned fuel mass at 3 CAD after spark timing, m fb,ST+3CAD . This short time after ignition is chosen to investigate the influence of the ignition and the very early flame kernel development. Even though the flame kernels are quite small at such a short time after ignition, their CCV can already be identified, which significantly influences the maximal incylinder pressure, p max , as will be shown in the following. Figure 13 shows the correlation between m fb,ST+3CAD and p max . The burned fuel mass is normalized by the averaged trapped fuel mass before combustion, ⟨ m f,tot ⟩ . For the unst case, the cycles with low m fb,ST+3CAD are highlighted in red and the misfire cycle is displayed with the rectangle, which will be discussed in Sect. 5.2.2. To quantify the correlation, the Bravais-Pearson correlation coefficients are displayed in Figs. 13, 14, 15 and 16, which are defined as where x and y are the standard deviations of variables x and y, respectively. xy is the covariance of x and y. |R| = 1 or |R| = 0 means that there is a perfect linear correlation or no correlation between the variables, respectively. The significance level is also illustrated in Figs. 13, 14, 15 and 16 by the p value for the H 0 -hypothesis test with t-distribution of the test statistic (Fisher 1992). A correlation is considered significant if p < 0.05.
A strong correlation is observed for both cases, indicating that the early flame kernel development plays a significant role in the CCV of p max . The relatively large scatter in both plots of Fig. 13 reveals that the combustion process varies also after the early flame kernel development. It is worth noting that the misfire cycle is not an extreme case in terms of the early flame kernel development (not the lowest m fb,ST+3CAD in Fig. 13b), which agrees with the experimental study by Peterson et al. (2011). They showed that misfiring is not a consequence of the failed ignition. CCV of the flame propagation will be analyzed based on the variations of p max excluding the influence of m fb,ST+3CAD , which is quantified by Fig. 13 Correlation between the burned fuel mass at 3 CAD after spark timing and the maximal in-cylinder pressure for case stab (a) and unst (b). For the unst case, the cycles with low m fb,ST+3CAD are highlighted in red and the misfire cycle is displayed with the rectangle where p f it max (m fb,ST+3CAD ) is the fitting line as shown in Fig. 13. In other words, Δp max is the difference between the observed pressure to the fitted line in Fig. 13. In the following, the correlations of these two quantities with the influencing variables introduced in Sect. 5.1 will be analyzed.

Early Flame Kernel Development
In this section, correlations of m fb,ST+3CAD and the influencing variables are analyzed. Figure 14 shows the correlation coefficients and the p values of the respective variables. As can be observed in Fig. 14, the importance of the influencing variables is different for both cases. For the stab case, a significant negative correlation ( p < 0.05 ) is found for the coherent vortex structure (12 mm), namely the minimum distance between the vortex center and the spark plug, d min = min d y<0 , d y>0 . This indicates that the early flame kernel at 3 CAD after the spark timing already starts to interact with the coherent vortex structure illustrated in Fig. 12. Consequently, a smaller distance leads to faster early flame kernel growth. However, the strength of the local rotational motion (LR) does not show a significant correlation. Additionally, no significant correlations can be identified either for the variables extracted in the small region close to spark plug ( r = 3 mm). For the unst case, a weak correlation is found for the average velocity in the y-direction close to spark plug, v sp , indicated by R ≈ 0.3 and p ≈ 0.08 . A larger v sp means that the ignited kernel will be pushed towards the spark plug (Fig. 11), which leads to higher heat loss and thus slow combustion. For this case, the flame kernels are smaller compared to the stab case and have not started to interact with the vortex at 3 CAD after spark timing. Thus, no significant correlations are observed for the coherent vortex structure ( r = 12 mm). For both cases, the influence of the large-scale variables considering the whole cylinder ( r = 43 mm) on the early flame kernel development

Flame Propagation
This section investigates sources for the CCV of the flame propagation, namely, Δp max . The correlation coefficients and the p values between Δp max and the influencing variables are shown in Fig. 15. For the stab case, the flame propagation is mainly influenced by the large-scale variables where correlations are higher and p values are lower. Significant correlations are found for the internal EGR concentration Y EGR and the positive swirl ratio SR+ . For the unst case, since the flame volume at 3 CAD after ignition is small, Δp max still correlates to the variables extracted from the region close to the spark plug and is strongly influenced by | | sp with |R| > 0.4 and p < 0.01 . Significant correlations between Δp max and k sgs sp or d are also observed, which is mainly because of their correlation with | | sp as shown in Fig. 16.
The cycles with small early flame kernels (low m fb,ST+3CAD ) are highlighted in red and the misfire cycle is displayed with a rectangle. Since the flame speed is small for the unst case, the flow close to the spark plug, especially the rotational motion, can stretch the flame and reduce the propagation speed. Therefore, larger | | sp , smaller d, or stronger LR leads to slower flame growth and ultimately lower Δp max . It can be observed that the misfire cycle has the highest | | sp , the smallest d, and the strongest LR among the cycles with small early flame kernels (red points in Fig. 16) indicating that misfire is related to the strong stretch of the flame kernel by the flow close to the spark plug and the coherent vortex structure. It is worth noting that the applied levelset model always assumes a positive turbulent flame speed (Eq. 4) and thus cannot predict the local flame extinction. The misfire cycle is obtained in the simulation because the flame kernel is stretched by the flow, generating flame segments with large curvatures, which reduces the propagation speed of the levelset front to even negative values via the first term on the right hand side of Eq. 1. Furthermore, for the large-scale variables, a significant correlation between Δp max and the negative cross-tumble ratio, CR− , is seen. CR− is also influenced by the strength and location of the coherent vortex and thus correlates to | | sp as shown in Fig. 16d.

Conclusions
This paper presents a combined numerical and experimental investigation of an optically accessible SI engine with port fuel injection. Two operating conditions are considered: a stable operating condition with continuous firing of a stoichiometric fuel-air mixture and an unstable operating condition in a skip-firing mode with homogeneous external EGR.
The LES is carried out with well-established turbulence, combustion, and ignition models. A detailed comparison between LES and experiments is performed. Even with some difference in the expansion stroke of the unst case, a good agreement of the incylinder pressure, especially the maximal in-cylinder pressure and its cycle-to-cycle variation, is obtained. It is also confirmed that LES is able to distinguish the stable and unstable operating conditions in terms of the flow field and flame propagation. Good agreements for the general flow structure, the phase averaged velocities, and their fluctuations during the intake and compression strokes are achieved. However, small discrepancies are observed for the location of the vortex center in the compression stroke. The spatial flame distribution is also well reproduced in the LES when considering all cycles of each case. For the stab case, the same trends as the experiments are also found in the LES for the flow fields and the spatial flame distribution conditioned on the high and low pressure cycles. Conversely, discrepancies are observed for the conditionally averaged quantities of the unst case, an observation that is likely attributed to the lack of a spark plasma and subsequent flame kernel growth in the ignition model.
A systematic quantitative correlation analysis between the influencing variables extracted from different scales of the 3D LES and the CCV of the early flame kernel development and the consequential flame propagation is conducted. The influencing variables include the spatially averaged mean temperature and pressure and the flow features in the spherical region ( r = 3 mm) centered between the spark plug electrodes, parameters for the three-dimensional coherent vortex structure ( r = 12 mm), and the large-scale parameters that consider the whole cylinder ( r = 43 mm). The most relevant influencing variables are different for different operating conditions and combustion phases. However, for both cases, the location of the coherent vortex center plays a significant role in the combustion process, which is also confirmed by the conditionally averaged flow fields in the measurement of the unst case. For the unst case, the flow direction in the region close to the spark plug is found to influence the early flame kernel development. A slower growth of the early flame kernel is observed if it is pushed towards the spark plug. For the stab case, the internal EGR is found to significantly correlate to the flame propagation. The large-scale flow features are also found to influence the flame propagation for both cases. Such main findings in terms of sensitivities regarding CCV are in close agreement with the experimental study by Welch et al. (2022). Finally, the flame kernel interactions with the flow close to the spark plug and the coherent vortex structure are found to be the cause for the misfire cycle in LES of the unst case, which indicates that proper control of the vortex location might be promising to avoid misfire.
Acknowledgements The authors gratefully acknowledge the funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Research Unit FOR2687. The authors gratefully acknowledge the computing time granted by the NHR4CES Resource Allocation Board and provided on the supercomputer CLAIX at RWTH Aachen University as part of the NHR4CES infrastructure. The simulations for the unst case were conducted with computing resources under the project p0020108. Simulations of the stab case were performed with computing resources granted by RWTH Aachen University under project rwth0501. Convergent Science provided CONVERGE licenses and technical support for this work.

Author Contributions
All authors contributed to the study conception and design. BB and AD contributed to the design and the conception of the experiments. Experiments were performed by CW and BB. Simulations of the stable operation condition were conducted by HC and HE. Simulations of the unstable operation condition were conducted by HC and SC. MD and HP contributed to the conception of the data analysis. Analysis was conducted by HC, HE, SC, and MD. The first draft of the manuscript was written by HC and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL.

Declarations
Competing interests A. Dreizler is on the editorial board of Flow, Turbulence and Combustion and B. Böhm is guest editor of the special issue: Cyclical Variations in Internal Combustion Engines; however, neither of them handled the review nor editorial process of this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.