© The Author(s), published by EDP Sciences, 2022

1 Introduction

There is a strong interest in the automotive industry to develop high efficient, low emission and low carbon fuels that could be used for both pure ICE based and electrified powertrains such as Hybrid Electric Vehicle (HEV), Plug-in Hybrid Electric Vehicle (PHEV) allowing to reach ambitious goals in the fields of emission reduction in the next three decades [1]. The use of Ultra-Lean Technology is a solution to increase combustion efficiency and to tackle pollutant formation in the design of modern spark-ignition piston engines [2] for light duty vehicles. This technology is currently of interest in the automotive research community, in combination with the use of alternative zero-emission fuels such as hydrogen [3]. One of the main technological issues that affects ultra-lean engines is its capability to achieve stable operating conditions. Ultra-lean means here being close to the flame lean extinction limit. In order to enhance flame propagation, lean-burn engines need high levels of turbulence and the mastering of the interactions of the flame with the turbulence is of crucial importance. In particular high levels of Combustion Cyclic Variability (CCV) can be encountered (IMEP covariance higher than 3%) which can deteriorate the engine operation. These are supposed to be related to the interaction between the turbulent structures and the flame kernel in the early stages of the combustion process, e.g. right after the onset of ignition [4]. For this purpose the use of high-fidelity Computational Fluid Dynamics (CFD) tools can help understand the interactions between the various physical phenomena at stake and to make relevant design choices. At the present time the RANS (Reynolds Average Navier–Stokes) modelling approach is daily used for the design and optimization of internal combustion engines. The whole turbulence spectrum is actually modelled and the approach implicitly assumes the existence of a mean phase-averaged flow. On the contrary, the Large-Eddy Simulation (LES) approach does not rely on the hypothesis of existence of a mean cycle as only a part of the turbulence spectrum is filtered leaving the possibility to directly resolve large eddies which are responsible for the CCV. It has been shown that LES is the only affordable CFD approach that allows to describe the CCV because unlike RANS, it resolves the most energetic turbulent structures which variability induces a deviation from the mean phase averaged cycle [5–7]. The understanding and mastering of these complex flow processes are at the centre of new engine technologies development that include both hydrogenated fuels and very lean mixtures. In this context, the objectives of the present work are two fold:

• The first objective is to validate an LES modelling framework, to predict CCV in spark-ignition engines, accounting for the effect of flame stretch on the laminar flame speed, which is important in diluted/lean flame conditions [8].

• Secondly we intend to apply analysis tools to help understand the underlying mechanisms at the origin of CCV in such conditions.

For this purpose the lean-burn TGR-E/FEV research engine configuration [9], operating at an Air/Fuel ratio equal to 1.83, is retained. For this engine, experimental data are available in cold gas-exchange conditions and for the reactive configuration with only slight differences in the geometry between cold and hot configurations. Some modelling work for these ultra-lean mixtures conditions has already been proposed in [4] and will be used as a basis for the present study which aims to go another step further into the combustion cyclic variability analysis with post-processing tools dedicated to complex turbulent reactive flows. After presenting the numerical modeling and set-up, the results of the study are organised as follows:

• Firstly, the LES numerical prediction of the flow field is validated in cold gas-exchange conditions by comparing the velocity field on the tumble mid-plane with the PIV (Particle Imaging Velocimetry) measurements.

• Secondly, the reactive configuration is investigated using ECFM-LES (Extended Coherent Flame Model for Large-Eddy Simulation) including the recently introduced effect of stretch on the laminar flame burning velocity [9]. Here, ECFM-LES is combined with the Adaptive Mesh Refinement (AMR) technique, available in CONVERGE CFD software, to help improving the interaction of the flame with the resolved scales of the flow. The capability of the LES model to reproduce CCV in very lean conditions is assessed by comparing experimental data and numerical results.

• Thirdly, additional innovative post-processing tools are applied to get insights into the correlation between CCV and the turbulent flow field variability, namely Empirical Mode Decomposition (EMD) technique [10] and a tumble shape characterisation analysis tool (Γ_3p)[11].

2 Engine configuration and numerical set-up

2.1 Engine configuration

The engine configuration investigated in this study is the single-cylinder research engine jointly developed by TGR-E (TOYOTA GAZOO Racing Europe GmbH) and FEV [9, 12]. Engine specifications are shown in Table 1. Among the three available engine operating points, the one characterized by leaner mixture (λ = 1.83) and having a transverse position of the spark plug is retained because it is characterized by higher level of CCV. For more details about the engine configuration the reader can refer to [4, 9, 12].

2.2 Physical and LES numerical models

The version 3.0 of the commercial CFD code CONVERGE [13] is used in all computations. A second-order central difference scheme is used for spatial discretization and second-order Crank–Nicolson for time advancement. The Dynamic Smagorinsky model [14] is retained to close the sub-grid scale stresses tensor. Spray injection is modeled with a Lagrangian formalism based on parcels [13]. Primary break-up is modeled prescribing, at the nozzle outlet, a Rosin–Rammler distribution for the droplet size imposing the Sauter mean radius. R3P7. The Sauter mean radius is chosen equal to 50 μm which corresponds to half the injection nozzle radius. Secondary break-up is modeled with KH (Kelvin–Helmholtz)–RT (Rayleigh–Taylor) [15] while droplet evaporation rate is computed with the Frössling correlation [16]. Fuel properties are modeled using the surrogate (C₆₇₆₇H_11.767O_0.07) already proposed for the same configuration in ref. [9]. The injection modeling set-up is a priori validated in a closed vessel computation by comparing liquid penetration, spray footprint and lateral spray views with numerical computations according with the approach proposed in [9]. R3P7. Spray numerical results along with their comparison with experimental data in the closed vessel are shown in Appendix A. The turbulent combustion model used in this work is ECFM-LES for flame propagation [17], while spark plug ignition is described by ISSIM-LES (Imposed Stretch Spark Ignition Model) [18].

