Flow and thermal field effects on cycle-to-cycle variation of combustion

• Cycle-to-cycle variation is investigated in a spark-ignited lean gas engine. • Flow ﬁ eld was shown to have a higher relative contribution to the cyclic variations. • Early variations near spark region are the major origins of the cyclic variations. • Burning rate and convection velocity ahead of ﬂ ame described the ﬂ ame asymmetry. • Increasing the spark size signi ﬁ cantly reduced the cyclic variations.


Introduction
Cycle-to-cycle variation (CCV) in internal combustion (IC) engines refers to the non-repeatability between different combustion cycles.During the past decades, numerous research efforts have been made to better understand and explain the phenomenological origin of CCVs.CCVs are known to lead to various undesirable effects such as reduced efficiency [1] and increased emissions [2].Various engine combustion concepts including homogeneous charge compression ignition (HCCI) [3], reactivity controlled compression ignition (RCCI) [4], spark ignition (SI) lean premixed combustion [5], stratified combustion [6], and engine downsizing [7] have been proposed to provide cleaner and more efficient combustion process [8].However, the operational range of such concepts may be significantly limited by CCVs [9].
Although CCVs are present in both compression ignition and SI engines, they are more commonly associated with SI concepts [10].Particularly, in lean burn SI engines, large cyclic variability has been shown to lead to considerably increased levels of emissions [10], increased fuel consumption of as high as 6% and decreased power output of even 10% [11].Yet, when compared to stoichiometric burn SI combustion, lean burn SI combustion has been noted to have various benefits such as increased thermal efficiency [12] and decreased NOx emissions [13].However, the operational range of lean combustion is strictly limited by CCVs [14,15] as the concept is known to be relatively sensitive to CCVs.An improved and detailed understanding of CCV sources would be a crucial step towards better understanding and, eventually reducing CCVs.
Experimental studies have provided invaluable insights into the CCV problem.For example, previous experimental studies have shown the connection of CCVs with exhaust gas recirculation [16], global fuelair ratio [17], in-cylinder turbulence level [12], engine speed [18], inlet air temperature [19], spark timing [20,21] and knock tendency [22].However, experimental engines are often prone to uncertainties in terms of exact boundary and operating conditions.In fact, as enabled by optical engine measurement, it is nowadays possible to measure various details of the in-cylinder flow field including planar images of flow field, temperature and mixture stratification [23].Attempts have also been made to optically access the full three-dimensional (3D) fields in a spark-ignition engine using temporally resolved multi-planar laser diagnostics [24].However, there are limitations regarding optical access in certain set-ups which restrict diagnostic possibilities and gaining full access to the 3D flow quantities within the turbulent combustion process.
Limitations in studying CCV sources are present in computational fluid dynamics (CFD) approaches as well.In the Reynolds Averaged Navier-Stokes (RANS) approach, all turbulent scales are modeled and consequently RANS approaches provide only information about a typical, average cycle.Hence, detailed phenomena such as flame-turbulence interaction cannot be studied in detail with the RANS approach.Such details are expected to be important when studying the unsteady characteristics and formation of CCVs [25,26].
On the other hand, large-eddy simulation (LES) has been shown to be able to quantitatively reproduce CCV through various studies e.g.[27][28][29].Most of the reacting LES CCV studies mainly concentrate on reproducing CCV levels in various conditions in order to establish the predictive capability of the numerical approach e.g.[30][31][32][33][34]. Thickened [30] and coherent [31] flame models have been employed with LES, successfully reproducing CCV in single-cylinder engine configurations.Curto-Risso et al. [32] investigated the effect of fuel-air ratio on CCV while reproducing heat release rate data from experiments.Later, Granet et al. [33] performed both LES and experiments of CCV in stable and unstable operating points of an SI engine.LES reproduced the flame behavior in both cases in terms of combustion progression and flame shape.In fact, quantitative prediction of CCV amplitude was demonstrated as well.In a more recent work, Ameen et al. [34] provided a parallel perturbation model to dissociate the long time-scale CCV problem into several shorter time-scale problems.The provided model was shown to be capable of predicting quantitative trends of CCV.We note that the studies above demonstrated how LES produces CCVs in certain engine configurations.Yet, the studies provided only limited insight to the interactive physical and chemical processes behind CCVs.Such phenomena are further investigated in the present paper in a more simple engine configuration.Some of the reacting LES studies concentrated on the causes of CCVs in SI engines [27][28][29][35][36][37].Vermorel et al. [27] used LES in SI fourvalve single cylinder engine fueled with a homogeneous propane-air mixture.It was suggested that variation in the coherent tumble motion, produced during the intake stroke, was the major triggering factor for CCVs.In the studied case, the effects of local or overall mixture variations were reported to be insignificant.Enaux et al. [28] carried out an LES study of CCVs in a propane-fueled SI engine.It was suggested that CCVs are essentially due to velocity fluctuations around the spark plug, which induce variations of the early flame kernel growth and of the overall combustion duration.A similar observation regarding the significance of the convection flow effect on the spark was recently reported by Truffin et al. [36].It was suggested that causes of CCV depend on the type of engine and its operation conditions.In a more recent work, Kodavasal et al. [29] used machine learning techniques to understand CCV causes.They suggested that machine learning techniques can implicitly learn complex relationships between different parameters defining the eventual outcome of each cycle.
Furthermore, Fontanesi et al. [35] performed a multi-cycle LES of a V-8 engine.They reported that inhomogeneity and variation in fuel distribution, turbulent energy, and velocity magnitude in the spark gap from cycle to cycle were major contributors to CCV.Koch et al. [37] performed an LES of multi-cycle SI engine using the G-equation, and Damköhler turbulent flame speed (s t ) closure.They reported that the cyclic variations are mainly related to the fluctuations in the sub-grid scale (SGS) kinetic energy rather than the thermodynamics conditions.However, it is noted that in such closure terms, the s t is usually tuned to represent the experimentally observed burning rate.Thereby, we note that any inaccuracy in the modeling of laminar flame speed may be compensated by tuning the s t for matching experimentally observed global combustion trends.A similar observation was recently reported by Burke et al. [38] who carried out a comparative study on different s t correlations.With relevance to s t correlation modeling, commonly used laminar flame speed (s L ) correlations, such as Gulder [39] and Metghalchi correlation [40], are known to be inaccurate at high pressure and temperature conditions.Thus, as part of this study, a detailed chemistry based model is developed for more accurate description of chemistry in such conditions.
Although lean combustion is known as a useful concept to reduce fuel consumption and emissions in SI engines, it has been noted to be prone to CCVs because the reduced flame speed in lean mixtures may not ensure a robust flame propagation [41,42].Such aspects make lean combustion an important topic of study especially in the CCV context.However, reactive LES studies only rarely concentrate on understanding the sources of CCVs in lean SI engines.Additionally, the impact of temperature on flame speed diminishes when the mixture is close to stoichiometric condition [43].This potentially contributes to the relatively low sensitivity of combustion to temperature stratification reported previously in the literature e.g.[37,28].
As a summary of the literature survey, we identify the following main research gaps on SI lean burn combustion.First, the role of turbulence and thermal field on local and global flame front propagation is presently not fully understood in lean conditions.Second, the role of initial spark kernel development and its connection to CCVs is presently unclear.In the present study a simple, spark-plug free configuration is investigated in order to assess the role of turbulence and thermal stratification on flame propagation.Third, it seems that much focus should be put on laminar flame speed modeling in high temperature and pressure conditions.This is important to better capture the effect of thermodynamical conditions on flame speed, especially for lean mixtures.Fourth, with relevance to thermal stratification, accurate wall heat transfer modeling has not been demonstrated previously in the context of CCV modeling.
The current numerical study concentrates on the effect of thermal and flow fields on CCV in a SI IC engine-like configuration with particular consideration of lean combustion conditions.The present computational setup can be considered to be an extension of a widely studied engine-like set-up [44,45] extended to a simple intakecompression configuration [46][47][48] with the expansion phase additionally included in the present work.The main objectives of the present study are the following: (1) Present an improved laminar flame speed correlation for lean mixtures based on detailed chemistry simulations for high temperature/pressure conditions, (2) Validate the present set-up in a non-reacting flow configuration with respect to (a) phase-averaged velocity and RMS velocity profiles gathered from 17 cycles during intake processes, and (b) net wall heat transfer and thermal stratification during compression, (3) Numerically investigate the occurrence of combustion CCV in the studied set-up, (4) Identify the sources of CCV in the present configuration and quantification of the temperature and flow field effects on CCV by computing local variables around the flame, (5) Investigate the asymmetric behavior of the flame in the system, (6) Provide numerical evidence on the sensitivity of the initial, modeled spark kernel size on the occurrence of CCVs.
The paper is organized in four sections.In Section 2, the numerical framework, geometry and case set-up are described.In Section 3, first the non-reacting and then the reacting results are presented.Moreover, the local flame analysis, asymmetric behavior of the combustion and kernel size effect are discussed in this part.In Section 4, applicability of present results in real applications are discussed.In Section 5, a summary of the findings is provided.

