Issue 
Sci. Tech. Energ. Transition
Volume 77, 2022
Dossier LES4ECE’21: LES for Energy Conversion in Electric and Combustion Engines, 2021



Article Number  20  
Number of page(s)  16  
DOI  https://doi.org/10.2516/stet/2022017  
Published online  25 October 2022 
Regular Article
ECFMLES modeling with AMR for the CCV prediction and analysis in leanburn engines
^{1}
IFP Energies nouvelles/Institut Carnot IFPEN Transports Energie, 1 et 4 avenue de BoisPréau, 92852 RueilMalmaison, France
^{2}
TOYOTA GAZOO Racing Europe GmbH, Chassis & Powertrain Development, Toyota Allee 7, 50858 Köln, Germany
^{*} Corresponding author: giampaolo.maio@ifpen.fr
Received:
26
November
2021
Accepted:
26
July
2022
A LargeEddy Simulation (LES) modeling framework, dedicated to ultralean sparkignition 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 gasexchange 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 cycletocycle 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 incylinder aerodynamic variations. R3P4. Finally two postprocessing tools, Empirical Mode Decomposition (EMD) and Γ_{3p} function, are proposed and combined to analyse for the first time the aerodynamic tumblebased incylinder 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.
Key words: Spark ignition engines / Combustion variability / Empirical Mode Decomposition
© The Author(s), published by EDP Sciences, 2022
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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), Plugin Hybrid Electric Vehicle (PHEV) allowing to reach ambitious goals in the fields of emission reduction in the next three decades [1]. The use of UltraLean Technology is a solution to increase combustion efficiency and to tackle pollutant formation in the design of modern sparkignition 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 zeroemission fuels such as hydrogen [3]. One of the main technological issues that affects ultralean engines is its capability to achieve stable operating conditions. Ultralean means here being close to the flame lean extinction limit. In order to enhance flame propagation, leanburn 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 highfidelity 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 phaseaveraged flow. On the contrary, the LargeEddy 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 sparkignition 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 leanburn TGRE/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 gasexchange conditions and for the reactive configuration with only slight differences in the geometry between cold and hot configurations. Some modelling work for these ultralean 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 postprocessing tools dedicated to complex turbulent reactive flows. After presenting the numerical modeling and setup, the results of the study are organised as follows:

Firstly, the LES numerical prediction of the flow field is validated in cold gasexchange conditions by comparing the velocity field on the tumble midplane with the PIV (Particle Imaging Velocimetry) measurements.

Secondly, the reactive configuration is investigated using ECFMLES (Extended Coherent Flame Model for LargeEddy Simulation) including the recently introduced effect of stretch on the laminar flame burning velocity [9]. Here, ECFMLES 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 postprocessing 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 setup
2.1 Engine configuration
The engine configuration investigated in this study is the singlecylinder research engine jointly developed by TGRE (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].
TGRE singlecylinder engine specifications.
2.2 Physical and LES numerical models
The version 3.0 of the commercial CFD code CONVERGE [13] is used in all computations. A secondorder central difference scheme is used for spatial discretization and secondorder Crank–Nicolson for time advancement. The Dynamic Smagorinsky model [14] is retained to close the subgrid scale stresses tensor. Spray injection is modeled with a Lagrangian formalism based on parcels [13]. Primary breakup 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 breakup 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_{6767}H_{11.767}O_{0.07}) already proposed for the same configuration in ref. [9]. The injection modeling setup 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 ECFMLES for flame propagation [17], while spark plug ignition is described by ISSIMLES (Imposed Stretch Spark Ignition Model) [18].
The ECFM model used is a modified version of the reference version, called speciesbased ECFM [19]. The progress variable is not transported but directly calculated as:(1)where is the filtered progress variable, the filtered species massfraction with superscript u for “unburned” and b for “burned”. The filtered species massfraction 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 subgridscale 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 autoignition 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 autoignition combustion mode contribution, is evaluated as:(4)where is the burned gases composition following [20], the postoxidation 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 isosurfaces 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 subgridscale turbulent transport, the subgridscale 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 cutcell 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 sparkplug, 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 nonnecessary 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 Progressvariable field on the tumble midplane overimposed 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 incylinder trapped mass and target Air/Fuel ratio estimated with the GTPower 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 Gasexchange cold flow
In a first phase, LES numerical modeling is validated for the prediction of the incylinder flow field. Indeed, the correct prediction of the incylinder 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 midplane. 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 incylinder 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 phaseaveraged mean velocity field on the tumble midplane 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 centreline (−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 gasexchange 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 multicycle traces it is observable that, despite the numerical pressure results show a high cycletocycle 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, autoignition, spark energy variation, etc.). Figure 5 shows also the evolution of the incylinder 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 multicycle incylinder 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 incylinder 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 incylinder 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 Multicycle 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. 
In the absence of autoignition 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 isosurface of the progress variable for a fastburning cycle (17) and for a slowburning 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 productmoment 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 isosurface 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 sparkplug 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 (ECHFSP), detailed in the following section.