The ECFM model used is a modified version of the reference version, called species-based ECFM [19]. The progress variable is not transported but directly calculated as:(1)where is the filtered progress variable, the filtered species mass-fraction with superscript u for “unburned” and b for “burned”. The filtered species mass-fraction are described by the following transport equations:(2)Where  and  denote Reynolds and Favre filter respectively. is the filtered density and X represents either u or b. ν is the molecular kinematic viscosity and is the subgrid-scale kinematic viscosity estimated at the combustion filter size  = nresΔ_x. Δ_x is a typical cell size and nres is set to 8 in the present study [4, 7]. The source term takes into account the reaction rate and is computed separately for burned and unburned species.

For unburned species, without considering the auto-ignition combustion mode contribution, the source term is computed as:(3)where ρ_u is the density of unburned gases, S_l the laminar flame speed and the Flame Surface Density (FSD).

For burned species, the source term , without considering the auto-ignition combustion mode contribution, is evaluated as:(4)where is the burned gases composition following [20], the post-oxidation reaction rate in burned gases defined by simplified kinetics [19].

In the above two equations, two terms remain to be defined: the FSD and the reaction rate . The last one is closed by the ISSIM model [18] accounting for the energy deposition due to spark ignition.

The resolution of the transport equation of the FSD is numerically unstable [17] and thus a modified FSD is proposed, being the local vector normal to the iso-surfaces of the progress variable. The corresponding transport equation is detailed in [4, 7, 17].

Following the work of Hawkes and Cant [21], the local flame normal direction is replaced by the average normal direction of the flame front defined as follows:(5)

Equation (5) takes into account the subgrid scale fluctuations of the vector orientation. The subgrid scale wrinkling is defined following:(6)with the laminar part of the FSD [17]. The factor 1/Ξ_sgs makes the laminar contribution disappear when the flow becomes very turbulent.

The transport equation for the filtered flame surface density is finally modelled as follows:(7)where is the displacement speed of the flame and the expansion factor. The contributions from the resolved part of the flame are the convection T_res, the laminar propagation P, the filtered curvature C_res and the filtered tangential strain rate S_res. The assumption of isotropic propagation directions of the flame front being not justified when the flame front is well resolved, the formulation of Hawkes and Cant [21] for the resolved strain rate S_res is used in the present study:(8)where(9)with . The term α_n makes the isotropic part of the orientation tensor disappear when the flow becomes laminar or very well resolved, while it becomes important when the flow is very turbulent or poorly resolved.

T_sgs, C_sgs and S_sgs are respectively the subgrid-scale turbulent transport, the subgrid-scale curvature and strain rate terms. The last term S_sgs includes a calibration parameter α_CFM and an efficiency function Γ developed by Bougrine et al. [22] which models the interaction between the resolved eddies of scales smaller than and the flame front.

The initial flame kernel is typically smaller than the local grid size. The ISSIM model [18] allows to model the transition from the SGS flame kernel to the resolved one with the help of a transition factor α(x, t) which is equals to 0 before ignition and gradually increases to unity when ignition is ended. The factor α(x, t) is itself a function of the local equivalent burned gases radius r_b obtained by a model transport equation (see [18] for a complete description of ISSIM). Equation (7) is then modified to integrate the ignition modeling:(10)where is the ISSIM modeled laminar stretch that replaces the resolved contributions C_res + S_res during the ignition phase [23].

The laminar flame speed S_l is computed from correlations as a function of the thermodynamic (pressure and temperature) and mixture engine conditions according to the kinetic model proposed in [24]. To account for the resolved flame stretch effects on the laminar flame speed [25] a correlation of Markstein number is built up according to the methodology proposed in [4]. This correction, lowering the mean laminar flame speed in lean conditions, allows to retrieve the experimental pressure signal without modifying the experimental spark advance, as already pointed out in [4].

2.3 Computational grid

In CONVERGE a Cartesian cut-cell approach is used for mesh generation. The computational grid size is controlled through the use of embedding to customize the cell size in the various engine parts and phases according to the following strategy:

• Cylinder: between 0.125 and 1 mm.

• Intake port: between 0.5 and 1 mm.

• Exhaust port: between 0.5 and 2 mm.

In the fired simulation, in addition to a local mesh refinement around the spark-plug, the Adaptive Mesh Refinement (AMR) is used to locally refine the computational grid by a factor two, leading to a 0.25 mm resolution, to improve flame wrinkling at the resolved level, without non-necessary refinement outside the reaction zone. The AMR criterion used is simply to add an AMR level when the progress variable lies between 1.e–5 and 0.9.

Figure 1 shows that the adaptive mesh refinement follows the flame front encompassing part of the fresh and burned gases.

Fig. 1 Progress-variable field on the tumble mid-plane over-imposed on the computational grid at −4 Crank Angle Degree (CAD) from the combustion top dead centre.

Figure 2 shows the evolution of the total cell number (in millions) along the whole engine cycle: it is bounded between 0.5 during the exhaust phase and 6 during the intake phase. It can be observed that the AMR for ECFM has a moderate impact on the total number of cells only during the combustion phase while it allows to significantly improve the flame front resolution.

