ECFM-LES modeling with AMR for the CCV
prediction and analysis in lean-burn engines

Giampaolo Maio, Zhihao Ding, Karine Truffin, Olivier Colin, Olivier Benoit,

Stéphane Jay

► To cite this version:

Giampaolo Maio, Zhihao Ding, Karine Truffin, Olivier Colin, Olivier Benoit, et al. ECFM-LES modeling with AMR for the CCV prediction and analysis in lean-burn engines. Science and Technology for Energy Transition, 2022, 77 (20), pp.16. 10.2516/stet/2022017. hal-03867855

HAL Id: hal-03867855

https://ifp.hal.science/hal-03867855

Submitted on 23 Nov 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ECFM-LES modeling with AMR for the CCV prediction and analysis in lean-burn engines

Giampaolo Maio<fnref n="1" />*, Zhihao Ding<fnref n="1" />, Karine Truffin<fnref n="1" />, Olivier Colin<fnref n="1" /> ID, Olivier Benoit<fnref n="2" /> ID, and Stéphane Jay<fnref n="1" /> ID

Received: 26 November 2021 / Accepted: 26 July 2022

**Abstract.** A Large-Eddy Simulation (LES) modeling framework, dedicated to ultra-lean spark-ignition engines, is proposed and validated in the present work. A direct injection research engine is retained as benchmark configuration. The LES model is initially validated using the cold gas-exchange conditions by comparing numerical results with PIV (Particle Imaging Velocimetry) experimental data. Then, the fired configuration is investigated, combining ECFM (Extended Coherent Flame Model) turbulent combustion model with Adaptive Mesh Refinement (AMR). The capability of the model to reproduce experimental pressure envelope and cycle-to-cycle variability is assessed. Within the major scope of the work, a particular focus on the Combustion Cyclic Variability (CCV) is made correlating them with the variability encountered in the in-cylinder aerodynamic variations. R3P4. Finally two post-processing tools, Empirical Mode Decomposition (EMD) and $\Gamma_{3p}$ function, are proposed and combined to analyse for the first time the aerodynamic tumble-based in-cylinder velocity field. Both tools make it possible to get deeply into the insight and visualization of the flow field and to understand the links between its cyclic variability and the combustion cyclic variability.

**Keywords:** Spark ignition engines, Combustion variability, Empirical Mode Decomposition.

## 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

<footnote n="1">¹IFP Energies nouvelles/Institut Carnot IFPEN Transports Energie, 1 et 4 avenue de Bois-Préau, 92852 Rueil-Malmaison, France</footnote>

<footnote n="2">²TOYOTA GAZOO Racing Europe GmbH, Chassis & Powertrain Development, Toyota Allee 7, 50858 Köln, Germany</footnote>

<footnote>* Corresponding author: giampaolo.maio@ifpen.fr</footnote>
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 (Γ₃ₚ)[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

<table><caption>**Table 1.** TGR-E single-cylinder engine specifications.</caption><tr><td>Bore/stroke</td><td>85.5/101.5 mm</td></tr><tr><td>Displacement</td><td>583 cm³</td></tr><tr><td>Geometrical CR</td><td>13.4</td></tr><tr><td>Ignition System Energy</td><td>Max 500 mJ</td></tr><tr><td>Fuel Direct Injection</td><td>7 holes</td></tr><tr><td>Injection pressure</td><td>Max. 500 bar</td></tr><tr><td>Fuel's octane quality</td><td>100 RON</td></tr></table>

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₁₁.₇₆₇O₀.₀₇) 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:

$$ \tilde{c} = \frac{\sum_{i} \tilde{Y}_{i}^{b}}{\sum_{i} \tilde{Y}_{i}^{u} + \sum_{i} \tilde{Y}_{i}^{b}}, \quad (1) $$

where $\tilde{c}$ is the filtered progress variable, $\tilde{Y}_i$ 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:
$$ \frac{\partial \bar{\rho} \tilde{Y}_i^X}{\partial t} + \nabla \cdot (\bar{\rho} \tilde{u}_i^X) = \nabla \cdot (\bar{\rho}(v + v_t) \nabla \tilde{Y}_i^X) + \bar{\rho} \tilde{\omega}_i^X, \quad (2) $$