Turbulence modeling: hybrid LES/RANS with wall treatment
The present study employs a zonal hybrid LES/RANS methodology with wall treatment (HLR-WT), which combines existing zonal hybrid [49] and wall modeling [50] approaches.The main idea is to use LES with an explicit SGS model in the core flow, while a low-Reynolds number RANS model ( −∊ k in the present setup) is adapted for the nearwall regions.Wall and core regions are separated with a fixed interface where continuity is imposed for modeled viscosity.Furthermore, to avoid substantial wall-normal grid requirements due to the very small viscous length scales (thin boundary layer), a one-dimensional wall model is used.The wall treatment developed by Nuutinen et al. [50] in the RANS context is specifically designed for engine-like boundary layers and engine wall heat transfer and attempts to account for nonequilibrium effects found in engine boundary layers.A recent measurement-based model analysis by Ma et al. [51] indicated that nonequilibrium models can improve the momentum and thermal boundary layer predictions compared to the empirical functions and equilibrium models.A thorough elaboration on the method and its implementation can be found in [52].Additionally, to close the filtered Navier-Stokes equations, the σ -model [53], recently advocated in engine-like flows [54], was used.

Premixed combustion model: G-equation
The G-equation model proposed by Williams [55] uses a level-set method to describe the evolution of the flame front as an interface between the unburned and burned gases.A non-reacting passive scalar G is introduced where the isosurface = G G 0 divides the domain into unburned and burned regions with < G G 0 and > G G 0 , respectively.The instantaneous and local G-equation can be derived by considering the instantaneous flame surface.A filtered equation is presented by Pitsch and Duchamp de Lageneste [56] for the LES context: where t ρ u , , , -and ∼ denote time, density, velocity vector in unburnt mixture, spatial filtering and spatial Favre filtering, respectively.An additional transport equation is solved for the variance G′ which is used to determine the turbulent flame brush thickness.In Eq. (1), s t refers to the turbulent flame speed.The numerous s t correlations available in the literature [57][58][59][60][61][62][63] are rarely validated for pressures above atmospheric conditions.In the present implementation, the empirical correlation from Muppala et al. [64] has been adapted: (3) while ′ u is the modeled fluctuation velocity related to the modeled turbulent kinetic energy (k mod ) by the assumption of isotropic turbulence: s L in Eq. ( 2) is representing the laminar flame propagation speed which generally is a function of composition and thermodynamic state of the mixture.In the context of the G-equation, s L has been usually predicted by adopting an empirical correlation such as the one introduced by Gulder [39].However, lack of laminar flame speed experiments at engine-like high pressure and temperature conditions makes it difficult to assess the correct closure for the G-equation.Recently, there have been studies concerning improving the prediction capability of s L for methane fuel such as [66,67], nonetheless the valid operating ranges of such correlations failed to cover the entire present case conditions.Therefore, s L has been estimated here from detailed chemical ki- netics simulations in one-dimensional laminar configurations.Detailed chemical kinetics can provide a consistent basis for studying the response of a laminar flame on changes in temperature (T), fuel-air equivalence ratio (ϕ) and pressure (P).However, at high T and P conditions, the laminar flame speed problem may become ill-posed and the deflagration wave can include an auto-ignition phenomenon, changing the flame speed [68].Recently, Krisman et al. [69] carried out DNS on laminar flame propagation at engine-like conditions and demonstrated that despite the ill-posed problem configuration, one-dimensional flame propagation solution from detailed chemistry could be used as a good approximation with respect to the DNS reference.
As the first objective of this study and prior to three dimensional (3D) computations, a new correlation function for laminar flame speed has been structured.The correlation structure is based on the general knowledge of laminar flame speed behavior with each parameter in the correlation.For instance, it is well established that laminar flame speed has a parabolic behavior with the equivalence ratio.In addition, the shape and applied algebraic relations are based on the known flame physics [68] and previously introduced correlations such as [39,67].Also, various versions have been tested to ensure a short and yet concise reproduction of detailed chemistry solutions in the desired conditions.Eq. ( 5) presents the function format with primary variables T P , and ϕ, and where = T T/1000 f and = P P/10 f 6 are normalization constants.Open-source chemical kinetics library Cantera [70] was used to solve a laminar free-flame problem with the GRI-3.0methane mechanism [71] at conditions spanned by T = 400-950 K, P = 1-50 bar and ϕ = 0.45-1.1.It is worth mentioning that the choice of chemical ki- netics could have a considerable impact on the laminar flame speed solutions [72].
A multi-variable least-squares fitting procedure for the s L data was carried out to estimate the constants (c 1 -c 9 ) in Eq. ( 5).Due to the nonlinear behavior of chemical kinetics, we found it beneficial to divide the T-P-ϕ domain into 6 regimes (Table 1), each having unique fit constants (Table A .5).The success of the fit was measured with respect to L2 error which was below 10% and we noted a maximum absolute error of 26% in the T-P-ϕ domain.The maximum error is related to regions at the domain boundaries and not on the mean conditions encountered mostly in the present LES (Table A.6). Fig. 1 represents a typical comparison between a fit-function and detailed data.The alternative approach of s L data tabulation was neglected here for the sake of use of the fit function which is easier to implement and employ in future studies., , , , , , , ], [1, log( ), log( ), , , , log( ), , ], /1000, /10 , , (5)