Fig. 2 Evolution of the total cell number (intake/exhaust manifold plus cylinder) along the cycle with and without progress variable based AMR (0 CAD corresponds to the combustion top dead centre).

2.4 Boundary conditions

Time varying intake pressure and temperature signals as well as wall temperatures are prescribed in accordance with previous work [4, 9]. The total injected fuel mass per cycle is equal to 4.88 mg. This set of boundary conditions allows to recover the correct in-cylinder trapped mass and target Air/Fuel ratio estimated with the GT-Power 1D model of the engine. The spark advance is imposed equal to the experimental one (−20.3 CAD after top dead centre).

3 Numerical results

3.1 Gas-exchange cold flow

In a first phase, LES numerical modeling is validated for the prediction of the in-cylinder flow field. Indeed, the correct prediction of the in-cylinder aerodynamic motion (average flow field and its RMS) is of primary importance for the consequent description of the CCV in reactive conditions. The numerical results collected over 10 cycles are compared with PIV experimental data collected over 38 cycles on the tumble mid-plane. The first numerical cycle was discarded to remove the effect of the initialization.

Figure 3 compares the numerical mean velocity field to the experimental one during the intake (−270 CAD from the combustion top dead centre) and during the compression strokes (−150 CAD from the combustion top dead centre). The 2D colormaps show a good reproduction of the in-cylinder motion including a precise prediction of the tumble centre location evolution. A slight overestimation of the velocity level is visible from the velocity colormaps, especially during the intake stroke. This overestimation is in accordance with a previous literature study conducted on the same configuration [9]. For −150 CAD, 1D mean velocity profiles are also extracted from the 2D field and compared in Figure 4 to confirm the quality of the LES results. RMS are also compared in the Figure 4 showing that, despite the low number of cycles makes the statistics poorly converged, the fluctuation level is globally well captured.

Fig. 3 Comparison of the phase-averaged mean velocity field on the tumble mid-plane between PIV measurements (left) and LES (right) at −270 and −150 crank angle degrees from the combustion top dead centre. 2D velocity magnitude is normalized by the mean piston velocity.

Fig. 4 1D Mean and RMS velocity profiles extracted at 4 vertical positions with respect to the centre-line (−30, −10, 10 and 30 mm).

3.2 Reactive calculation

After the validation of the flow field prediction, the fired configuration was computed in LES. As for the gas-exchange calculation, a first cycle is conducted to initialize the field. Subsequently, 30 consecutive cycles were run. Even if not enough to reach the statistical convergence, the 30 cycles are sufficient to start the CCV analysis in accordance to what is classically done in literature [26, 27]. The computational time for each cycle is about 48 h on 180 CPUs. This computational cost is strongly mitigated by the use of AMR in CONVERGE. Figure 5 compares the experimental and numerical pressure envelopes. Numerical pressure signals are correctly located in the experimental envelope except for a slight earlier start of combustion. From the multi-cycle traces it is observable that, despite the numerical pressure results show a high cycle-to-cycle variability, the extreme experimental events are not reproduced. However, the occurrence of such extreme cycles is rare. In fact, conducting a statistical analysis on the pressure traces, it is possible to observe in Figure 5 that mean and standard deviation (σ) are very well reproduced by numerical simulations. The mean pressure peak is correctly identified except for a very small shift of its location, while the amplitude of the standard deviation is slightly underestimated. This can be either due to the low number of cycles or to some physical phenomena that are not accounted in the present computation (temperature variations, auto-ignition, spark energy variation, etc.). Figure 5 shows also the evolution of the in-cylinder peak pressure value along the 30 cycles in order to demonstrate that each cycle is completely uncorrelated with the previous and the following one.

Fig. 5 Upper figure: LES multi-cycle in-cylinder pressure traces compared to the experimental ones: the max and min of the experimental envelope are also highlighted. Middle figure: Mean and Mean ± σ of the pressure traces for LES and experiment; numerical pressure statistics are computed over 30 LES cycles while the experimental ones over 1000 cycles. Lower figure: maximum cylinder pressure for the 30 consecutive LES cycles.

3.2.1 CCV analysis

To further analyse the CCV and the performance of the LES model with respect to the experiments, Figure 6 shows the scattering of the CA50 (middle of combustion) versus the CA2 (start of combustion). The numerical scatter exhibits the same linear behaviour as observed in the experiments and the amplitude of the scattering is similar. As already observed on the pressure traces, the cycles showing the earliest start of combustion (CA2) in the LES are roughly 2 CAD in advance compared to the equivalent cycles in the experiment. Furthermore, Figure 6 displays the Matekunas diagram [28] which scatters the maximum in-cylinder pressure (P_max) for each cycle as a function of the Crank Angle (CA) at which P_max occurs. LES results follow as the experiments a globally linear trend, but the experimental distribution shows two linear bands, an upper one and a lower one, while the LES only reproduce the latter. This difficulty to reproduce extreme events is in line with the underestimation of the in-cylinder pressure standard deviation observed in the previous section. The same data are binned, in Figure 6, to obtain a PDF distribution of P_max. Despite the low number of samples for the LES, it confirms the ability of the proposed LES approach to reproduce the global shape of the distribution, but as stated before, also the difficulty to reproduce the extreme cycles which are at the tails of the distribution.

Fig. 6 Multi-cycle analysis on global engine quantities: Upper plot: scattering of the CA50 (middle of combustion) versus CA2 (start of combustion). Middle plot: Matekunas diagram scattering Pmax as a function of the Crank Angle (CA) at which Pmax occurs. Lower plot: Pmax probability density: each bar displays the bin’s raw count divided by the total number of counts and the bin width; the pressure range 70, 120 is divided in 10 bins equally spaced.