The velocity field in the X direction (UX).

The mean Equivalence Ratio around the Spark Plug (ERSP).
For SparkPlug (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 (ECHFSP 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 incylinder 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 ECHFSP, 0.49 for UX and 0.34 for ERSP. 
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 cutplane 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 midplane. 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 postprocess 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 datadriven method for analyzing the nonlinear and nonstationary 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 (NAMEMD) was thus developed [35]. Another improvement was made in Adaptiveprojection intrinsically transformed MEMD (APITMEMD) [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 NAAPITMEMD. It was applied to a velocity field constructed with a largescale organised synthetic vortex and relatively smallscale isotropic homogeneous turbulence. Subsequently, the methodology was also applied to highspeed 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 largescale organized motion based only on the instantaneous velocity field, following:(11)with U_{LF} the largescale motion of lowfrequency and U_{HF} the turbulent fluctuations of highfrequency.
In [10], the EMD method has been adapted to the analysis of PIV fields in complex, timedependent 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 multicycle 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 incylinder 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 cyclebased 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 LowFrequency (LF) part which represents the largescale motion. The sum of the remaining IMFs is considered as the HighFrequency (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:

Apply the Ensemble NAAPITMEMD to each row (or column) of the bidimensional instantaneous velocity field U(x, y), i.e., decompose along one dimension. Sum the last IMF and the residual as “pseudo LF”.

Apply the Ensemble NAAPITMEMD 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.

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 speedup 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 incylinder 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 incylinder flow field from 3D LES in order to separate the largescale 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 sparkplug) 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 wellorganized tumble motion is still visible, while for cycle 17 a less wellstructured tumble is found. The tumble breakdown leads to an increase in the level of turbulent fluctuations near the sparkplug, 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 productmoment correlation coefficients between the spatial average of the high frequency energy around the spark plug at the sparktiming 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 nonorganized 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 productmoment 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 productmoment 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 Γ_{3} function is hereafter applied to the flow field, giving access to a 3D visualisation of the tumble evolution before ignition.
4.3 Γ_{3} function
The tumble flow is the main large scale rotational motion in the combustion chamber of a sparkignition engine. Algorithms called Γ_{1} and Γ_{2} 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 Γ_{2} criterion is Galileaninvariant because local convection is taken into account. A comparison between these two criteria was performed in [44] and showed that Γ_{1} was much more appropriate for locating vortex centres while Γ_{2} could help identify larger vortex regions. The Γ_{1} criterion has often been applied to incylinder flows, although using only 2D instantaneous velocity fields to detect and track the location of the tumble centre (see for example in [10]).
For threedimensional flows like the tumble flow, Γ_{1} is limited due to its 2D nature and thus has been extended in 3D. A very simple extension is to compute Γ_{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 Γ_{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 Γ in 3D was realized in [47], named Γ_{3} and was applied on external flows to identify vortical structures in 3D volume. Γ_{3} was then applied to the incylinder flow as presented in [11] with slight modifications: Γ_{3} 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.
Γ_{3} is first evaluated inside the combustion chamber using the following expression:(13)with φ the crank angle, V a subvolume centred at the point x and r(x_{0}) the distance vector points from x to x_{0}. (x_{0}, φ) is the local velocity u(x_{0}, φ) at x_{0} without the contribution from local convection 〈u〉_{V} (x, φ):(14)〈u〉_{V} (x, φ) is the average velocity inside V.
The formulation of Γ_{3} allows a reliable estimation of the rotation vector at point x by a simple normalization of Γ_{3}:(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 nonexistence 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 (Yaxis 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 Γ_{3} 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 Yaxis and this behavior varies from cycle to cycle.
Fig. 13 Γ_{3p} norm isosurface at 0.5 coloured by its Y direction magnitude. The local direction is also shown over the isosurface 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 midplane and 3D vector field inside its 3D isosurface 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 isosurface Γ_{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 Yaxis of cycle 20 at −22 CAD, unlike cycle 17, and its location on the central midplane is close to the sparkplug.
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 ECFMLES 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 ultralean spark ignition TGREconfiguration. 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 cycletocycle 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 incylinder 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 incylinder aerodynamic analysis, were applied to postprocess 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 sparkplug, 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 midplane 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 postprocessing 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 TGRE (TOYOTA GAZOO Racing Europe GmbH) for sharing the experimental configuration and measurements.
References
 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. [Google Scholar]
 Nagasawa T., Okura Y., Yamada R., Sato S., Kosaka H., Yokomori T., Iida N. (2021) Thermal efficiency improvement of superlean burn spark ignition engine by stratified water insulation on piston top surface, Int. J. Engine Res. 22, 5, 1421–1439. [CrossRef] [Google Scholar]
 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 leanburn natural gas SI engine, Int. J. Hydrogen Energy 45, 35, 17935–17952. [CrossRef] [Google Scholar]
 Benoit O., Truffin K., Jay S., van Oijen J., Drouvin Y., Kayashima T., Adomeit P., Angelberger C. (2022) Development of a LargeEddy Simulation methodology for the analysis of cycletocycle combustion variability of a lean burn engine, Flow Turbul. Combust. 108, 2, 559–598. [CrossRef] [Google Scholar]
 Enaux B., Granet V., Vermorel O., Lacour C., Pera C., Angelberger C., Poinsot T. (2011) LES study of cycletocycle variations in a spark ignition engine, Proc. Combust. Inst. 33, 2, 3115–3122. [Google Scholar]
 Truffin K., Angelberger C., Richard S., Pera C. (2015) Using largeeddy simulation and multivariate analysis to understand the sources of combustion cyclic variability in a sparkignition engine, Combust. Flame 162, 12, 4371–4390. [CrossRef] [Google Scholar]
 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 multicycle LES, Combust. Flame 156, 8, 1525–1541. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 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 largeEddy simulation. Technical report, SAE Technical Paper. [Google Scholar]
 Sadeghi M., Truffin K., Peterson B., Böhm B., Jay S. (2021) Development and application of Bivariate 2DEMD for the analysis of instantaneous flow structures and cycletocycle variations of incylinder flow, Flow Turbul. Combust. 106, 1, 231–259. [CrossRef] [Google Scholar]
 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. [Google Scholar]
 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 SymposiumAutomotive and Engine Technology, Springer Vieweg, Wiesbaden, pp. 205–224. [CrossRef] [Google Scholar]
 Senecal P. (2017) Converge (v2.4), Convergent Science Inc., Madison, WI. [Google Scholar]
 Germano M., Piomelli U., Moin P., Cabot W.H. (1991) A dynamic subgridscale Eddy viscosity model, Phys. Fluids A Fluid Dyn. 3, 7, 1760–1765. [CrossRef] [Google Scholar]
 Senecal P.K., Pomraning E., Richards K.J., Briggs T.E., Choi C.Y., McDavid R.M., Patterson M.A. (2003) Multidimensional modeling of directinjection diesel spray liquid length and flame liftoff length using CFD and parallel detailed chemistry, in: SAE Transactions, pp. 1331–1351. [Google Scholar]
 Frössling N. (1938) Gerlands beitr, Geophysics 52, 170–175. [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Colin O., Truffin K. (2011) A spark ignition model for large eddy simulation based on an FSD transport equation (ISSIMLES), Proc. Combust. Inst. 33, 2, 3097–3104. [CrossRef] [Google Scholar]
 Colin O., Chevillard S., Bohbot J., Senecal P., Pomraning E., Wang M. (2018) Development of a SpeciesBased Extended Coherent Flamelet Model (SBECFM) for Gasoline Direct Injection engine (GDI) simulations, in: Internal Combustion Engine Division Fall Technical Conference, Vol. 51999, American Society of Mechanical Engineers, p. V002T06A016. [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Hawkes E.R., Cant R.S. (2000) A flame surface density approach to largeEddy simulation of premixed turbulent combustion, Proc. Combust. Inst. 28, 1, 51–58. [CrossRef] [Google Scholar]
 Bougrine S., Richard S., Colin O., Veynante D. (2014) Fuel composition effects on flame stretch in turbulent premixed combustion: Numerical analysis of flamevortex interaction and formulation of a new efficiency function, Flow Turbul. Combust. 93, 2, 259–281. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Bounaceur R., Herbinet O., Fournet R., Glaude P.A., BattinLeclerc 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. [Google Scholar]
 Poinsot T., Veynante D. (2005) Theoretical and numerical combustion, RT Edwards Inc. [Google Scholar]
 Granet V., Vermorel O., Lacour C., Enaux B., Dugué V., Poinsot T. (2012) LargeEddy simulation and experimental study of cycletocycle variations of stable and unstable operating points in a spark ignition engine, Combust. Flame 159, 4, 1562–1575. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Matekunas F.A. (1983) Modes and measures of cyclic combustion variability, in: SAE Transactions, pp. 1139–1156. [Google Scholar]
 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 nonstationary time series analysis, Proc. R. Soc. Lond. Ser. A: Math. Phys. Eng. Sci. 454, 1971, 903–995. [CrossRef] [MathSciNet] [Google Scholar]
 Tanaka T., Mandic D.P. (2007) Complex empirical mode decomposition, IEEE Signal Proc. Lett. 14, 2, 101–104. [CrossRef] [Google Scholar]
 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 ProcessingICASSP’07, Vol. 3, IEEE, p. III–1009. [Google Scholar]
 Rilling G., Flandrin P., Goncalves P., Lilly J.M. (2007) Bivariate empirical mode decomposition, IEEE Signal Proc. Lett. 14, 12, 936–939. [CrossRef] [Google Scholar]
 Ur Rehman N., Mandic D.P. (2009) Empirical mode decomposition for trivariate signals, IEEE Trans. Signal Proc. 58, 3, 1059–1068. [Google Scholar]
 Ur Rehman N., Mandic D.P. (2011) Filter bank property of multivariate empirical mode decomposition, IEEE Trans. Signal Proc. 59, 5, 2421–2426. [CrossRef] [Google Scholar]
 Ur Rehman N., Park C., Huang N.E., Mandic D.P. (2013) EMD Via MEMD: Multivariate noiseaided computation of standard EMD, Adv. Adapt. Data Anal. 5, 2, 1350007. [CrossRef] [MathSciNet] [Google Scholar]
 Hemakom A., Goverdovsky V., Looney D., Mandic D.P. (2016) Adaptiveprojection intrinsically transformed multivariate empirical mode decomposition in cooperative brain–computer interface applications, Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci. 374, 2065, 20150199. [CrossRef] [PubMed] [Google Scholar]
 Sadeghi M., Foucher F., AbedMeraim K., MounamRousselle C. (2019) Bivariate 2D empirical mode decomposition for analyzing instantaneous turbulent velocity field in unsteady flows, Exp. Fluids 60, 8, 1–26. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Voisine M., Thomas L., Borée J., Rey P. (2011) Spatiotemporal structure and cycle to cycle variations of an incylinder tumbling flow, Exp. Fluids 50, 5, 1393–1407. [CrossRef] [Google Scholar]
 Abraham P., Liu K., Haworth D., Reuss D., Sick V. (2014) Evaluating LargeEddy Simulation (LES) and highspeed Particle Image Velocimetry (PIV) with phaseinvariant Proper Orthogonal Decomposition (POD), Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 69, 1, 41–59. [CrossRef] [Google Scholar]
 Cao J., Ma Z., Li X., Xu M. (2019) 3d proper orthogonal decomposition analysis of engine incylinder velocity fields, Meas. Sci. Technol. 30, 8, 085304. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Graftieaux L., Michard 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. [CrossRef] [Google Scholar]
 Bücker I., Karhoff D.C., Klaas M., Schröder W. (2012) Stereoscopic multiplanar piv measurements of incylinder tumbling flow, Exp. Fluids 53, 6, 1993–2009. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
 Nicollet F. (2019) Analysis of cyclic phenomena in a gasoline direct injection engine of flow and mixture formation using LargeEddy Simulation and highspeed particle image velocimetry. [Google Scholar]
 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. [CrossRef] [Google Scholar]
 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. [CrossRef] [Google Scholar]
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 Progressvariable field on the tumble midplane overimposed 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 phaseaveraged mean velocity field on the tumble midplane 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 centreline (−30, −10, 10 and 30 mm). 

In the text 
Fig. 5 Upper figure: LES multicycle incylinder 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 Multicycle 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. 

In the text 
Fig. 7 Flame kernel development identified by the progress variable isosurface 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 ECHFSP, 0.49 for UX and 0.34 for ERSP. 

In the text 
Fig. 10 Equivalence ratio distribution at −22 CAD on the tumble and spark midplane. 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 productmoment 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 isosurface at 0.5 coloured by its Y direction magnitude. The local direction is also shown over the isosurface 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 midplane and 3D vector field inside its 3D isosurface 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.