Ignition treatment
In the current study, ignition is treated numerically by initializing a spherical flame kernel with a predefined diameter of 0.2 mm during one time step at the time of the ignition.After initialization, the G field undergoes a preliminary process to enforce the condition that the function G is a signed distance function.Such ignition treatments have been employed in similar studies such as [33,37,29], while more elaborate ignition treatments have been employed in [73,27,74,75], and specifically for G-equation in [76].It is worth mentioning that the transition phase from laminar to turbulent flame speed has been neglected in this implementation, so from the beginning, the flame is propagating with turbulent flame speed similar to the approach adopted by Koch et al. [37].The spark timing is at 330 CAD for all the studied reacting cases.

Geometry and charge preparation
The present computational set-up is based on a three stage process to provide the initial charge for non-reacting cases (according to DNS by Schmitt et al. [48] utilized for validation purposes), and for reacting cases (to generate conditions relevant to lean-burn gas engines) as illustrated in Fig. 2 and elaborated in what follows.
Stage I (cold flow simulation): Stage I is an isothermal incompressible flow case corresponding to the experimental study of Morse et al. [44].In this stage the valve is constantly open while the piston is moving between top dead center (TDC) and bottom dead center (BDC) positions.This stage runs for 17 cycles and provides 17 naturally generated flow fields which will be used for later stages.The dynamic interactions of the incoming hollow jet flowing through the valve with the flow at TDC lead to CCVs in such reciprocating engines [45].
Stage II (intake simulation): A single naturally generated flow field from Stage I is used in Stage II, and the fuel composition and temperature are uniformly set in the intake channel part (Fig. 2, marker 1).Furthermore, a fully burnt mixture of air and residual gases with T = 900 K is set in the initial field below the valve level (see Fig. 2, marker 2).
Stage III (compression/combustion): In this stage, we have a fuel-air mixture in a closed valve cylinder.The flow field in Stage III is provided by the intake process in Stage II which was originally initialized with Stage I flow field.Additionally, the fuel-air composition, initial temperature and boundary conditions were set already at Stage II.The mixture is compressed, ignited and expanded at this stage.
It is worth noting that cycles 1-3 chosen from the Stage I and used in the reacting configuration are not consecutive.The three Stage II initializations (Stage I TDC flow fields) have been selected on the basis of the most differing azimuthally averaged flow fields in the Stage I intake stroke following [45].These differences are not assumed to directly imply combustion variations, rather, the procedure represents a simple way of selecting cycles with meaningful flow field differences.
For the non-reacting simulations, DNS studies such as [46][47][48] have been used as the reference data to validate the implemented models and therefore have been precisely followed when setting the case parameters.However, in order to provide a more interesting case regarding CCV in reacting cases, certain parameters such as equivalence ratio, fuel and wall temperature have been modified.Table 2 describes the Stage III geometry specifications.In the studied engine-like geometry, no connecting rod is available and the piston position is calculated from a sinusoidal function following [46][47][48].Also, all the studied cases are summarized in Tables 3 and 4.