In the absence of auto-ignition events, CCV is mainly related to the development of the flame kernel in the earlier crank angles after the spark timing. Figure 7 illustrates the initial flame kernel development through an iso-surface of the progress variable for a fast-burning cycle (17) and for a slow-burning one (20). The peak pressure, shown in Figure 5, is used as a criterion to distinguish fast and slow burning cycles. At CAD = −10 (9.7 CAD after the ignition time), it can be observed that for cycle 17 the total flame surface is larger than cycle 20. To generalize the observation and quantify the correlation between the flame surface magnitude and the pressure peak, Pearson product-moment correlation coefficients are computed and shown in Figure 8 for three different CAD (−10, −5 and 0). The qualitative observation made in Figure 7 is quantitatively supported in Figure 8 by the magnitude of the correlation coefficients.

Fig. 7 Flame kernel development identified by the progress variable iso-surface at 0.5 coloured by the local velocity magnitude. Upper figure: fast burning cycle. Lower figure: slow burning cycle.

Fig. 8 Scattering of the total flame surface versus the peak pressure for all LES cycles. The total flame surface is computed at −10, −5 and 0 CAD from the combustion top dead centre. The cycle number is shown for each dot in the plot. The correlation coefficient between the two variables is shown in the plot title (r).

The development of the initial flame kernel can be affected mainly by two phenomena:

• The flow field, and consequently, the turbulence intensity around the spark-plug at the ignition timing.

• The mixture distribution around the spark plug.

By applying the multivariate linear regression model proposed by Truffin et al. [6] to simultaneously evaluate the effect of multiple variables, it is found that the main source of variability for the present engine configuration, comes from the flow field. In particular, using the flame surface magnitude at −5 CAD as the response variable, three independent variables are found in the best fitting model:

• The high frequency kinetic energy around the spark plug coming from EMD analysis (EC-HF-SP), detailed in the following section.

• The velocity field in the X direction (UX).

• The mean Equivalence Ratio around the Spark Plug (ER-SP).

For Spark-Plug (SP) quantities, a volume average operation in a sphere of 3 mm radius around the spark plug is considered. Figure 9 shows that the main impact on the response quantity comes from the aerodynamic variables (EC-HF-SP and UX). In particular, as will be further explained in the next section with the EMD analysis, it is found that the high frequency kinetic energy around the spark plug and the in-cylinder velocity field have higher impact on CCV than mixture distribution.

Fig. 9 Correlation coefficients computed with the linear regression model proposed in [6]. The first bar represents the performance of the best fitting reduced model to reproduce the response variable (flame surface at −5 CAD). The correlation coefficient between flame surface and the best fitting model is equal to 0.76. The other bars show the correlation coefficients between flame surface and each variable independently from the other variables: 0.72 for EC-HF-SP, 0.49 for UX and 0.34 for ER-SP.

This result was expected as for the present engine configuration, the mixture stratification around the spark plug is moderate since the injection happens early (end of injection at −269 CAD); consequently fuel droplets have enough time to evaporate and homogenize in the cylinder volume. To support this analysis, Figure 10 shows the equivalence ratio distribution on the spark cut-plane for a fast and a slow burning cycle. The ϕ variation in the whole chamber is limited to ±0.08 around the mean value (ϕ = 0.55). Furthermore the two cycles show similar ϕ values in the spark plug vicinity.

Fig. 10 Equivalence ratio distribution at −22 CAD on the tumble and spark mid-plane. The colorbar is clipped between the maximum and minimum value found in the ϕ field.

To further investigate the role of the flow field on the CCV, the following section proposes new complementary analysis tools that can be used to post-process 3D CFD LES engine results.

4 Flow analysis tools

4.1 2D multivariate Empirical Mode Decomposition

The Empirical Mode Decomposition (EMD), first introduced by Huang et al. [29] for univariate time series, is a data-driven method for analyzing the nonlinear and non-stationary signals. It decomposes a signal into a finite set of Amplitude and/or Frequency Modulated (AM/FM) components, called Intrinsic Mode Functions (IMFs), and a residual component. It was extended for multivariate signals, firstly for complex signals [30, 31] and then for Bivariate (BEMD) [32], Trivariate (TEMD) [33] and finally Multivariate (MEMD) [34]. The basic idea of all these approaches is to project the multivariate signal in certain uniformly sampled directions onto a circle (BEMD), a sphere (TEMD) or an hypersphere (MEMD), and then estimate an envelope of the signal which is required in the EMD algorithm.

The mode mixing problem was also inherited from the original EMD to the MEMD, and a Noise Assisted MEMD (NA-MEMD) was thus developed [35]. Another improvement was made in Adaptive-projection intrinsically transformed MEMD (APIT-MEMD) [36] in which the imbalance between different components of the multivariate signal was taken into account. The 2D multivariate EMD of the present study is a combination of the former two methods that was proposed by Sadeghi et al. [37] and noted as Ensemble NA-APIT-MEMD. It was applied to a velocity field constructed with a large-scale organised synthetic vortex and relatively small-scale isotropic homogeneous turbulence. Subsequently, the methodology was also applied to high-speed PIV measurements in the symmetry plane of the tumble of an optically accessible engine [10]. Both results confirmed that the turbulent part could be separated from the large-scale organized motion based only on the instantaneous velocity field, following:(11)with U_LF the large-scale motion of low-frequency and U_HF the turbulent fluctuations of high-frequency.