where $\cdot$ and $\sim$ denote Reynolds and Favre filter respectively. $\bar{\rho}$ is the filtered density and $X$ represents either $u$ or $b$. $v$ is the molecular kinematic viscosity and $\hat{v}_t$ is the subgrid-scale kinematic viscosity estimated at the combustion filter size $\Delta = n_{res} \Delta_x$. $\Delta_x$ is a typical cell size and $n_{res}$ is set to 8 in the present study [4, 7]. The source term $\tilde{\omega}_i^X$ 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 $\tilde{\omega}_i^u$ is computed as:

$$ \tilde{\omega}_i^u = - \frac{\tilde{Y}_i^u}{1 - \tilde{c}} \left( \frac{\rho_u}{\bar{\rho}} S_l \bar{\Sigma} + \tilde{\omega}_c^{\text{ign}} \right), \quad (3) $$

where $\rho_u$ is the density of unburned gases, $S_l$ the laminar flame speed and $\bar{\Sigma}$ the Flame Surface Density (FSD).

For burned species, the source term $\tilde{\omega}_i^b$, without considering the auto-ignition combustion mode contribution, is evaluated as:

$$ \tilde{\omega}_i^b = Y_i^{b*} \left( \frac{\rho_u}{\bar{\rho}} S_l \bar{\Sigma} + \tilde{\omega}_c^{\text{ign}} \right) + (1 - \tilde{c}) \tilde{\omega}_i^{\text{post}}, \quad (4) $$

where $Y_i^{b*}$ is the burned gases composition following [20], $\tilde{\omega}_i^{\text{post}}$ 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 $\bar{\Sigma}$ and the reaction rate $\tilde{\omega}_c^{\text{ign}}$. 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 $\bar{\Sigma} = |\nabla c|$ is numerically unstable [17] and thus a modified FSD $\bar{\Sigma}_{\tilde{c}} = \bar{\Sigma} + \nabla \cdot ((\bar{c} - \tilde{c})\bar{N})$ is proposed, $\bar{N} = -\nabla \tilde{c} / |\nabla c|$ 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 $\bar{N}$ is replaced by the average normal direction of the flame front defined as follows:

$$ \langle \bar{N} \rangle_s = - \frac{\nabla \tilde{c}}{|\nabla c|} = \frac{1}{\Xi_{\text{sgs}}} \bar{N}. \quad (5) $$

Equation (5) takes into account the subgrid scale fluctuations of the vector orientation. The subgrid scale wrinkling $\Xi_{\text{sgs}}$ is defined following:

$$ \Xi_{\text{sgs}} = \frac{\bar{\Sigma}_{\tilde{c}}}{\bar{\Sigma}_{\tilde{c}}^{\text{lam}}} \quad (6) $$

with $\bar{\Sigma}_{\tilde{c}}^{\text{lam}} = |\nabla \tilde{C}| + (\bar{c} - \tilde{C})\nabla \cdot \bar{N}$ the laminar part of the FSD [17]. The factor $1/\Xi_{\text{sgs}}$ makes the laminar contribution disappear when the flow becomes very turbulent.

The transport equation for the filtered flame surface density $\bar{\Sigma}_{\tilde{c}}$ is finally modelled as follows:

$$ \begin{aligned} \frac{\partial \bar{\Sigma}_{\tilde{c}}}{\partial t} = & -\underbrace{\nabla \cdot (\tilde{u} \bar{\Sigma}_{\tilde{c}})}_{T_{\text{res}}} - \underbrace{\nabla \cdot (S_d (\bar{N})_s \bar{\Sigma}_{\tilde{c}})}_{P} \\ & + \underbrace{S_d \nabla \cdot (\langle \bar{N} \rangle_s) \bar{\Sigma}_{\tilde{c}}}_{C_{\text{res}}} + \underbrace{(\nabla \cdot \tilde{u} - \bar{n} \tilde{n} : \nabla \tilde{u}) \bar{\Sigma}_{\tilde{c}}}_{S_{\text{res}}} + \underbrace{\nabla \cdot ((v + \hat{v}_t) \nabla \bar{\Sigma}_{\tilde{c}})}_{T_{\text{sgs}}} \\ & - \underbrace{S_l \frac{(\bar{\Sigma}_{\tilde{c}} - \bar{\Sigma}_{\text{lam}}) \bar{\Sigma}_{\tilde{c}}}{1 - \tilde{c}} + \frac{2}{3} S_l (1 + \bar{\tau}) \frac{1 - \tilde{c}}{\bar{c}} (\bar{\Sigma}_{\tilde{c}} - \hat{\Sigma}_{\text{lam}}) \bar{\Sigma}_{\tilde{c}}}_{C_{\text{sgs}}} \\ & + \underbrace{\alpha_{\text{CFM}} \Gamma \left( \frac{\hat{u}'}{S_l^0}, \hat{\Delta} \right) \frac{\hat{u}'}{\hat{\Delta}} \bar{\Sigma}_{\tilde{c}}}_{S_{\text{sgs}}}, \end{aligned} \quad (7) $$

where $S_d = \rho_u S_l / \bar{\rho}$ is the displacement speed of the flame and $\bar{\tau}$ the expansion factor. The contributions from the resolved part of the flame are the convection $T_{\text{res}}$, the laminar propagation $P$, the filtered curvature $C_{\text{res}}$ and the filtered tangential strain rate $S_{\text{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_{\text{res}}$ is used in the present study:

$$ S_{\text{res}} = (\delta_{ij} - n_{ij} n_{ij}) \frac{\partial \tilde{u}_i}{\partial x_j} \bar{\Sigma}_{\tilde{c}}, \quad (8) $$

where

$$ n_{ij} n_{ij} = (\langle \bar{N} \rangle_s)_i (\langle \bar{N} \rangle_s)_j - \frac{1}{3} \alpha_n \delta_{ij} \quad (9) $$

with $\alpha_n = 1 - |\langle \bar{n} \rangle_s|^2$. The term $\alpha_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_{\text{sgs}}$, $C_{\text{sgs}}$ and $S_{\text{sgs}}$ are respectively the subgrid-scale turbulent transport, the subgrid-scale curvature and strain rate terms. The last term $S_{\text{sgs}}$ includes a calibration parameter $\alpha_{\text{CFM}}$ and an efficiency function $\Gamma$ developed by Bougrine et al. [22] which models the interaction between the resolved eddies of scales smaller than $\Delta$ 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 $\alpha(\mathbf{x}, t)$ which is equals to 0 before ignition and gradually increases to unity when ignition is ended. The factor $\alpha(\mathbf{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:

$$ \begin{aligned} \frac{\partial \bar{\Sigma}_{\tilde{c}}}{\partial t} = & T_{\text{res}} + T_{\text{sgs}} + S_{\text{sgs}} + \alpha C_{\text{sgs}} + \alpha (C_{\text{res}} + S_{\text{res}}) \\ & - \nabla \cdot \left( \alpha S_d \langle \bar{N} \rangle_s \bar{\Sigma}_{\tilde{c}} \right) + (1 - \alpha) K_{\text{ign}} \bar{\Sigma}_{\tilde{c}} + \bar{\omega}_{\Sigma}^{\text{ign}}, \end{aligned} \quad (10) $$
<figure><img src="image_2.png" /><figcaption>**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.</figcaption></figure>

where $K_{\text{ign}} = 2r_b^{-1}(1+\tau)\Xi_{\text{sgs}}S_l$ is the ISSIM modeled laminar stretch that replaces the resolved contributions $C_{\text{res}} + S_{\text{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.

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

<figure><img src="image_3.png" /><figcaption>**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).</figcaption></figure>

during the combustion phase while it allows to significantly improve the flame front resolution.

## 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
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

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.

## 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 ($\sigma$) 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

<figure><img src="image_5.png" /><figcaption>**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).</figcaption></figure>

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.

### 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
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

**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.

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

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

**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 $P_{max}$ as a function of the Crank Angle (CA) at which $P_{max}$ occurs. Lower plot: $P_{max}$ 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.

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
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

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.

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.

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.

<figure><img src="image_13.png" /><figcaption>**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).</figcaption></figure>

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.

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 $\phi$ variation in the whole chamber is limited to ±0.08 around the mean value ($\phi$ = 0.55). Furthermore the two cycles show similar $\phi$ values in the spark plug vicinity.

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

<figure><img src="image_14.png" /><figcaption>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.</figcaption></figure>

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:

$$U = U_{LF} + U_{HF} \quad (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.
<figure><img src="image_15.png" /><figcaption>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.</figcaption></figure>

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.

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]:

$$EC_{HF} = 0.5(\overline{U_{xHF}^2} + \overline{U_{zHF}^2}) \quad (12)$$

$\overline{U_{xHF}^2}$ and $\overline{U_{zHF}^2}$ 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
<figure><img src="image_16.png" /><figcaption>**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.</figcaption></figure>

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_y_{HF}$ (normal to the tumble plane) does not alter the above conclusions.

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

### 4.3 $\Gamma_3$ function

The tumble flow is the main large scale rotational motion in the combustion chamber of a spark-ignition engine. Algorithms called $\Gamma_1$ and $\Gamma_2$ have been developed to identify such a vortical structure by Graftiaux et al. [43]. The criterion is not Galilean invariant but allows a simple and robust method to identify rotational structures, while the $\Gamma_2$ criterion is Galilean-invariant because local convection is taken into account. A comparison between these two

criteria was performed in [44] and showed that $\Gamma_1$ was much more appropriate for locating vortex centres while $\Gamma_2$ could help identify larger vortex regions. The $\Gamma_1$ 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, $\Gamma_1$ is limited due to its 2D nature and thus has been extended in 3D. A very simple extension is to compute $\Gamma_1$ 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 $|\Gamma_1|$ 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 $\Gamma_1$ in 3D was realized in [47], named $\Gamma_3$ and was applied on external flows to identify vortical structures in 3D volume. $\Gamma_3$ was then applied to the in-cylinder flow as presented in [11] with slight modifications: $\Gamma_3$ was first computed to
<figure><img src="image_17.png" /><figcaption>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.</figcaption></figure>

obtain the axis of rotation, and then recomputed using the quantities projected in a plane normal to the axis of rotation, denoted as $\Gamma_{3p}$. In this work, a modified version of $\Gamma_{3p}$ will be computed to identify the tumble structure in the compression stroke.

$\Gamma_3$ is first evaluated inside the combustion chamber using the following expression:

$$ \Gamma_3(\mathbf{x}, \varphi) = \frac{1}{V} \int_{\mathbf{x}_0} \frac{\mathbf{r}(\mathbf{x}_0) \times \hat{\mathbf{u}}(\mathbf{x}_0, \varphi)}{\|\mathbf{r}(\mathbf{x}_0)\| \|\hat{\mathbf{u}}(\mathbf{x}_0, \varphi)\|} d\mathbf{x}_0 \quad (13) $$

with $\varphi$ the crank angle, V a sub-volume centred at the point $\mathbf{x}$ and $\mathbf{r}(\mathbf{x}_0)$ the distance vector points from $\mathbf{x}$ to $\mathbf{x}_0$. $\hat{\mathbf{u}}(\mathbf{x}_0, \varphi)$ is the local velocity $\mathbf{u}(\mathbf{x}_0, \varphi)$ at $\mathbf{x}_0$ without the contribution from local convection $\langle \mathbf{u} \rangle_V(\mathbf{x}, \varphi)$:

$$ \hat{\mathbf{u}}(\mathbf{x}_0, \varphi) = \mathbf{u}(\mathbf{x}_0, \varphi) - \langle \mathbf{u} \rangle_V(\mathbf{x}, \varphi), \quad (14) $$

$\langle \mathbf{u} \rangle_V(\mathbf{x}, \varphi)$ is the average velocity inside V.

The formulation of $\Gamma_3$ allows a reliable estimation of the rotation vector at point $\mathbf{x}$ by a simple normalization of $\Gamma_3$:

$$ \mathbf{e}_{\Gamma_3}(\mathbf{x}, \varphi) = \frac{\Gamma_3(\mathbf{x}, \varphi)}{\|\Gamma_3(\mathbf{x}, \varphi)\|} \quad (15) $$

Using Equation (15), projections of the local distance and velocity vectors on the rotational plane are obtained as follows:

$$ \mathbf{r}_p(\mathbf{x}_o) = \mathbf{r}(\mathbf{x}_o) - \mathbf{e}_{\Gamma_3}(\mathbf{x}, \varphi) \cdot (\mathbf{r}(\mathbf{x}_o) \cdot \mathbf{e}_{\Gamma_3}(\mathbf{x}, \varphi)), \quad (16) $$

$$ \mathbf{u}_p(\mathbf{x}_o, \varphi) = \hat{\mathbf{u}}(\mathbf{x}_o, \varphi) - \mathbf{e}_{\Gamma_3}(\mathbf{x}, \varphi) \cdot (\hat{\mathbf{u}}(\mathbf{x}_o, \varphi) \cdot \mathbf{e}_{\Gamma_3}(\mathbf{x}, \varphi)), \quad (17) $$

where $\mathbf{r}_p(\mathbf{x}_o)$ and $\mathbf{u}_p(\mathbf{x}_o, \varphi)$ are the projected distance and velocity, respectively. The final formulation of $\Gamma_{3p}$ is then given by:

$$ \Gamma_{3p}(\mathbf{x}, \varphi) = \frac{1}{V} \int_{\mathbf{x}_0} \frac{\mathbf{r}_p(\mathbf{x}_0) \times \mathbf{u}_p(\mathbf{x}_0, \varphi)}{\|\mathbf{r}_p(\mathbf{x}_0)\| \|\mathbf{u}_p(\mathbf{x}_0, \varphi)\|} d\mathbf{x}_0. \quad (18) $$

Compared to the definition in [11], $\mathbf{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 $\Gamma_{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 $\Gamma_{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 $\Gamma_{3p}$ reaches its maximum. All the local centres form a Tumble Centre Line (TCL) around which the tumble rotates.

$\Gamma_{3p}$ does not provide the rotation intensity. For that purpose a local tumble ratio was defined in [11] along the TCL at each centre $\mathbf{x}_c$ using a local average vorticity. In this study we propose another formulation to quantify the angular velocity around the TCL:

$$ \mathbf{TR}(\mathbf{x}_c, \varphi) = \frac{1}{\omega_{\text{ref}}} \frac{\int_{\mathbf{x}_0} \mathbf{r}_p(\mathbf{x}_0) \times \rho(\mathbf{x}_0, \varphi) \mathbf{u}_p(\mathbf{x}_0, \varphi) d\mathbf{x}_0}{\int_{\mathbf{x}_0} \rho(\mathbf{x}_0, \varphi) \|\mathbf{r}_p(\mathbf{x}_0)\|^2 d\mathbf{x}_0} \quad (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 $\mathbf{x}_c$. The computed angular velocity is normalized by the engine rotational speed $\omega_{\text{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 $\Gamma_{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 $\Gamma_3$ results

$\Gamma_{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
<figure><img src="image_18.png" /><figcaption>Fig. 13. $\Gamma_{3p}$ norm iso-surface at 0.5 coloured by its Y direction magnitude. The local direction is also shown over the iso-surface by $\Gamma_{3p}$ vector field. Cycle 17 and cycle 20 at three different CAD ($-65, -30$ and $-22$) are considered.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 14. $\Gamma_{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.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**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.</figcaption></figure>

obtained results: Γ₃ₚ 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.

In agreement with the EMD analysis, it is also possible to notice that for cycle 20, the iso-surface |Γ₃ₚ| = 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 Γ₃ₚ, 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.

## 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 Γ₃ₚ 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, Γ₃ₚ 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

1. Yugo M., Gordillo V., Shafiei E., Megaritis A. (2021) A look into the life cycle assessment of passenger cars running on advanced fuels, in: *Proceedings of the SIA Powertrains & Power Electronics Conference, Paris, France*, pp. 9–10.

2. Nagasawa T., Okura Y., Yamada R., Sato S., Kosaka H., Yokomori T., Iida N. (2021) Thermal efficiency improvement of super-lean burn spark ignition engine by stratified water insulation on piston top surface, *Int. J. Engine Res.* **22**, 5, 1421–1439.

3. Wang J., Duan X., Liu Y., Wang W., Liu J., Lai M.-C., Li Y., Guo G. (2020) Numerical investigation of water injection quantity and water injection timing on the thermodynamics, combustion and emissions in a hydrogen enriched lean-burn natural gas SI engine, *Int. J. Hydrogen Energy* **45**, 35, 17935–17952.

4. Benoit O., Truffin K., Jay S., van Oijen J., Drouvin Y., Kayashima T., Adomeit P., Angelberger C. (2022) Development of a Large-Eddy Simulation methodology for the analysis of cycle-to-cycle combustion variability of a lean burn engine, *Flow Turbul. Combust.* **108**, 2, 559–598.

5. Eaux B., Granet V., Vermorel O., Lacour C., Pera C., Angelberger C., Poinsot T. (2011) LES study of cycle-to-cycle variations in a spark ignition engine, *Proc. Combust. Inst.* **33**, 2, 3115–3122.

6. Truffin K., Angelberger C., Richard S., Pera C. (2015) Using large-eddy simulation and multivariate analysis to understand the sources of combustion cyclic variability in a spark-ignition engine, *Combust. Flame* **162**, 12, 4371–4390.

7. Vermorel O., Richard S., Colin O., Angelberger C., Benkenida A., Veynante D. (2009) Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle LES, *Combust. Flame* **156**, 8, 1525–1541.

8. Duva B.C., Chance L.E., Toulson E. (2020) Dilution effect of different combustion residuals on laminar burning velocities and burned gas Markstein lengths of premixed methane/air mixtures at elevated temperature, *Fuel* **267**, 117153.

9. Benoit O., Luszcz P., Drouvin Y., Kayashima T., Adomeit P., Brunn A., Jay S., Truffin K., Angelberger C. (2019) *Study of ignition processes of a lean burn engine using large-Eddy simulation*. Technical report, SAE Technical Paper.

10. Sadeghi M., Truffin K., Peterson B., Böhm B., Jay S. (2021) Development and application of Bivariate 2D-EMD for the analysis of instantaneous flow structures and cycle-to-cycle variations of in-cylinder flow, *Flow Turbul. Combust.* **106**, 1, 231–259.

11. Buhl S., Gleiss F., Köhler M., Hartmann F., Messig D., Brücker C., Hasse C. (2017) A combined numerical and experimental study of the 3D tumble structure and piston boundary layer development during the intake stroke of a gasoline engine, *Flow Turbul. Combust.* **98**, 2, 579–600.

12. Luszcz P., Takeuchi K., Pfeilmaier P., Gerhardt M., Adomeit P., Brunn A., Kupiek C., Franzke B. (2018) Homogeneous lean burn engine combustion system development – concept study, in: *18th Stuttgart International Symposium-Automotive and Engine Technology*, Springer Vieweg, Wiesbaden, pp. 205–224.

13. Senecal P. (2017) *Converge (v2.4)*, Convergent Science Inc., Madison, WI.

14. Germano M., Piomelli U., Moin P., Cabot W.H. (1991) A dynamic subgrid-scale Eddy viscosity model, *Phys. Fluids A Fluid Dyn.* **3**, 7, 1760–1765.

15. Senecal P.K., Pomraning E., Richards K.J., Briggs T.E., Choi C.Y., McDavid R.M., Patterson M.A. (2003) Multi-dimensional modeling of direct-injection diesel spray liquid length and flame lift-off length using CFD and parallel detailed chemistry, in: *SAE Transactions*, pp. 1331–1351.

16. Frössling N. (1938) *Gerlands beitr. Geophysics* **52**, 170–175.

17. Richard S., Colin O., Vermorel O., Benkenida A., Angelberger C., Veynante D. (2007) Towards large Eddy simulation of combustion in spark ignition engines, *Proc. Combust. Inst.* **31**, 2, 3059–3066.

18. Colin O., Truffin K. (2011) A spark ignition model for large eddy simulation based on an FSD transport equation (ISSIM-LES), *Proc. Combust. Inst.* **33**, 2, 3097–3104.

19. Colin O., Chevillard S., Bohbot J., Senecal P., Pomraning E., Wang M. (2018) Development of a Species-Based Extended Coherent Flamelet Model (SB-ECFM) for Gasoline Direct Injection engine (GDI) simulations, in: *Internal Combustion Engine Division Fall Technical Conference*, Vol. 51999, American Society of Mechanical Engineers, p. V002T06A016.

20. Colin O., Benkenida A., Angelberger C. (2003) 3D modeling of mixing, ignition and combustion phenomena in highly stratified gasoline engines, *Oil Gas Sci. Technol. - Rev. IFP Energies nouvelles* **58**, 1, 47–62.

21. Hawkes E.R., Cant R.S. (2000) A flame surface density approach to large-Eddy simulation of premixed turbulent combustion, *Proc. Combust. Inst.* **28**, 1, 51–58.

22. Bougrine S., Richard S., Colin O., Veynante D. (2014) Fuel composition effects on flame stretch in turbulent premixed combustion: Numerical analysis of flame-vortex interaction
and formulation of a new efficiency function, *Flow Turbul. Combust.* **93**, 2, 259–281.

23 Robert A., Richard S., Colin O., Martinez L., De Francqueville L. (2015) LES prediction and analysis of knocking combustion in a spark ignition engine, *Proc. Combust. Inst.* **35**, 3, 2941–2948.

24 Bounaceur R., Herbinet O., Fournet R., Glaude P.-A., Battin-Leclerc F., Pires Da Cruz A., Yahyaoui M., Truffin K., Moreac G. (2010) *Modeling the laminar flame speed of natural gas and gasoline surrogates*. Technical Report, SAE Technical Paper.

25 Poinsot T., Veynante D. (2005) *Theoretical and numerical combustion*, RT Edwards Inc.

26 Granet V., Vermorel O., Lacour C., Enaux B., Dugué V., Poinsot T. (2012) Large-Eddy simulation and experimental study of cycle-to-cycle variations of stable and unstable operating points in a spark ignition engine, *Combust. Flame* **159**, 4, 1562–1575.

27 Ameen M.M., Yang X., Kuo T.-W., Som S. (2017) Parallel methodology to capture cyclic variability in motored engines, *Int. J. Engine Res.* **18**, 4, 366–377.

28 Matekunas F.A. (1983) Modes and measures of cyclic combustion variability, in: *SAE Transactions*, pp. 1139–1156.

29 Huang N.E., Shen Z., Long S.R., Wu M.C., Shih H.H., Zheng Q., Yen N.-C., Tung C.C., Liu H.H. (1998) The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis, *Proc. R. Soc. Lond. Ser. A: Math. Phys. Eng. Sci.* **454**, 1971, 903–995.

30 Tanaka T., Mandic D.P. (2007) Complex empirical mode decomposition, *IEEE Signal Proc. Lett.* **14**, 2, 101–104.

31 Altaf M.U.B., Gautama T., Tanaka T., Mandic D.P. (2007) Rotation invariant complex empirical mode decomposition, in: *2007 IEEE International Conference on Acoustics, Speech and Signal Processing-ICASSP'07*, Vol. 3, IEEE, p. III-1009.

32 Rilling G., Flandrin P., Goncalves P., Lilly J.M. (2007) Bivariate empirical mode decomposition, *IEEE Signal Proc. Lett.* **14**, 12, 936–939.

33 Ur Rehman N., Mandic D.P. (2009) Empirical mode decomposition for trivariate signals, *IEEE Trans. Signal Proc.* **58**, 3, 1059–1068.

34 Ur Rehman N., Mandic D.P. (2011) Filter bank property of multivariate empirical mode decomposition, *IEEE Trans. Signal Proc.* **59**, 5, 2421–2426.

35 Ur Rehman N., Park C., Huang N.E., Mandic D.P. (2013) EMD Via MEMD: Multivariate noise-aided computation of standard EMD, *Adv. Adapt. Data Anal.* **05**, 02, 1350007.

36 Hemakom A., Goverdovsky V., Looney D., Mandic D.P. (2016) Adaptive-projection intrinsically transformed multi-variate empirical mode decomposition in cooperative brain-computer interface applications, *Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci.* **374**, 2065, 20150199.

37 Sadeghi M., Foucher F., Abed-Meraim K., Mounam-Rous-selle C. (2019) Bivariate 2D empirical mode decomposition for analyzing instantaneous turbulent velocity field in unsteady flows, *Exp. Fluids* **60**, 8, 1–26.

38 Fogleman M., Lumley J., Rempfer D., Haworth D. (2004) Application of the proper orthogonal decomposition to datasets of internal combustion engine flows, *J. Turbul.* **5**, 1, 023.

39 Voisine M., Thomas L., Borée J., Rey P. (2011) Spatio-temporal structure and cycle to cycle variations of an in-cylinder tumbling flow, *Exp. Fluids* **50**, 5, 1393–1407.

40 Abraham P., Liu K., Haworth D., Reuss D., Sick V. (2014) Evaluating Large-Eddy Simulation (LES) and high-speed Particle Image Velocimetry (PIV) with phase-invariant Proper Orthogonal Decomposition (POD), *Oil Gas Sci. Technol.* - *Rev. IFP Energies nouvelles* **69**, 1, 41–59.

41 Cao J., Ma Z., Li X., Xu M. (2019) 3d proper orthogonal decomposition analysis of engine in-cylinder velocity fields, *Meas. Sci. Technol.* **30**, 8, 085304.

42 Rulli F., Fontanesi S., d'Adamo A., Berni F. (2021) A critical review of flow field analysis methods involving proper orthogonal decomposition and quadruple proper orthogonal decomposition for internal combustion engines, *Int. J. Engine Res.* **22**, 1, 222–242.

43 Graftieux L., Richard M., Grosjean N. (2001) Combining piv, pod and vortex identification algorithms for the study of unsteady turbulent swirling flows, *Meas. Sci. Technol.* **12**, 9, 1422.

44 Bücker I., Karhoff D.-C., Klaas M., Schröder W. (2012) Stereoscopic multi-planar piv measurements of in-cylinder tumbling flow, *Exp. Fluids* **53**, 6, 1993–2009.

45 Janas P., Wlokas I., Böhm B., Kempf A. (2017) On the evolution of the flow field in a spark ignition engine, *Flow Turbul. Combust.* **98**, 1, 237–264.

46 Nicollet F. (2019) Analysis of cyclic phenomena in a gasoline direct injection engine of flow and mixture formation using Large-Eddy Simulation and high-speed particle image velocimetry.

47 Gohlke M., Beaudoin J.-F., Amielh M., Anselmet F. (2008) Thorough analysis of vortical structures in the flow around a yawed bluff body, *J. Turbul.* **9**, N15.

48 Schmitt M., Hu R., Wright Y.M., Soltic P., Boulouchos K. (2015) Multiple cycle LES simulations of a direct injection natural gas engine, *Flow Turbul. Combust.* **95**, 4, 645–668.
# 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.

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>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.</figcaption></figure>

# 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.

<figure><img src="image_23.png" /><figcaption>Fig. B.1. Comparison between the original (upper figures) and the accelerated EMD methodology (lower figure) for cycle 2 and CAD = -30.</figcaption></figure>