Mesh
In this section, the computational meshes employed in the simulations are discussed.For the non-reacting case, two meshes (coarse and fine) have been used (Table 3).The core resolution of the non-reacting case is 1.0 and 0.6 mm for the coarse and fine mesh, respectively.
On the other hand, the three meshes (Table 4) adopted in the reacting cases have slightly more cells in the core flow compared to the non-reacting ones.In addition, in all the reacting cases, there is a cylinder-shaped refinement region starting from the cylinder head with the height of 4 mm and diameter of 44 mm to accurately capture the spark region.Moreover, the core resolution of coarse, medium and fine reacting cases are set as 0.75 mm, 0.6 mm and 0.4 mm, respectively.In the current implementation, the cell layers are removed/added during the compression/expansion and hence the core resolution is maintained relatively constant.

Table 1
Division of temperature and composition fields for the present s L correlation.Fig. 3 shows the mesh structure for the reactive cases.A cut-plane of the geometry illustrates the refinement region from the side view.In the refinement region, the grid size is the smallest (0.05 mm) close to the top of the cylinder, while it gradually increases towards the bottom of the refinement region up to 0.25 mm.
The mesh sensitivity is demonstrated by Fig. 4 where the temperature and velocity magnitude probability density functions (PDF) are displayed.It is observed that for the velocity PDFs, there is a slight deviation visible between the cases as the time goes on, however the cases are still statistically very similar to each other.The PDFs suggest that even the coarsest mesh could have been used for simulation of the present configuration with a good accuracy, however the intermediate Fig. 2. Description of the non-reacting (according to [48]) and reacting cases: cold flow cycles (Stage I), single intake stroke (II) and compression/combustion cycle (III).System dimensions are expressed in millimeters.mesh (medium) has been adapted for the rest of the study as it represent the fine mesh case PDFs slightly better.

Non-reacting flow before ignition
Following the second objective of the present study, the purpose of the non-reactive cases is to validate the utilized numerical approach.Focal aspects therein include (1) engine-like flow reproduction in a multi-cycle simulation, (2) thermal stratification during compression, and (3) the compression-associated turbulent heat transfer mechanism.Fig. 2 presents the outline of these simulations starting from a reciprocating cold flow case corresponding to the experimental study of Morse et al. [44], which was previously utilized in DNS [45] and LES studies [77] as well.
With respect to typical mean flow and fluctuation profiles, functionality of the present case (with a cell count of × 2.3 10 6 cells at BDC and × 1.1 10 6 cells at TDC) appears to be in line with the finer LES of [77] as shown in Fig. 5. Modifications of the original experimental configuration were introduced by Schmitt et al. [45] to generate engine-like intake (Stage II) and compression (Stage III) processes, and this rationale is replicated here (Fig. 2).
In Fig. 6 (left), HLR-WT (coarse: × 0.43 10 6 cells at BDC; fine:  × 2.3 10 6 cells at BDC) cases show robust functionality with mild grid sensitivity.In addition to changing core resolution (fine: 0.6 mm; coarse: 1.0 mm), scaled wall normal spacing at TDC + Δ w TDC , reaches ca. 10 (fine) and 21 (coarse).Thermal stratification (Fig. 6, right) increases dramatically after mid-compression.It should be noted that reaching such predictions is highly non-trivial for wall models -indeed, simulations with a standard linear/logarithmic law for temperature (not shown here) lead to substantial deviations.
The turbulent wall heat transfer, driven by impinging and ejecting vortical streams, is visualized in instantaneous snapshots in Fig. 7. Fig. 5. Stage I: Phase-averaged mean axial velocity (left) and RMS velocity (right) profiles at 90 °after TDC.The present method (HLR-WT) is referred against DNS [45] and LES [77].Fig. 6.Total heat flux (left) [46] and temperature PDFs (right) [48] throughout the compression stroke of the non-reactive case.
While the smallest scales of the DNS simulation are obviously out of reach for the present computations, the turbulent heat transfer mechanism and the influence of large scales appear to be well captured.Importantly, the strong influence of flow impingement is relatively well reproduced by the wall treatment in the HLR-WT method.