In [10], the EMD method has been adapted to the analysis of PIV fields in complex, time-dependent geometries such as those encountered in engine configurations where the size of the velocity field changes with the piston position. R2P1. As confirmed in [10], EMD is suitable for the analysis multi-cycle flow fields in engines as it is not limited by a statistical convergence criterion. Other decomposition methodologies exist in the literature. For example, the Proper Orthogonal Decomposition (POD) is a statistical method that has been widely applied to in-cylinder velocity data, obtained from LES and mainly from PIV measurements by a number of groups [38–42]. It allows the velocity field to be decomposed into a set of orthogonal basis functions. However, the method needs numerous instantaneous velocity fields to obtain a statistical convergence. On the contrary, the EMD technique, capable of decomposing each instantaneous field, can yield meaningful cycle-based analysis. It is therefore more suitable for the analysis of the few cycles obtained by LES.

The EMD method used in this study is modified and improved in terms of the mode combination strategy. In [10], an imposed number of modes N_IMF is predefined for the decomposition in both horizontal and vertical directions and thus the velocity field needs to be interpolated into a uniformly sampled square grid. Once the decomposition is completed, pseudo IMFs are combined into N_IMF IMFs, and finally only the last IMF and the residual are retained in the Low-Frequency (LF) part which represents the large-scale motion. The sum of the remaining IMFs is considered as the High-Frequency (HF) part that is related to turbulent fluctuations. To obtain LF and HF parts, the complete decomposition is in fact not necessary and the method can be therefore accelerated as summarized below:

1. Apply the Ensemble NA-APIT-MEMD to each row (or column) of the bi-dimensional instantaneous velocity field U(x, y), i.e., decompose along one dimension. Sum the last IMF and the residual as “pseudo LF”.

2. Apply the Ensemble NA-APIT-MEMD to each column (or row) of the “pseudo LF” resulted from the last step, i.e., decompose along another dimension. Sum again the last IMF and the residual as LF part U_LF(x, y) of the input signal.

3. The HF part U_HF(x, y) can be extracted directly by subtracting U_LF(x, y) from U(x, y).

Using this approach, the last pseudo IMF and the residual are used in the decomposition in each direction in order to get the LF part. The HF part is then simply deduced subtracting the LF field from the instantaneous field, which allows a speed-up of a factor (N_IMF + 1)/2. Moreover, the number of IMFs N_IMF is not a priori defined, as we only focus on the construction of the LF part instead of all the individual IMF modes. Consequently, the mode combination strategy with predefined number of modes used in [10] is not anymore necessary. Such an improvement in the computational cost of the EMD method allows its massive application to the analysis of in-cylinder flows during the engine cycle. R3C3. The accuracy and the CPU performances of the modified EMD methodology are evaluated and compared versus the original one for a single engine cycle. The results of this comparison are considered satisfactory and the CPU performance is in line with the expectation as shown in Appendix B.

In the present work, the modified EMD tool is applied for the first time to analyse the in-cylinder flow field from 3D LES in order to separate the large-scale organized motion (low frequency part) from the small scale turbulent fluctuations (high frequency part). The objective is to find a relationship between the variability of the low and high frequency parts of the flow field and the variability of the combustion speed.

4.2 EMD results

Bivariate 2D EMD is here applied to analyse the 2D velocity field within the central tumble (and spark-plug) plane. The study focuses on the last part of the compression stroke which characterizes the flow before ignition takes place. Both velocity components (x and z), belonging to the tumble plane, are here considered. Figure 11 shows, for a fast combustion cycle (17) and a slow one (20), the low frequency part of the velocity field and the high frequency part obtained by EMD. The results are extracted at −30 CAD (10 CAD before ignition) and at −22 CAD (2 CAD before ignition) to analyse the flow field before ignition and the resulting flame perturbation. It can be observed that the low frequency part, characterizing the global features of the flow, shows a variability between the two cycles. At −30 CAD, for cycles 17 and 20, an organized tumble motion is identified even though the location of the tumble structure is not the same: more centred with respect to position of the spark for cycle 20. With respect to the high frequency part at CAD −30, similar structures are observed for both cycles. On the contrary, when considering the instants closer to the spark timing (e.g. at −22 CAD), for cycle 20 a well-organized tumble motion is still visible, while for cycle 17 a less well-structured tumble is found. The tumble breakdown leads to an increase in the level of turbulent fluctuations near the spark-plug, which is noticeable in the high frequency part of the EMD.

Fig. 11 Bivariate 2D EMD results for two different CAD (−30 and −22) and for two combustion cycles, a fast one (17) and a slow one (20). The Low Frequency (LF) and the High Frequency (HF) parts are shown. x is the horizontal coordinate and z the vertical one. For the LF part a vector field representation is chosen to better visualize the organized tumble motion.

To generalize the observation to the whole set of computed cycles, the Pearson product-moment correlation coefficients between the spatial average of the high frequency energy around the spark plug at the spark-timing and the flame surface at −10 CAD are computed. The spatial average is performed on a 2D disk of radius 3 mm around the spark plug. The kinetic energy of the HF part is calculated according to [10]:(12)

and are the x and z squared velocity components averaged in a disk of radius 3 mm around the spark position. EC_HF represents the kinetic energy contribution of the high frequency EMD mode averaged over the same area. Physically, it describes the kinetic energy contribution of small non-organized turbulent scales. Figure 12 shows the dispersion of the flame surface at −10 CAD and the HF kinetic energy computed at −22 CAD. The Pearson product-moment correlation coefficient confirms that the variables are correlated and that the magnitude of the turbulent fluctuations around the spark plug has a significant impact on the consequent combustion speed. Furthermore, Figure 12 shows that the addition of the high frequency part of the third velocity component U_yHF (normal to the tumble plane) does not alter the above conclusions.

Fig. 12 Scatter plots of the total flame surface value computed at −10 CAD as a function of the high frequency kinetic energy spatially averaged on a 2D disk of radius 3 mm around the spark plug. For each point the corresponding cycle number is indicated in the figure. The Pearson product-moment correlation coefficients of the two variables are also shown in the title. Upper figure: EC is computed considering the two velocity components within the central tumble plane. Lower figure: the third velocity component (normal to the central tumble plane) is also added in the EMD analysis.

To complete the 2D EMD analysis on the origin of the CCV coming from the flow field, the Γ₃ function is hereafter applied to the flow field, giving access to a 3D visualisation of the tumble evolution before ignition.

4.3 Γ₃ function

The tumble flow is the main large scale rotational motion in the combustion chamber of a spark-ignition engine. Algorithms called Γ₁ and Γ₂ have been developed to identify such a vortical structure by Graftieaux et al. [43]. The criterion is not Galilean invariant but allows a simple and robust method to identify rotational structures, while the Γ₂ criterion is Galilean-invariant because local convection is taken into account. A comparison between these two criteria was performed in [44] and showed that Γ₁ was much more appropriate for locating vortex centres while Γ₂ could help identify larger vortex regions. The Γ₁ criterion has often been applied to in-cylinder flows, although using only 2D instantaneous velocity fields to detect and track the location of the tumble centre (see for example in [10]).

For three-dimensional flows like the tumble flow, Γ₁ is limited due to its 2D nature and thus has been extended in 3D. A very simple extension is to compute Γ₁ separately over several parallel 2D planes throughout the combustion chamber, as presented in [44] on PIV data and more recently in [45] on LES data. A main drawback of this type of extension is that the axis of rotation cannot be obtained, and an improvement was made in [46] in which the author determined the local axis of rotation using a sweeping method: all the possible directions were tested and the one giving the maximal average value of |Γ₁| over the corresponding 2D rotation plane was retained. However this approach being not parallelized by the author was very time consuming and consequently only applied on two instantaneous 3D velocity fields. A real extension of Γ in 3D was realized in [47], named Γ₃ and was applied on external flows to identify vortical structures in 3D volume. Γ₃ was then applied to the in-cylinder flow as presented in [11] with slight modifications: Γ₃ was first computed to obtain the axis of rotation, and then recomputed using the quantities projected in a plane normal to the axis of rotation, denoted as Γ_3p. In this work, a modified version of Γ_3p will be computed to identify the tumble structure in the compression stroke.

Γ₃ is first evaluated inside the combustion chamber using the following expression:(13)with φ the crank angle, V a sub-volume centred at the point x and r(x₀) the distance vector points from x to x₀. (x₀, φ) is the local velocity u(x₀, φ) at x₀ without the contribution from local convection 〈u〉_V (x, φ):(14)〈u〉_V (x, φ) is the average velocity inside V.

The formulation of Γ₃ allows a reliable estimation of the rotation vector at point x by a simple normalization of Γ₃:(15)

Using Equation (15), projections of the local distance and velocity vectors on the rotational plane are obtained as follows:(16) (17)where r_p(x_o) and u_p(x_o, φ) are the projected distance and velocity, respectively. The final formulation of Γ_3p is then given by:(18)

Compared to the definition in [11], u_p is used both in the denominator and numerator in the current study in order to keep a more consistent expression. The magnitude of Γ_3p indicates whether an organized rotation exists: the unit value is related to a perfect vortex without any deformation while a small value close to zero corresponds to the non-existence of rotational structures. On the basis of Γ_3p, the tumble centre can be easily extracted as in the 2D case: along the main axis of the tumble motion (Y-axis in this study), a local centre can be determined where the magnitude of Γ_3p reaches its maximum. All the local centres form a Tumble Centre Line (TCL) around which the tumble rotates.

Γ_3p does not provide the rotation intensity. For that purpose a local tumble ratio was defined in [11] along the TCL at each centre x_c using a local average vorticity. In this study we propose another formulation to quantify the angular velocity around the TCL:(19)with V_r representing the volume including only points very close (distance less than 0.5 mm) to the rotational plane of the centre x_c. The computed angular velocity is normalized by the engine rotational speed ω_ref in rad/s. This formulation is quite similar to that of the global tumble ratio except that only points close to the rotational plane are used with their projected quantities and that the tumble ratio is computed around the local centre instead of the global centre of gravity. TR here is also a vector and thus provides the rotation intensity in different directions.

Thanks to the Γ_3p criterion, the tumble centre line and the corresponding local Tumble Ratio (TR), we are able to describe precisely the tumble structure in the chamber. Its evolution during the compression can provide information about the tumble deformation or breakdown.

4.4 Γ₃ results

Γ_3p vector field is computed for cycle 17 and cycle 20 during the compression stroke to visualize the 3D evolution of the tumble structure and its orientation. Figure 13 shows the obtained results: Γ_3p function allows to visualize the evolution of the main tumble vortex in 3D. The main vortex reduces its size approaching toward the top dead centre. It is also possible to observe that the position of the tumble centre and its local axis of rotation are not perfectly aligned with those of a perfect tumble rotating around the Y-axis and this behavior varies from cycle to cycle.