Visualization of premixed combustion
In this part, the occurrence of CCV and an overview of the flame development in the reacting configuration is presented.The premixed flame development in the present configuration is provided in Fig. 8 for Next, similarities and discrepancies from the three chosen cycles (named cycle 1, 2 and 3) will be explored from the viewpoint of the present numerical model (objective 3).Figs. 9 and 10 indicate the temporal temperature field evolution of the three cycles from top and side views, respectively.It is noted that already at 20 CAD AST there are significant discrepancies in the initial spark evolution observed on both cross-sections.Such an observation signifies, in particular, the role of turbulence and temperature stratifications already at earlier times.A Fig. 9. Instantaneous temperature for the three cycles at the middle plane from top view.similar observation has been previously reported by Enaux et al. [28] and Ozdor et al. [11].Also, an early visual inspection of the figures reveals that cycle 2 not only burns the fastest in the beginning but also continues to burn the fastest throughout the cycle.On the other hand, it can be seen that cycle 3 is the slowest burning cycle among the three.It is also interesting to note that such relatively high CCV can be seen even without a geometrical model for the spark plug.

Global analysis of flow statistics
Fig. 11 shows the density-weighted mean of fuel mass fraction, incylinder temperature and pressure, and filtered heat release rate for the three cycles.It is noted that, consistent with Figs. 9 and 10, cycle 2 burns the fastest while cycle 3 the slowest.Clear differences in fuel consumption arise before 10 CAD AST similar to the work by Enaux et al. [28].Consistently, cycle 2 poses the fastest heat release rate (HRR).Considering the geometrically simple set-up, it is remarkable to note such relatively large discrepancies in peak pressures (30-35%) and maximum temperatures (10-15%) in the present numerical simulations.
As an attempt to justify such discrepancies, global statistics of the 3D turbulent flow field are first considered.Fig. 12 shows the global (cell-volume weighted) PDFs of temperature, resolved velocity magnitude, and modeled velocity fluctuation for the unburnt end-gas of the three cycles.The presented PDFs indicate that there is no significant global statistical variation between the cycles which could clearly stand out as the explanation for the noted CCV.However, as Figs. 9 and 10 indicated, the early spark conditions may have a strong role in the observed CCV.Thereby, we next extend the global analysis to more local analysis with the focus on the early times AST.

Local flame front analysis
The focus in this section is to elaborate the occurrence of CCV and to identify its major contributor in the present configuration (objective 4).In the present set-up, the fuel-air stratification remains relatively minute and thereby only thermal stratification and turbulence effects Fig. 10.Instantaneous temperature for the three cycles at the middle plane from side view.
11.The density-weighted mean cylinder CH 4 mass fraction, temperature, pressure and filtered HRR of the three studied cycles.
exist.Also, it has been noted that the combustion rates of these cycles start to deviate from each other rather quickly (less than 10 CADs) after the spark time (t spark = 330 CAD).Therefore, the local instantaneous temperature and flow fields around the spark location at the early times after the spark time may become significant from the viewpoint of CCV.

Flow conditions around the initial spark
Fig. 13 depicts the instantaneous temperature and velocity magnitude fields at the spark plane (wall distance ≈1.02 mm).The three cy- cles clearly pose very different local temperature and velocity values.In cycle 1, the ignition kernel is in the middle of a hot temperature pocket, while in cycle 2 (fastest cycle), the kernel is located between a hot and cold temperature region.On the other hand, the ignition kernel in cycle 3 (slowest cycle) is mostly in rather cold parts of the flow in contrast to cycles 1 and 2. Based on the observed local temperature close to spark, initially the s L of cycle 1 and 3 is expected to be the highest and lowest, respectively.
Next, it is observed in Fig. 13 that at 330 CAD the ignition kernel in cycle 2 (fastest cycle) is exposed to a higher velocity stream compared to the other two cycles.In particular, cycle 3 (slowest cycle) presented the lowest velocity magnitude close to the spark at 330 CAD.Although velocity magnitude is not a direct measure of turbulence, it can be directly linked to the local velocity gradients and hence also to the local strain rates and sub-grid kinetic energy.Consequently, the noted high and low velocities near the initial spark in cycles 2 and 3 are expected to yield higher and lower modeled turbulence fluctuations around the ignition kernel, respectively.Based on this measure, cycle 3 seems to be slowest in terms of both turbulence and s L consistent with the global combustion rates of the three cycles in Fig. 11.However, quantitative measures would be required to confirm these qualitative observations.