Fig. 13 Γ3p norm iso-surface at 0.5 coloured by its Y direction magnitude. The local direction is also shown over the iso-surface by Γ3p vector field. Cycle 17 and cycle 20 at three different CAD (−65, −30 and −22) are considered.

Fig. 14 Γ3p norm on the spark mid-plane and 3D vector field inside its 3D iso-surface at 0.6. Cycle 17 and cycle 20 at CAD −22 are considered.

In agreement with the EMD analysis, it is also possible to notice that for cycle 20, the iso-surface |Γ_3p| = 0.5 still retains an organised structure at −22 CAD whereas for cycle 17, tumble destruction occurs. The faster tumble breakdown generates locally more turbulent kinetic energy near the top dead centre, which favours the growth and propagation of the flame kernel. Moreover, a main tumble structure still exists along the Y-axis of cycle 20 at −22 CAD, unlike cycle 17, and its location on the central mid-plane is close to the spark-plug.

Using the function Γ_3p, the local tumble ratio along the TCL is also extracted. Results are shown in Figure 15 for cycle 17 and 20 at −55 CAD and −30 CAD. As the main axis of rotation is aligned with the negative Y coordinate, the value of the tumble ratio is always negative. Along the TCL, at −55 CAD, the tumble intensity is higher for cycle 17 than for cycle 20 in absolute value, while for −30 CAD the opposite tendency is observed. For the same cycles, Figure 15 also shows the evolution of the average local TR in the central zone corresponding to the spark plug coordinate along the compression stroke. This average value is compared to the global TR computed in the whole cylinder volume. By observing the evolution of the local average TR, a crossing point around −35 CAD is visible which justifies the observation made on the evolution of the local TR along the TCL. For cycle 17, between −30 CAD and −20 CAD, it is possible to observe a discontinuity in the plot, for the local average TR, which corresponds to the tumble breakdown observed in Figures 13 and 14. The local tumble evolution is globally in good agreement with the global tumble evolution, but as expected, the absolute value is slightly different. The intensity of the tumble decay may also play an important role on the tumble breakdown and consequently on the combustion speed as already shown in the literature [48]. In order to generalize this qualitative observation to the whole set of cycles, from the analysis of the local and global TR, a correlation coefficient between the slope of the global TR and the value of the flame surface at TDC for all cycles is also computed and is equal to 0.62. This value confirms the impact of the tumble decay on the combustion cyclic variability.

Fig. 15 Upper/middle figure: evolution of the local tumble ratio TR, at −55 CAD and −30 CAD for cycles 17 and 20, along the TCL identified with its Y coordinate. Lower figure: comparison between the local mean TR around the spark plug (SP) with the global tumble ratio for cycle 17 and 20 along the compression stroke between −95 and −20 CAD. The mean operation is made considering an enclosing volume of ±3 mm in Y direction.

5 Conclusion

A LES modelling framework, combining the ECFM-LES approach and AMR is proposed, validated and exploited with specific analysis tools in this work. Firstly, LES results were validated for flow field prediction and secondly for combustion prediction for the ultra-lean spark ignition TGR-E-configuration. For this engine, detailed experimental results were available. The proposed LES modeling approach allows to achieve a good compromise between accuracy and computational cost by using AMR. Furthermore, a satisfactory statistical representation of the cycle-to-cycle variation is obtained with very encouraging perspectives to make this approach suitable to perform parametric variations on realistic configurations. Using a linear reduced order model it was found that, for the analysed configuration, the in-cylinder aerodynamic variations constitute the main origin of CCV. To further analyse these variations two original data analysis tools (Bivariate EMD2D and Γ_3p function), dedicated to the in-cylinder aerodynamic analysis, were applied to post-process LES engine results. It was found that the tumble structure evolution and its subsequent breakdown before ignition have a primary importance on combustion velocity and the proposed tools help to find the origin of this transition and to bring quantitative elements of analysis. In particular, the tumble destruction promotes the development of small scale turbulent structures that boost the flame kernel development and consequently the combustion velocity. EMD2D allows to distinguish the difference, among the cycles, in the large scale and small scale structures of the turbulence. In particular, it was found that the small scale structures located around the spark-plug, have a primary importance for the combustion velocity. In addition to EMD, Γ_3p allows to access a 3D visualization of the tumble motion allowing to compute the local tumble intensity decay and visualize the 3D breakdown. However, in the present work, the EMD analysis was restricted to the 2D tumble mid-plane and it is of interest for future work to proceed with the development of a 3D EMD methodology in order to improve the identification of the CCV sources. The combination of dedicated models for complex turbulence/chemistry interactions in lean conditions and specific post-processing tools for the analysis of the effect of turbulence on the flame development would help identify the origins of CCV. This opens new perspectives in the mastering of these effects for modern engine development operating under lean conditions with strong levels of turbulence such as for hydrogen engines.

Acknowledgments

The authors acknowledge TGR-E (TOYOTA GAZOO Racing Europe GmbH) for sharing the experimental configuration and measurements.

References

Appendix A

Closed vessel spray numerical results

The spray numerical set up is a priori validated in a closed vessel configuration for which injector experimental characterization is available.

Figure A.1 shows the comparison between the absolute penetration of each of the seven holes numerical jet and the maximum experimental one and the comparison between experimental side view shadowgraphy and numerical spray projected area.

Fig. A.1 Upper figure: 1D plot comparison between experimental spray penetration (black dots) and numerical penetration issuing the 7 holes injector. Lower figure: 2D colormap comparison between experimental spray Shadowgraph and spray projected area 1 ms after the 1st droplet gets out of the nozzle.

Appendix B

Modified EMD tool

Figure B.1 shows that the agreement between the original and accelerated EMD methodologies is overall good. Very small discrepancies can be observed between the two field for both LF and HF fields keeping the global error below 4%. Furthermore, it was verified that the computational time speeds up by a factor of 7 and this result is in line with the expectations.

Fig. B.1 Comparison between the original (upper figures) and the accelerated EMD methodology (lower figure) for cycle 2 and CAD = −30.

All Tables

All Figures

	Fig. 1 Progress-variable field on the tumble mid-plane over-imposed on the computational grid at −4 Crank Angle Degree (CAD) from the combustion top dead centre.
In the text

	Fig. 2 Evolution of the total cell number (intake/exhaust manifold plus cylinder) along the cycle with and without progress variable based AMR (0 CAD corresponds to the combustion top dead centre).
In the text

	Fig. 3 Comparison of the phase-averaged mean velocity field on the tumble mid-plane between PIV measurements (left) and LES (right) at −270 and −150 crank angle degrees from the combustion top dead centre. 2D velocity magnitude is normalized by the mean piston velocity.
In the text

	Fig. 4 1D Mean and RMS velocity profiles extracted at 4 vertical positions with respect to the centre-line (−30, −10, 10 and 30 mm).
In the text

Fig. 5 Upper figure: LES multi-cycle in-cylinder pressure traces compared to the experimental ones: the max and min of the experimental envelope are also highlighted. Middle figure: Mean and Mean ± σ of the pressure traces for LES and experiment; numerical pressure statistics are computed over 30 LES cycles while the experimental ones over 1000 cycles. Lower figure: maximum cylinder pressure for the 30 consecutive LES cycles.

In the text

Fig. 6 Multi-cycle analysis on global engine quantities: Upper plot: scattering of the CA50 (middle of combustion) versus CA2 (start of combustion). Middle plot: Matekunas diagram scattering Pmax as a function of the Crank Angle (CA) at which Pmax occurs. Lower plot: Pmax probability density: each bar displays the bin’s raw count divided by the total number of counts and the bin width; the pressure range 70, 120 is divided in 10 bins equally spaced.

In the text

	Fig. 7 Flame kernel development identified by the progress variable iso-surface at 0.5 coloured by the local velocity magnitude. Upper figure: fast burning cycle. Lower figure: slow burning cycle.
In the text

	Fig. 8 Scattering of the total flame surface versus the peak pressure for all LES cycles. The total flame surface is computed at −10, −5 and 0 CAD from the combustion top dead centre. The cycle number is shown for each dot in the plot. The correlation coefficient between the two variables is shown in the plot title (r).
In the text

Fig. 9 Correlation coefficients computed with the linear regression model proposed in [6]. The first bar represents the performance of the best fitting reduced model to reproduce the response variable (flame surface at −5 CAD). The correlation coefficient between flame surface and the best fitting model is equal to 0.76. The other bars show the correlation coefficients between flame surface and each variable independently from the other variables: 0.72 for EC-HF-SP, 0.49 for UX and 0.34 for ER-SP.

In the text

	Fig. 10 Equivalence ratio distribution at −22 CAD on the tumble and spark mid-plane. The colorbar is clipped between the maximum and minimum value found in the ϕ field.
In the text

	Fig. 11 Bivariate 2D EMD results for two different CAD (−30 and −22) and for two combustion cycles, a fast one (17) and a slow one (20). The Low Frequency (LF) and the High Frequency (HF) parts are shown. x is the horizontal coordinate and z the vertical one. For the LF part a vector field representation is chosen to better visualize the organized tumble motion.
In the text

Fig. 12 Scatter plots of the total flame surface value computed at −10 CAD as a function of the high frequency kinetic energy spatially averaged on a 2D disk of radius 3 mm around the spark plug. For each point the corresponding cycle number is indicated in the figure. The Pearson product-moment correlation coefficients of the two variables are also shown in the title. Upper figure: EC is computed considering the two velocity components within the central tumble plane. Lower figure: the third velocity component (normal to the central tumble plane) is also added in the EMD analysis.

In the text

	Fig. 13 Γ3p norm iso-surface at 0.5 coloured by its Y direction magnitude. The local direction is also shown over the iso-surface by Γ3p vector field. Cycle 17 and cycle 20 at three different CAD (−65, −30 and −22) are considered.
In the text

	Fig. 14 Γ3p norm on the spark mid-plane and 3D vector field inside its 3D iso-surface at 0.6. Cycle 17 and cycle 20 at CAD −22 are considered.
In the text

Fig. 15 Upper/middle figure: evolution of the local tumble ratio TR, at −55 CAD and −30 CAD for cycles 17 and 20, along the TCL identified with its Y coordinate. Lower figure: comparison between the local mean TR around the spark plug (SP) with the global tumble ratio for cycle 17 and 20 along the compression stroke between −95 and −20 CAD. The mean operation is made considering an enclosing volume of ±3 mm in Y direction.

In the text

	Fig. A.1 Upper figure: 1D plot comparison between experimental spray penetration (black dots) and numerical penetration issuing the 7 holes injector. Lower figure: 2D colormap comparison between experimental spray Shadowgraph and spray projected area 1 ms after the 1st droplet gets out of the nozzle.
In the text

	Fig. B.1 Comparison between the original (upper figures) and the accelerated EMD methodology (lower figure) for cycle 2 and CAD = −30.
In the text