Conditional averages in the flame vicinity
Fig. 14 shows the resolved flame area beside the volume averaged temperature, velocity magnitude, modeled velocity fluctuation ( ′ u mod ), s L and s t within a 4 mm band ahead of the flame front in the unburnt gas.The aim of the following analysis is to explain the globally observed combustion discrepancies between the three cycles by the local fluid dynamical and thermodynamic conditions ahead of the flame.The focus in these analyses is the early time instances after the spark as the first changes in the cycles occur during that period.
The average temperature and velocity magnitude plot coincide with the qualitative findings of Fig. 13 regarding the early behavior of the cycles.In Fig. 14, cycle 1 (intermediate cycle) and 2 (fastest cycle) present the highest average temperature and velocity magnitude in the 4 mm unburnt band at the spark time (t spark = 330 CAD), respectively.However, cycle 3 (slowest cycle) indicated the lowest average temperature and velocity magnitude in the chosen band.The initial lower average temperature and velocity magnitude in cycle 3 implies that this is the slowest cycle among the three.However, to justify the behavior of the other two cycles, the role of turbulence must be investigated further.
Next, the modeled velocity fluctuation ( ′ u mod ) graph in Fig. 14 il- lustrates that cycle 2 (fastest cycle) close to the spark time has a higher value in the chosen band compared to the other two cycles consistent with the velocity magnitude plot.This indicates that at the spark time, cycle 2 presents the highest turbulence effect on the modeled burning rate (s t ).On the other hand, higher average temperature in cycle 1 compared to the other cycles at early times AST results in the highest s L values in this cycle.The s t figure exhibits that the higher turbulence in cycle 2 ultimately overcomes the slight temperature deficiency compared to cycle 1.Hence, cycle 2 presents slightly higher s t values compared to cycle 1 at early times AST.It is worth taking into consideration that s t represents only the modeled part of the total burning rate (TBR) of the system.The source term (the term displacing the levelset) in G-equation is governed by both the modeled part (s t ) and the resolved part (the flame area) of the simulation.
Moreover, the flame area plot in Fig. 14 indicates that the deviation between cycles begins before 340 CAD, and cycle 2 and 3 initially represent the fastest and slowest area growth rate, respectively.At the time of ignition, the cycles start with identical ignition kernels (therefore identical flame areas), however the initial s t and resolved velocity variations between the cycles initiate the first changes in the flame area (resolved part of source term) of these cycles.The increased flame area then contributes to TBR even further by increasing the source term.In fact, a larger kernel (with larger flame area) will burn faster compared to a smaller one despite of having identical modeled turbulent velocity (s t ).Therefore, the cycle that initially burns faster, tends to burn faster till the end.
In addition, it is evident that the spark is initially more susceptible to the local field variations, as even a small difference in the fields can represent a substantial change in the unburnt conditions seen by the spark.Moreover, the resolved (large-scale) flow structures are initially larger than the spark kernel, however as time advances and the burnt region grows, the cases become more similar statistically.Such a trend can be seen in Fig. 14 as the graphs generally become more similar as the kernel is growing.
It is also interesting to see how the spark location can affect the thermal and flow fields.In Fig. 15, the instantaneous velocity magnitude and temperature fields from cross-sections at the spark plane (left column) and at the middle plane (right column) are shown.It is evident that closer to the cylinder head (and near-wall regions), thermal stratification increases due to turbulent wall heat transfer.Therefore, moving the spark position farther away from the cylinder head will further reduce the thermal field impact on the CCV.This is an important consideration especially for the engines with deep penetrating spark plugs.
For the velocity field however, the wall-tangential velocity root mean squares (RMS) reach a maximum somewhere near the wall, whereas the wall-normal RMS increases towards the core (not shown).Such a trend has been also observed in the DNS study by Schmitt et al. [47].Thus, no general conclusion can be made regarding the relation between the spark position and magnitude of flow field variations in the present set-up.Fontanesi et al. [35] indicated that even in the presence of a spark plug, the velocity direction does not introduce statistically significant combustion variations, while the velocity magnitude plays a major role in the promotion or hindrance of the flame progress.On the other hand, a counter argument has been provided by Truffin et al. [36].They reported that the velocity direction plays a significant role in the flame propagation process.They reported that the velocity direction can promote or hinder the flame propagation at early times AST by pushing the flame out or inside of the spark plug cavity, respectively.This issue however is out of the scope of the current work as there is no physical spark plug present in the configuration.
In summary, global statistical analysis was found to be a somewhat non-descriptive measure of CCV, while near spark flow conditions and local flame front analysis at early times largely explained the CCV formation.Such observations indicate that local fields, in particular, at early times can have a significant impact on CCVs.

Asymmetric combustion and quadrant analysis
In this part, the asymmetric character of combustion will be discussed and the previously presented global statistical analysis is extended to more localized quadrant analysis (objective 5).Fig. 16 shows a cut-plane of the temperature field (left) and the burnt volume in each of the four quadrants of cycle 1.It is evident that the flame in sector 2 reaches the cylinder liner the fastest.Fig. 17 shows the PDF of modeled velocity fluctuation ( ′ u mod ), re- solved velocity magnitude and s t for the four sectors at 330 and 335 CAD as well as the unburnt temperature PDF at 181, 280, 300, 330, 340 and 370 CAD.The s t PDFs are obtained from a narrow 4 mm chosen band in the unburnt region ahead of the flame, while all the other PDFs are calculated from the whole unburnt region.The ′ T u , mod and s t PDFs fail to explain the clear trends of the burnt volume in all the quadrants.However, the velocity magnitude PDFs show that the velocity magnitude in sector 2 has slightly higher values compared to the other sectors consistent with the burnt volume graph in Fig. 16.This suggests that velocity magnitude may have a significant role in the asymmetric flame behavior.Hence, in the next step, the convection effect on the flame formation is investigated further.Fig. 18 depicts the temperature field in addition to the narrow band ahead of the flame colored with the sum of the s t and the convection velocity normal to flame front (U n ) which hereafter we denote it with α ( = + α s U t n ).In Fig. 18, α is illustrated for 370 and 378 CAD at 1.02 mm away from the cylinder head.At the flame front, s t indicates the radial growth rate of the burnt volume normal to the flame area.The role of convection velocity normal to the flame surface (U n ) is to deform the flame front.Note that U n does not directly increase the burnt volume, however it will cause the wrinkling of the flame and overall increase of the total flame surface area.
In Fig. 18, it is observed that the flame areas with lower α values (see black markers 1-4) at 370 CAD are retarded at 378 CAD compared to the flame surfaces with higher values (see red markers 5-8).As an example, consider marker 2 at 370 CAD which represents a very low α value due to the low temperature region and also convection velocity moving towards the flame surface in the unburnt mixture at this time.As a result, in the next shown time instance (Fig. 18 surface results in the observed asymmetry.
The previous qualitative analysis indicated that the flame shape can be largely understood by α.Next, to provide quantitative insight, it is useful to consider the correlation between α and radial flame distance to the initial spark location.Fig. 19 demonstrates the flame front distance to the initial spark position as a function of α value at 370 and 378 CAD.It is observed that the cells with higher α values are generally farther away from the initial spark location.It should be noted that the linear least-square fit is obtained from the whole dataset, while the scatter data is shown for every one in 25th cell data.
In summary, it has been noted that in the present set-up, the convection velocity has a higher contribution to the α values (therefore asymmetric flame propagation) simply due to the overall higher magnitudes of convection velocity compared to the modeled s t .Similar behaviors have been seen in planes farther from the cylinder head, but for the brevity reasons they are not shown.

Sensitivity analysis of the modeled spark kernel size
The previous simulations were carried out with a rather small initial spark kernel radius.Next, brief insight to the occurrence of CCVs will be discussed in the present numerical model.To perform such analysis, the  the combustion progress arises essentially from the total flame area.For a small kernel, the initial flame area growth rate may occur slowly if the kernel remains in a low turbulence location.In contrast, for a large kernel the flame area growth rate starts rather rapidly for two reasons.First, the larger kernel has initially a larger flame area indicating already faster combustion.Second, as the kernel spans typically always high and low turbulence regions, the wrinkling onset takes place relatively fast for the large kernels leading to fast surface area growth and quite fast combustion.The provided discussion explains why small kernels are more prone to CCVs than the large kernels in the present numerical model.We note that similar observations on spark size effects were reported in the previous two-dimensional DNS study by Pera et al. [78].

Discussion
The present numerical model has been shown to be consistent and we have been able to better explain various new aspects on CCVs.In particular, the link between turbulent flow and CCVs can now be consistently explained in the present setup.Hence, we consider that the present model is also useful in drawing certain conclusions on the practical relevance of the results.Since real engine design often relies on modeling, we next discuss the relevance of the results both for modelers and for real engine designers.
For engine CFD modelers, the provided numerical evidence indicates several important aspects on premixed combustion modeling in the LES context.First, the initial flame kernel size is clearly an important parameter which can be used to either amplify or damp the  observed CCVs in scale-resolving simulations.Second, CCVs were shown to occur even in a simplified engine geometry without a spark plug.Most likely, the spark plug complicates the situation even more due to combined effects of vortex shedding, flame-spark plug interaction, and surface heat transfer.Third, special attention needs to be put on laminar flame speed modeling for any combustion model in the lean burn regime.Fourth, here the turbulent flame speed model by Muppala et al. [64] was used.It would be interesting to find deeper insight to the model performance in engine-like conditions using e.g.DNS studies.Fifth, accurate wall heat transfer modeling may be needed in order to capture the in-cylinder temperature stratification accurately.
For real engine designers, the practical relevance of the results consists of the following aspects.First, it was shown that turbulence may largely dictate the flame propagation and CCVs, although temperature stratification was also indicated to have some influence on CCVs.Additionally, higher turbulence (higher RPM) may also further enhance temperature stratification which may decrease local flame speed near the cylinder walls.However, since only a single engine speed (RPM) was studied herein, the overall impact of RPM, and hence turbulence levels, on CCVs remains as an open question.Second, to reduce the CCVs originated from the cyclic variation of thermal field, it is suggested to move the spark farther away from the wall.This is due to the fact that level of thermal stratification turns out to be more pronounced near the cold walls.Third, the initial spark strength was noted to be a key quantity in controlling CCVs in SI engines.Small spark development is highly dependent on any change in the local flow and thermal conditions.
In addition, as perhaps the most important practical observation, we note that it would be essential to understand and to stabilize the early flow conditions near the spark plug for controlling CCVs due to early flame development.Such understanding would require detailed understanding of the engine geometry, valve positioning, tumble and swirl which all influence the local flow conditions around the spark.Since control of turbulent flow can be very challenging, we note that the following strategy could turn out to be successful in reducing CCVs in SI context.If, on a given cycle, the local turbulence level around the spark is large, then keep the ignition timing fixed.In contrast, if on a given cycle, the local turbulence level is small, then ignite the mixture earlier.However, the suggested strategy would require sensing technology which may be challenging to implement in practice.

Conclusions
Here, premixed, lean burn SI combustion was investigated in a simple engine geometry using LES and G-equation combustion model.The results have addressed the noted research gaps and provided numerical insight to CCV formation.The main findings of the paper, as related to the objectives of the paper, are summarized as follows: 1.Although a simple cylindrical geometry has been used, considerable cycle-to-cycle variation was observed for the three cycles which had similar global statistics in terms of temperature and velocity PDFs. 2. It was observed that generally local thermal and flow fields around the spark at early times after the spark timing substantially influence the cycle-to-cycle variation.The local variations are much more influential when the flame kernel is small, since even a small field variation is seen by the spark as a major change in the unburnt mixture conditions.3. The initial differences in the s t and resolved velocity between the cycles lead to the early changes in the flame area which then will further contribute to the total burning rate difference between the cycles.In other words, the cycle that initially burns faster (larger flame area), tends to burn faster throughout the cycle.4. In the present numerical study with relatively low composition stratification, both temperature and flow fields contributed to the observed cycle-to-cycle variation, while the flow field impact appeared to be more dominant.5.In the present configuration, the thermal stratification increases as we move towards the cylinder head.This implies that moving the spark position away from the cylinder head will reduce the thermal field impact on the cycle-to-cycle variation.6. Local analysis of the flow field revealed that the flame asymmetry 20.Pressure and fuel mass fraction plots for the fast and slow burning cycles (cycle 2 and 3) for spark radius 0.1 and 1 mm.can be attributed to the combined effect of s t and convection velo- city.Furthermore, the convective flow effect has been observed to be dominant in the present set-up due to the higher magnitudes compared to the modeled s t .7. Increasing the initial kernel size led to shorter combustion duration, higher peak pressures and smaller cycle-to-cycle variations.
Increasing the initial kernel size in the studied conditions made the kernel more resistant to flow and thermal condition variations leading to lower CCVs in the considered example cases.
The present study dealt with a relatively premixed system, therefore the effect of composition stratification has been neglected.It is suggested that more work is required to understand the role of composition stratification in similar lean configurations.In addition, further investigation is needed to understand the influence of engine parameters such as RPM on the contribution of local fields contribution to cycle-tocycle variation.Moreover, it is important to understand the effect of coherent flow structures such as tumble and swirl on cycle-to-cycle variation in the lean gas engine context.

Fig. 3 .Fig. 4 .
Fig. 3. Mesh structure for the reacting cases from top view (left) and from side view (cut-plane) at BDC (right).

Fig. 12 .
Fig. 12.The PDFs taken from the three chosen cycles: modeled velocity fluctuation ′ u mod [m/s] (top left), velocity magnitude (top right) and temperature (bottom).

Fig. 13 .
Fig. 13.Temperature and velocity magnitude of the three cycles at the spark time (wall distance ≈1.02 mm) at 330 CAD from the top view as even a small difference in the fields can represent a substantial change in the unburnt conditions seen by the spinterestingark.

Fig. 14 .
Fig. 14.Weighted average temperature, velocity magnitude, modeled velocity fluctuation, laminar and turbulent flame speed and flame area obtained from a 4 mm band ahead of flame in the unburnt gas for the three cycles.Cycle 2 and 3 are the fastest and slowest cycles, respectively.

Fig. 15 .
Fig. 15.Temperature and velocity magnitude stratification of cycle 1 at the spark plane (1.02 mm) and middle plane (6.25 mm) at 330 CAD.

Fig. 18 .
Fig. 18.Temperature and = + α s U t n fields at ≈1.02 mm away from the cylinder head.The black markers (1-4) indicate flame areas with lower α values compared to the red markers (5-8).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 19 .
Fig. 19.Distance to the initial spark position as a function of α.The scatter data is shown for every one in twenty fifth cell values, while the linear fit represents the whole data.

Table 2
Engine specifications.

Table 3
Non-reacting compression case description (Stage III).

Table A . 6
The present laminar flame speed performance in different regions.The present s L coefficients for different regions.