Karine Truffin
email: karine.truffin@ifpen.fr

Christian Angelberger

Stéphane Richard

Cécile Pera

Using Large-Eddy Simulation and multivariate analysis to understand the sources of combustion cyclic variability in a spark-ignition engine

Keywords: Large-Eddy Simulation, Cyclic combustion variability, spark-ignition engine, turbulent combustion, multivariate analysis, correlation

come

Using large-eddy simulation and multivariate analysis to understand the sources of combustion cyclic variability in a spark-ignition engine

fueled with a premixture of gaseous propane and air. The unstable points are obtained from the reference by reducing significantly the equivalence ratio and by an important dilution by nitrogen respectively. A LES methodology is proposed and shown to be able to reproduce the experimental findings concerning phase-averaged mean and statistical variations around it of a number of key engine combustion parameters. The CCV and factors causing it are first illustrated by comparing typical slow and fast burning cycles in combination with simple correlation plots of major engine parameters, this allows qualitatively showing how local and global sources concur to generate CCV. In a second step, single parameter and multivariate regressions build from the LES results allow quantifying the relative importance of different local and global CCV sources. Finally, the comparison of the obtained findings as to the relative importance of major parameters on CCV are compared with qualitative summary from an extensive experimental survey by Ozdor et al. The presented LES results overall confirm

Introduction

The operation of spark-ignition engines (SIE) is characterised by a nonrepeatability of instantaneous combustion rate between different cycles at nominally identical engine operating conditions, commonly referred to as cycle-tocycle variation [START_REF] Ozdor | Cyclic variability in spark-ignition engines. a literature survey[END_REF][START_REF] Heywood | Internal combustion engine fundamentals[END_REF] or cyclic combustion variability (CCV). CCV inevitably appears over the whole engine operation range, due in particular to the unsteady, cyclic and turbulent nature of flow and combustion in SIE. It is common practice to consider that a CCV amplitude (measured in terms of standard deviation of Indicated Mean Effective Pressure (IMEP), normalised by its mean value) of less than around 5 % is acceptable [START_REF] Chen | Impact of ignition energy phasing and spark gap on combustion in a homogenous direct injection gasoline si engine near the egr limit[END_REF]. As long as the CCV amplitude is sufficiently small, engine simulation softwares can predict with sufficient accuracy combustion using statistical approaches as RANS (Reynolds-Averaged Navier-Stokes) that neglect cyclic variations and aim at reproducing a phase averaged, statistically most probable, cycle. For higher CCV amplitudes however, individual cycles will behave very differently from this statistically most probable cycle. As a result, predictions of fuel consumption or emissions over a number of cycles may differ quite substantially from the one based on the mean engine cycle. These differences increase with the CCV amplitude and may reach extreme levels in the case of misfires or extreme knocking cycles [START_REF] Robert | Les prediction and analysis of knocking combustion in a spark ignition engine[END_REF][START_REF] Robert | Les study of deflagration to detonation mechanisms in a downsized spark ignition engine[END_REF].

In a context of increasingly stringent constraints on fuel consumption, CO 2 production, and pollutant emissions from road transport, it becomes crucial to be able to predict and control individual engine cycles, and thus to address the occurrence and effects of CCV. Engine technologies as downsizing [START_REF] Rudloff | Analysis of preignition in highly charged si engines: Emphasis on the auto-ignition mode[END_REF][START_REF] Zaccardi | Development of specific tools for analysis and quantification of pre-ignition in a boosted si engine[END_REF], direct injection (DI) [START_REF] Adomeit | Analysis of cyclic fluctuations of charge motion and mixture formation in a disi engine in stratified operation[END_REF] or controlled auto-ignition (CAI) [START_REF] Thirouard | Nature of cai combustion and air/fuel ratio stratification effects[END_REF][START_REF] Adomeit | Operation strategies for controlled auto ignition gasoline engines[END_REF] are examples of technologies presently explored in order to reduce the CO 2 emissions from future SIE. Yet the occurrence under certain operating conditions of excessive CCV when implementing these technologies is one of the factors limiting their practical performance or range of operation. Being able to predict CCV in early design phases based on an improved knowledge of their sources and effects could effectively contribute exploiting the full potential of these promising SI technologies under real operation.

The understanding of how the complex combination of different sources leads to the occurrence of CCV for a specific engine design or mode of operation is still limited. CCV is indeed the result of a complex combination of different flow phenomena. These phenomena can be classified into two main categories, depending on the related spatial and temporal scales:

• Global phenomena [START_REF] Bedford | Measurements of cycle to cycle variability of the inlet flow of fuel injectors using lda[END_REF][START_REF] Persson | Cylinderto-cylinder and cycle-to-cycle variations at hcci operation with trapped residuals[END_REF], which are related to global operating characteristics: trapped mass, intake mass flow rate, tumble ratio, overall equivalence ratio, mean cylinder pressure, mean intake charge temperature, exhaust gas recirculation (EGR) rate, etc. Their spatial scales are of the order of some characteristic dimension of the engine, and their time scales of the order of the crank angle and up to the duration of an engine cycle;

• Local phenomena [START_REF] Mehrani | Modeling the effects of mixture composition on cyclic variablity[END_REF][START_REF] Poinsot | Theoretical and Numerical Combustion, 2nd Edition[END_REF][START_REF] Aleiferis | Flame chemiluminescence studies of cyclic variations and a/f ratio of the reacting mixture in a si engine[END_REF], which are related to local flow variables: temperature, pressure, mixture composition, turbulence, flow velocity, strain rate, etc. Their spatial scales range from some micrometers up to some millimeters, their time scales from some microseconds up to the crank angle.

Both scales can exhibit cycle-to-cycle variations and do strongly interact, making the understanding of CCV complex. This is especially true in the combustion chamber, where the non-linear response of combustion is a key contributor to CCV.

Experimental studies [START_REF] Bedford | Measurements of cycle to cycle variability of the inlet flow of fuel injectors using lda[END_REF][START_REF] Persson | Cylinderto-cylinder and cycle-to-cycle variations at hcci operation with trapped residuals[END_REF] certainly allow identifying CCV for a given engine concept, and they can help testing control strategies to limit its occurrence and impact. Aleiferis [START_REF] Aleiferis | Cyclic variations of initial flame kernel growth in a honda vtec-e lean-burn spark-ignition engine[END_REF] conducted research on a pent-roof single-cylinder port fuel injection engine fueled with a stratified mixture of iso-octane/air. He concluded that the flame area was highly correlated to the crank angle at which 5 % fuel mass fraction was burnt, while the flame volume was highly correlated to spark energy. In a recent paper [START_REF] Chen | Impact of ignition energy phasing and spark gap on combustion in a homogenous direct injection gasoline si engine near the egr limit[END_REF] the impact of the ignition system on the combustion in a direct-injection engine was studied, showing that larger spark gaps, longer spark life durations and multi-spark ignitions were beneficial to the extension of the dilution limit and to the combustion stability. Such experiments only come into play once a prototype has actually been built, in design phases where modifications to the basic concept are difficult to achieve because of cost and time constraints. Furthermore, the achievable understanding is limited by the impossibility to have access to all necessary thermodynamic and flow quantities on an instantaneous, cycle resolved basis. Nevertheless this is required to gain a basic understanding of the sources of CCV. The understanding gained from such studies is thus only partial and rarely valid beyond the specific studied case.

Large-Eddy Simulation (LES) is a 3D-CFD technique that has the potential to address the numerous coupled phenomena influencing engine combustion, and potentially gives access to any quantity needed to understand and characterize CCV. Its ability to deal with turbulent flows [START_REF] Banaeizadeh | Large-eddy simulations of turbulent flows in internal combustion engines[END_REF][START_REF] Huijnen | Study of turbulent flow structures of a practical-steady engine head flow using large-eddy simulations[END_REF][START_REF] Celik | Prediction of in-cylinder turbulence for ic engines[END_REF][START_REF] Haworth | Large eddy simulation of in-cylinder flows[END_REF] and reactive flows [START_REF] Tatschl | Les simulation of flame propagation in a direct-injection si-engine to identify the causes of cycle-to-cycle combustion variations[END_REF][START_REF] Goryntsev | Analysis of misfire processes in realistic direct injection spark ignition engine using multi-cycle large eddy simulation[END_REF][START_REF] Drake | Advanced gasoline engine development using optical diagnostics and numerical modeling[END_REF] in SI engines has been demonstrated. But only a few LES studies were dedicated to the prediction of CCV [START_REF] Granet | Largeeddy simulation and experimental study of cycle-to-cycle variations of stable and unstable operating points in a spark ignition engine[END_REF][START_REF] Enaux | Les study of cycle-to-cycle variations in a spark ignition engine[END_REF][START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle les[END_REF] using either the Thickened Flame (TFLES) [START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF], or the Extended Coherent Flame (ECFM-LES) [START_REF] Richard | Towards large eddy simulation of combustion in spark ignition engines[END_REF][START_REF] Hawkes | A flame surface density approach to large-eddy simulation of premixed turbulent combustion[END_REF][START_REF] Candel | Flame stretch and the balance equation for the flame surface area[END_REF] turbulent combustion models and the AVBP code1 [START_REF] Gourdain | High performance parallel computing of flows in complex geometries part 1: methods[END_REF][START_REF] Gourdain | High performance parallel computing of flows in complex geometries part 2: applications[END_REF]. Vermorel et al. [START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle les[END_REF] have applied the LES models developed by Richard et al. [START_REF] Richard | Towards large eddy simulation of combustion in spark ignition engines[END_REF] to explore the origins of CCV in an SIE fueled with gaseous propane using AVBP. They have demonstrated that computing 10 consecutive complete four-stroke engine cycles allowed reproducing experimental findings on CCV.

They also illustrated how analysing this LES could be used to identify sources of CCV. In the studied case they were found to be stochastic in nature and related to the turbulence of the intake flow, and its coupling to spark-ignition and flame propagation. Lacour & al. [33] have acquired the SGEmac experimental database dedicated to a detailed study of CCV in an optical access single cylinder SIE fueled with a homogeneous mixture of propane and air. Low CCV operating points were acquired in order to fully characterize the engine and allow model validation. Two operating points were then explored in order to characterize high CCV values resulting either from an important dilution by nitrogen of the fresh gas mixture (EGR emulation), or by reducing the fuel/air equivalence ratio. Available data comprise crank angle resolved measurements of pressures and temperatures at different locations in the engine set-up, as well as visualisations of velocities and combustion progress using optical diagnostics. This database has served for validating the prediction of CCV using LES.

In [START_REF] Enaux | Large eddy simulation of a motored single-cylinder piston engine: Numerical strategies and validation[END_REF], LES of the flow under motored (without combustion) conditions were performed with AVBP and were shown to yield an accurate reproduction of the flow field inside the combustion chamber and of the acoustics in the intake and exhaust ducts. In [START_REF] Granet | Largeeddy simulation and experimental study of cycle-to-cycle variations of stable and unstable operating points in a spark ignition engine[END_REF] and [START_REF] Enaux | Les study of cycle-to-cycle variations in a spark ignition engine[END_REF] the same numerical set-up around AVBP was used for LES of fired operations. The simulation of 25 consecutive full cycles of the stable and 50 cycles of an unstable operating point in [START_REF] Granet | Largeeddy simulation and experimental study of cycle-to-cycle variations of stable and unstable operating points in a spark ignition engine[END_REF] using the TFLES combustion model [START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF] demonstrated the ability of the employed LES approach to quantitatively predict the CCV levels observed experimentally. Qualitative analysis of the LES was performed to understand the causes for incomplete combustion which occurred in some cycles of the high-CCV case, but no quantitative analysis was proposed that could have allowed correlating CCV with global or local phenomena.

The present work aims i/ at complementing these studies of CCV with an alternative LES methodology based on AVBP, ii/ at both reproducing experimentally observed CCV in the SGEmac engine, and iii/ at proposing a systematic way to explore and quantify its sources. While the above mentioned studies were based on simulating the whole engine set-up with LES, the computational domain including the whole intake and exhaust lines between the intake and exhaust plena, the present work is based on an alternative LES methodol-ogy [START_REF] Pera | Large eddy simulation of a motored single-cylinder engine using system simulation to define boundary conditions: Methodology and validation[END_REF], exploiting system simulation of the intake and exhaust ducts to impose unsteady boundary conditions for the LES of the combustion chamber and its immediate neighborhood. Furthermore, premixed combustion is addressed using ECFM-LES [START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle les[END_REF][START_REF] Richard | Towards large eddy simulation of combustion in spark ignition engines[END_REF], a LES formulation of the widely used Coherent Flame Model.

The objective of the present work is threefold:

• apply the LES methodology developed and validated for motored operation in [START_REF] Pera | Large eddy simulation of a motored single-cylinder engine using system simulation to define boundary conditions: Methodology and validation[END_REF] to fired operation based on the ECFM-LES combustion model, and demonstrate its ability to reproduce experimentally observed CCV levels;

• acquire a basic understanding of the origins of CCV and of their impact, as a result of both global and local phenomena;

• propose a systematic method for analysing LES of CCV in order to identify major parameters affecting it, and providing a quantification of their relative importance.

The SGEmac single cylinder engine experiment

The SGEmac engine test bench [START_REF] Lacour | An experimental database dedicated to the study and modelling of cyclic variability in spark-ignition engines with les[END_REF] was specifically designed to allow a detailed validation of the ability of LES to predict CCV. It consists of a single cylinder spark-ignition engine equipped with a 4-valve pentroof combustion chamber and with a flat piston. The main engine characteristics are summarised in Table 1. The whole SGEmac set-up [START_REF] Granet | Largeeddy simulation and experimental study of cycle-to-cycle variations of stable and unstable operating points in a spark ignition engine[END_REF] includes: [START_REF] Ozdor | Cyclic variability in spark-ignition engines. a literature survey[END_REF] an intake plenum to which air is fed via a sonic nozzle in order to perfectly control mass flow; [START_REF] Heywood | Internal combustion engine fundamentals[END_REF] a mixing plenum into which gaseous propane is injected far upstream the combustion chamber for creating a perfect reactive pre-mixture; (3) the combustion chamber with multiple optical accesses; (4) an exhaust plenum in order to perfectly control and define boundary conditions. The volumes of both plena were chosen sufficiently large to damp acoustic perturbations. The characteristics of the three simulated operating points from the SGEmac database [START_REF] Lacour | An experimental database dedicated to the study and modelling of cyclic variability in spark-ignition engines with les[END_REF] are summarized in Table 2. The stab ref operating point is the stable reference point at stoichiometry with no external dilution. The unstable operating points unstab dil and unstab lean result from respectively diluting the fresh gases with additional N 2 , and air (by reducing the fuel/air equivalence ratio). Both of these unstable points were defined in order to keep the mean indicated effective pressure (IMEP) as close as possible to that of stab ref, resulting in necessary adaptations of the spark timing.

LES methodology

Numerical approach

The present LES were performed using the parallel AVBP solver [START_REF] Gourdain | High performance parallel computing of flows in complex geometries part 1: methods[END_REF][START_REF] Gourdain | High performance parallel computing of flows in complex geometries part 2: applications[END_REF] which is dedicated to the simulation of reactive unsteady flows and explicitly solves the compressible multi-species Navier-Stokes equations on unstructured meshes with a cell-vertex finite-volume formulation. A second order centered Lax-Wendroff convective scheme is used in this study, in combination with an explicit time-advancement in an attempt to achieve a reasonable compromise between stability and numerical dissipation [START_REF] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF].

Computational domain and boundary conditions

In order to avoid simulating the full engine set-up using LES, the methodology developed in [START_REF] Pera | Large eddy simulation of a motored single-cylinder engine using system simulation to define boundary conditions: Methodology and validation[END_REF] was adopted. First, a sketch of the full engine set-up is built and simulated using the IFP-Engine system simulation (or 1D CFD) library [START_REF] Malbec | Modelling turbochard spark-ignition engines: Towards predictive real time simulators[END_REF][START_REF] Lafossas | Application of a new 1d combustion model to gasoline transient engine operation[END_REF] within LMS-Imagine.Lab solver. From this system simulation, crank-angle resolved pressure and temperature evolutions at the locations corresponding to the boundaries of the 3D computational domain are extracted, and prescribed as boundary conditions for the LES. Species mass fractions of the gaseous mixture of air and propane are also imposed. Only one-way coupling between the system simulation and LES is hereby accounted for, so that no feedback from the LES to the system simulation is addressed. As the cyclic variability of the flow in the intake and exhaust ducts had been experimentally

shown to be very small [START_REF] Lacour | An experimental database dedicated to the study and modelling of cyclic variability in spark-ignition engines with les[END_REF], the same time evolutions of pressure and temperature were imposed for all simulated, consecutive LES cycles.

At the intake and exhaust boundaries, a multi-species extension [START_REF] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF] of the Navier-Stokes Characteristic Boundary Conditions (NSCBC) [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF] were imposed via a linear relaxation method (LRM) [START_REF] Selle | The actual impedance of non-reflecting boundary conditions: implications for the computation of resonators[END_REF]. This method allows separating waves entering and leaving the domain, thus minimizing spurious reflexions at the boundaries. The cut-off frequency of this low-pass filter is directly proportional to a relaxation coefficient that ensures that the imposed pressure relaxes towards the target time signal coming from system simulation. Finally a mixed inflow /outflow formulation was selected in order to handle the occurrence of flow reversals during certain phases of the engine cycle, a nominal inflow becoming an outflow and vice versa. To this purpose, the boundary conditions are switching from an inlet to an outlet formulation according to the instantaneous direction of the local flow velocity. Fig. 1(a) shows the extent of the computational domain used in the LES, comprising the combustion chamber and its immediate neighborhood. During fired operation, hot burnt products from the combustion chamber can be pushed back into the intake ducts and reach the inflow boundary if a too short portion of the intake duct is included in the LES. In order to avoid numerical problems resulting from the composition change at the boundary, 1D CFD results were used to estimate the backflow distance, the length of the intake duct included in the mesh was set to 70 mm. The outlet boundary conditions was set at the location of an experimental pressure sensor at the junction of the two exhaust ducts coming from the cylinder. No attempt was made to impose turbulent flow conditions at the inflow boundaries of the LES domain. As a consequence of the reciprocating nature of the engine, the flow inside the intake and exhaust ducts is strongly oscillating, even during the closed valve phases. Therefore, the turbulent flow inside the ducts cannot be considered as corresponding to a fully established duct flow, despite the high flow Reynolds numbers and the important duct lenghts. Any attempt to impose realistic turbulent flow patterns on inflow boundaries would therefore be a complex task, that would require a detailed study of the unsteady 3D flow in the engine ducts. In the present LES, turbulence in the ducts was therefore only a result from the shear generated by the wall boundary layers. Neglecting inflow turbulence appears all the more justified as the key flow feature responsible for turbulence during combustion is the strong flow detachment generated around the valves during intake that induces a tumbling flow, the compression of which by the piston ultimately leads to small scall turbulence. An estimation of the time scales of any turbulence entering the cylinder from the ducts shows that they are much shorter than the duration of the compression phase. Such turbulence is therefore dissipated very fast during intake and does not have any impact on combustion. The obtained accurate predictions of the tumbling flow and its compression in the present LES a posteriori justify this simplifying assumption. Convective wall heat transfer was modeled using a logarithmic law of the wall as described in [START_REF] Schmitt | Large-eddy simulation and experimental study of heat transfer, nitric oxide emissions and combustion instability in a swirled turbulent high-pressure burner[END_REF] using the Diwakar model [START_REF] Diwakar | Assessment of the ability of a multidimensional computer code to model combustion in a homogeneous-charge engine[END_REF], which was specifically dedicated to the prediction of wall heat flux in spark-ignition engines. Finally, the spark plug was not discretized in the present simulations. This appeared justified as the electrodes protrude very few into the combustion chamber. Attempts (not shown here) to account for the spark plug in the computational mesh demonstrated that this does not fundamentally affect conclusions drawn on the on-set and early development of combustion described in the present paper.

Moving mesh management

The engine cycle is divided into 41 mesh phases, each characterized by a constant connectivity, which changed from phase to phase to adapt to the important displacement of the piston and valves. Each mesh phase contains an initial and a final grid generated before the simulation using an external grid generator. Within a mesh phase, the position of each node is interpolated in time and space between its initial position in the initial mesh and its final position in the target mesh. The quality of the deforming tetrahedral mesh is controlled through mesh smoothing [START_REF] Moureau | Towards large eddy simulation in internal-combustion engines: simulation of a compressed tumble flow[END_REF]. In order to restrict the grid deformation while maintaining enough spatial resolution, a second-order interpolation is used [START_REF] Forsythe | Validation of a fluid-structure interaction model for a bileaflet mechanical heart valve[END_REF]. In the present computations, the grid deformation maintains an average spatial resolution close to 0.7 mm in the cylinder during the engine cycle, with down to 0.2 mm around the spark plug position and during combustion, as well as in the valve curtains. The smallest cell sizes of 0.04 mm are used for the residual valve lift of 0.35 mm. From the cylinder-head outwards into the ducts, meshes are gradually coarsened, largest cells sizes in the ducts reaching values of about 4 mm. The smallest grid is found at TDC and consists of 2 million tetrahedral cells, while the largest grid of 8.2 million tetrahedral cells is found at the intake BDC.

Sub-grid scale modelling

Unclosed sub-grid stresses in the momentum equations are modeled using the Sigma model [START_REF] Nicoud | Using singular values to build a subgrid-scale model for large eddy simulations[END_REF]. It is specifically designed so that turbulent viscosity vanishes in different canonical laminar flows (pure shear, pure rotation, isotropic expansion and contraction, stagnation point) for which no sub-grid scale (sgs) activity is expected. This model also properly propagates a vortex ring along the wall after the impingement of a pulsating jet [START_REF] Toda | Assessment of subgrid-scale models with a les-dedicated experimental database: pulsatile impinging jet in turbulent cross flow[END_REF], a phenomenon commonly encountered in piston engines.

The ECFM-LES turbulent premixed combustion model [START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle les[END_REF] is used to describe the flame front propagation. It combines two approaches: first, the CFM-LES model [START_REF] Richard | Towards large eddy simulation of combustion in spark ignition engines[END_REF] is an adaptation of the original CFM model [START_REF] Candel | Flame stretch and the balance equation for the flame surface area[END_REF][START_REF] Pope | The evolution of surfaces in turbulence[END_REF] to LES, based on a transport equation for the filtered flame surface density (FSD) Σ c . Second, ECFM [START_REF] Colin | Angelberger, 3d modeling of mixing, ignition and combustion phenomena in highly stratified gasoline engines[END_REF] is an extension of the CFM model to allow the description of the turbulent combustion of stratified mixtures under variable pressure and temperature conditions. Full details of the ECFM-LES formalism and of the transport equations for the filtered mass densities species ρ Y i , the fuel tracer Y T F , the mean and fresh gases enthalpy can be found in [START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle les[END_REF][START_REF] Richard | Towards large eddy simulation of combustion in spark ignition engines[END_REF] and are therefore not given here. For the present calculations, only major species are considered:

Fuel, O 2 , N 2 , CO 2 , H 2 O, CO and H 2 .
The composition of the burnt gases is obtained considering two reactions:

α[C x H y + (x + y 4 )O 2 → xCO 2 + y 2 H 2 O] (1 -α)[C x H y + x 2 O 2 → xCO + y 2 H 2 ] (1)
The parameter α is defined following:

if φ < 1 : α = 1 else α = 4(x+ y 4 ) φ -2x 2x + y ( 2
)
where φ is the equivalence ratio. The species reaction rate ωk and the heat release ωhs are then modeled following:

ωk = ρ u Y T F S L Σ c αν ′ k1 + (1 -α)ν ′ k2 (3) ωhs
= k ∆h 0 f,k ωk (4)
ν ′ k1 and ν ′ k2 are the coefficients of reactions 1 and 2 of species k in Eq. ( 1). The burnt gas enthalpy is computed from the average and fresh gases enthalpy [START_REF] Vermorel | Towards the understanding of cyclic variability in a spark ignited engine using multi-cycle les[END_REF] following

h b s = ( h s -(1 -c) h u s )/ c.
The superscripts u and b refer to fresh and burnt gases conditions, respectively. c is the Favre filtered progress variable defined using the fuel mass fraction according to c = 1 -Y u F u /Y T F . Y T F is the fuel mass fraction before the onset of combustion [START_REF] Colin | A spark ignition model for large eddy simulation based on a fsd transport equation (issim-les)[END_REF]. The burnt gas temperature and density are therefore deduced knowing the composition and enthalpy of the burnt gases. The post-oxidation in the burnt gases is omitted in the present study as the fuel/air mixture is nearly homogeneous and lean or diluted.

Spark-ignition is modeled using the ISSIM-LES model [START_REF] Colin | A spark ignition model for large eddy simulation based on a fsd transport equation (issim-les)[END_REF] which includes the description of the electrical circuit and 3D modeling of the flame kernel growth during ignition. An electrical circuit model, based on the description provided by AKTIM [START_REF] Duclos | Arc and kernel tracking ignition model for 3d sparkignition engine calculations[END_REF][START_REF] Richard | Development of a one-dimensional computational fluid dynamics modeling approach to predict cycle-to-cycle variability in spark-ignition engines based on physical understanding acquired from large-eddy simulation[END_REF] defines the electrical power P ign (t) and energy E ign (t) = t t spark P ign (t ′ )dt ′ transferred to the gas during the spark life (breakdown and glow phase), as well as the initial burnt gases profile and reaction rates. Throughout the ignition phase, the flame kernel evolution is directly controlled by the flame surface density equation, the source terms of which are modified during early ignition to accurately represent the flame surface growth. The burnt gas temperature is computed as described previously. ISSIM-LES thus accounts for the effect of mixture stratification, convection and turbulent stretch on the early flame kernel growth. As long as the flame radius r b remains smaller than the combustion filter size ∆ = 5∆ x where ∆ x is the cell size and c is smaller than one, the flame kernel is described at the SGS level using ISSIM-LES. The burnt gases r b is given by a transport equation to allow multi-sparks ignitions.

When r b approaches ∆, the progress variable c reaches unity and the transition between ISSIM-LES and ECFM-LES is progressively reached using a coefficient α = 0.5(1 + tanh(r b / ∆) -0.75/0.15) which is one during ignition and zero afterwards. This enables a smooth transition to the freely propagating turbulent flame. The impact of this criteria has not been examined here but it should be weak in this situation where the mesh size is small (0.2mm) so that transition occurs fast (a few degrees after spark-ignition).

The parameters of ISSIM-LES are set to reproduce the spark-ignition system used in the SGEmac experiments. It consists of a coil ignition system (Sagem type BAE 04). The spark gap is 1.2 mm, the total energy stored in the coil is sufficiently high (about 600 mJ) to prevent misfires. The secondary inductance is around 30 H and the secondary resistance around 10 kΩ. The spark duration can vary from cycle to cycle as it depends on the available energy in the secondary spark and on the energy transferred to the gas, the evolution of which is given by an ODE and incorporates the time evolution of the spark length with the convection [START_REF] Richard | Development of a one-dimensional computational fluid dynamics modeling approach to predict cycle-to-cycle variability in spark-ignition engines based on physical understanding acquired from large-eddy simulation[END_REF]. Overall, the spark life varies here between 500 and 600 µs, that is around 4 CAD.

In Eq. ( 3), the local laminar flame speed S L is obtained from a correlation of Metghalchi and Keck [START_REF] Metghalchi | Laminar burning velocity of propane-air mixtures at high temperature and pressure[END_REF][START_REF] Granet | Largeeddy simulation and experimental study of cycle-to-cycle variations of stable and unstable operating points in a spark ignition engine[END_REF] as a function of local thermodynamic conditions, and with a correction to account for the dilution rate [START_REF] Colin | Angelberger, 3d modeling of mixing, ignition and combustion phenomena in highly stratified gasoline engines[END_REF]. It takes the following form:

S L = S 0 L T u T 0 α p p 0 β (1 -2.1 X dil ) (5)
with α = 2.18 -0.8( φ -1), β = -0.16 + 0.22( φ -1), T 0 = 298 K, p 0 = 1 bar.

In the ECFM-LES and ISSIM-LES approaches fuel chemistry is not explicitly solved for when determining the heat released by the flame, but is rendered implicitly via the laminar flame speed and its effect on spark-ignition, flame surface creation and destruction, and heat release rate.

Computing time

The computation of a full engine cycle of 720 CAD took about 55 hours on a number of cores between 256 and 512 of a SGI Altix ICE 8200. The number of cores was adapted to the varying mesh size during the different mesh phases in order to optimise the exploitation of parallelism during the full engine cycle.

LES of the stable operating point stab ref

In order to validate the developed LES methodology, it was first applied to the simulation of twelve consecutive full engine cycles of the stable reference operating point stab ref.

The first two cycles were discarded (as they were in all other presented LES), as they are considered being polluted by the initial conditions and the results were compared with experimental findings obtained by acquiring 100 consecutive cycles.

The ability of the present methodology to predict the intake flow observed experimentally has already been reported in [START_REF] Pera | Large eddy simulation of a motored single-cylinder engine using system simulation to define boundary conditions: Methodology and validation[END_REF] and will not again be detailed here.

Fig. 2 illustrates a typical flame propagation based on an individual LES realisation of stab ref at four different instants after spark-ignition (SI). At SI + 5 CAD, the flame is not spherical anymore, and exhibits a preferential propagation towards the exhaust side (negative x) under the effect of the mean flow velocity at the spark resulting from the tumble flow created during the intake stroke. Although it interacts earlier with the cylinder head on the exhaust side, it nevertheless keeps propagating at later instants faster towards the exhaust than towards the intake side, resulting in an overall elliptical shape.

A similar behaviour was reported in the LES by Granet et al. [START_REF] Granet | Largeeddy simulation and experimental study of cycle-to-cycle variations of stable and unstable operating points in a spark ignition engine[END_REF] using the TFLES combustion model.

In order to assess the quality of the LES in terms of combustion, Fig. 3 shows the time evolution of the phase averaged ratio between the resolved flame surface S res = V |∇c|dV and the total flame surface S tot = V Σ c dV as defined in [START_REF] Wang | Large eddy simulation of a growing turbulent premixed flame kernel using a dynamic flame surface density model[END_REF], where c is the mass-weighted filtered progress variable such as ρ c = ρc.

A ratio S res /S tot close to 0 indicates a RANS-type simulation, where the flameturbulence interaction is entirely modeled, while a ratio close to 1 would correspond to a DNS. For the stab ref case, the spatial resolution of the mesh during the first 20 CAD is sufficiently high to resolve around 80 % of the flame surface.

In later stages, the spatial resolution becomes lower due to the mesh movement, yet still around 50 % of the flame surface is resolved. The same resolution is obtained in the two unstable cases as shown with the unst lean case.

pcyl (CAD) -2σ(CAD) < pcyl (CAD) < pcyl (CAD) + 2σ(CAD), ( 6
)
where pcyl is the phase averaged cylinder pressure at a given crank angle CAD

(normalised for all cycles between 0 and 720), and σ the standard deviation at a given crank angle resulting from a statistical phase averaged analysis of all individual engine cycles. For the sake of clarity, the symbol ˆwill be omitted thereafter.

Table 3 summarises the cycle averaged mean values and the variability of different variables resulting from a statistical analysis of the experiments and the present LES.

The variability is evaluated via a coefficient of variation (COV), defined as the standard deviation of the considered quantity X divided by the mean value, and expressed in %. It is computed for any quantity X as :

cov(X) = 1 N N i=1 X - N i=1 X/N 2 N i=1 X/N . ( 7
)
where N is the number of realisations.

The IMEP is the cylinder pressure integrated as a function of the cylinder volume and it is used as a measure of indicated work output per cycle of a cylinder:

IM EP = 1 V d p cyl dV, (8)
where V d is the cylinder displacement volume. Table 3 confirms that the LES accurately reproduces the experimentaly observed mean trapped mass, IMEP and maximum cylinder pressure.

Comparisons of the COV of these quantities confirm that the CCV levels are under-predicted by our LES, by 18 % for the IMEP and up to 50 % for the trapped mass. Considering that the trapped mass is not a direct experimental observation, but results from an a posteriori analysis of individual experimental cycles that introduces additional uncertainty, and that its variability is very small, the present results can be considered as overall validating the ability of the proposed LES methodology to predict combustion in the SGEmac engine under stable operating conditions. Note that a very similar level of reproduction of experimental findings for the stable operating point were reported in [START_REF] Granet | Largeeddy simulation and experimental study of cycle-to-cycle variations of stable and unstable operating points in a spark ignition engine[END_REF],

based on analysing 25 LES cycles.

LES of the unstable operating points

The developed LES methodology was applied to the study of combustion and its variability in the two unstable cases of the SGEmac database: unstab dil, resulting from the above reference operating point by an important dilution of the fresh gases by N 2 (mimicking an increase in the external EGR rate), and unstab lean, obtained from the reference point by a strong reduction of the overall fuel/air equivalence ratio.

These modifications of the fresh gases' composition have a strong effect on the laminar flame speed. Fig. 5 shows time evolutions of the statistical mean laminar flame speed averaged over the flame surface as a function of the crank angle.

This information was extracted from the LES discussed in detail hereafter following:

< S L > s = S L (x i , t)Σ(x i , t)dV Σ(x i , t)dV ( 9
)
where dV is the volume of the computational cell. As compared to the reference case, its value is strongly decreased for the unstable cases. Furthermore, its values for unstab dil are 10 to 30 % smaller than those encountered for the unstab lean case during the ignition and propagation phases. An opposite trend is apparent in the late phases of combustion. These strong reductions of the laminar burning characteristics (laminar flame speed, but also the laminar flame thickness not shown here) do strongly modify the response of the flame to the turbulent in-cylinder flow during the ignition and combustion phases, and are the key factor for increasing the CCV as compared to the reference case. Overall, the combustion proceeds in a qualitatively similar way to that observed for the stable reference case, with a generally faster propagation towards the exhaust side (negative x direction) as a result of the tumbling flow. However, an important level of variability is apparent for the two unstable cases, and important differences can be observed between slow and fast cycles. For the unstab dil case, this is especially apparent when comparing cycle 17, the slowest of the shown cycles, and cycle 22 which is the fastest. The flame shape for the former is very elongated towards the negative x direction and has hardly propagated in the y direction, while the shape for the faster cycle is more isotropic, despite a preferential propagation towards the exhaust side. At the shown instant, the flame of cycle 22 occupies around 25 % of the chamber volume, while this value only reaches around 12 % for the slowest cycle 17. Fig. 7 exhibits a similar difference between the slowest shown cycle 15 of unstab lean, occupying a small volume of the chamber and being very elongated towards the exhaust side, and the fastest cycle 16 that exhibits a more isotropic shape. A similar behaviour was reported in [START_REF] Granet | Largeeddy simulation and experimental study of cycle-to-cycle variations of stable and unstable operating points in a spark ignition engine[END_REF], and shown to result from interactions of the early flame kernel with the spark plug cavity and from characteristics of the mean flow field. A detailed analysis will be presented below.

Illustration of the combustion variability

Reproduction of the experimental CCV

In order to quantify the CCV of the two unstable cases, 30 consecutive full engine cycles were simulated for both cases, using the developed methodology.

The first two cycles of each case are discarded and the results are compared with experimental findings obtained by acquiring 200 consecutive cycles.

Figs. 8 and9 compare LES results with experimental findings in terms of cylinder pressure and of FMBe. The latter stands for fuel mass fraction burnt (energetic) and is obtained from a combustion analysis, performed using an inhouse analysis tool. This tool performs an instantaneous energy balance, during the combustion phase, based on the cylinder pressure trace, piston law, thermodynamic state of the trapped mixture, and on the wall heat losses estimated using a Woschni correlation [START_REF] Woschni | A universally aplicable equation for the instantaneous heat transfer coefficient in the internal combustion engine[END_REF]. The resulting overall instantaneous heat release rate is non-dimensionalized with the energy content of the trapped fuel mass, and integrated over time to yield the cumulative heat released that allows tracing overall combustion progress as well as its timing and duration. Both the experimental and simulated pressure traces were analyzed using this analysis tool. More quantitative comparisons are provided in the following Section.

Quantitative validation of the CCV predictions

The combustion process is usually divided in four main stages [START_REF] Ozdor | Cyclic variability in spark-ignition engines. a literature survey[END_REF]. First, the sparking and flame initiation correspond to the formation of a flame kernel during the spark life (the first 500 µs in the present situation), which can be successful or not, depending on the energy deposit and flow conditions in the vicinity of the spark. In our case, the flame kernel initiation is always successful.

Then the initial flame kernel development is defined as the period during which the burnt gas mass fraction reaches 1 to 2 %. It ends when the flame kernel reaches a size at which it is influenced by large eddies. These two stages correspond here to the ignition phase. The third phase is the turbulent propagation, followed by the last phase that corresponds to the end of combustion when the flame extinguishes at the combustion chamber wall.

In order to yield a more quantitative assessment of the quality of reproduction of the experimentally observed CCV by the present LES, Table 4 the COV of the maximum cylinder pressure is slightly under-estimated, while the COV of CA2 and CA50 are rather over-estimated. A major discrepancy appears for the COV of trapped mass. The reason for this deviation could not be identified, but one can note that in any case the experimental COV of trapped mass is small (0.3 %), so that it is not a major source of CCV.

For the diluted case, the mean and variability of the trapped mass are very accurately reproduced by the LES. The error of prediction of mean IMEP and maximum pressure are of the same order than for the lean case. This also holds for the mean CA2, the combustion progress being predicted slightly in advance, while the CA50 is quite well retrieved. In terms of variability, the COV levels of all parameters are slightly over-estimated by LES, yet satisfactory. The highest discrepancy in terms of COV is found for IMEP, where the LES over-estimates the experimental value by 20 %.

Analysis and quantification of CCV sources

Comparison of slow and fast cycles

In order to compare slow and fast cycles of the two unstable cases, the flame propagation is visualised for different CA degrees after the spark ignition using the same iso-surface as in Section 4 and vertical central cuts showing the reaction rate. The two cycles that are systematically compared are characterized by very different combustion speeds, but have similar trapped masses and gas composition. predicted by the LES of the diluted case. During the early ignition the flame is first driven into the spark plug cavity for both cycles, and both exhibit a similar initial flame kernel growth rate. Later, at 35 CAD after ignition, the flame of cycle 9 however propagates faster into the rest of the chamber than for cycle 17.

These differences coincide with the evolution of the flow velocity in the flame vicinity. The latter is characterized by spatially averaging the fresh gas velocity components along the flame front as:

< u(t)| c=0.001 >= u(x i , t)| c=0.001 dV δ( c -0.001)dV (10)
Fig. 11 right, shows the crank angle evolution of these axial (top) and vertical (bottom) fresh gases velocities. During early ignition, the axial velocity of cycle 17 is small compared to that for cycle 9, and there is a quite strong positive (i.e. upward) vertical component. This combination tends to counteract the propagation of the flame out of the spark cavity. In opposition, cycle 9 exhibits a clearly higher (and negative) axial fresh gases velocity, combined with a smaller and even negative vertical component, which both favour the propagation downward and away from the spark plug. The sudden increase of the axial velocity (in amplitude) that occurs 20 to 30 CAD after ignition for both cases coincides with the flame propagation outside the cavity in regions with higher velocities.

Fig. 12 proposes a comparable analysis for the lean case. In the slower cycle 21, the flame kernel can hardly propagate out of the spark plug cavity as clearly seen at 17 CA degrees after spark timing. Later, the flame is cut in two, but the flame propagating out of the cavity is small and a fraction is quenched at the walls. On the contrary, in cycle 18, the flame can propagate outside the spark plug cavity from the very beginning of the ignition, and grows much faster into the main part of the chamber. The comparisons of the fresh gas velocity components initially show similar axial velocities at ignition timing, but later on the vertical velocity for cycle 21 is positive and large thus counteracting the flame propagation away from the spark. On the contrary, the high negative vertical fresh gases velocity for cycle 18 greatly facilitates the propagation away from the spark plug. For both cases, the fresh gas velocity seen by the flame approaches zero when the flame front reaches the wall.

Fig. 13 presents plots of the CA2 versus CA50 for the LES and experiments of the stable and the two unstable cases. These curves show that on average, the smaller the CA2 (i.e. the faster the initial burning phase), the smaller the CA50 (the faster the main combustion phase). The standard deviation of the scatter plots is close to 1 degree in the stable case, 1.25 degree in the lean one and 2 degrees in the diluted case. The LES reproduce the overall experimental tendencies. It is interesting to note that some extreme cycles contradict the average correlation, both in LES and experiments, which will be studied below.

This first analysis of the CCV in the unstab dil and unstab lean cases reveal that the flame is strongly perturbed by the flow motion in particular during the early combustion phase, that in turn impacts the ensuing flame development and overal combustion duration.

CCV sources and correlations

The above presented analysis allows studying 3D phenomena as the local interactions of the flame with the flow and chamber geometry, but it does not help identifying overall dependencies. In order to achieve the latter in what follows, quantities describing the combustion efficiency and speed are correlated with local and global parameters in order to attempt identifying statistical correlations.

Thereafter all local quantities are spatial averages over a radius of 1 mm at the spark plug center at ignition timing, and are identified with the subscript spk .

To quantify the correlation, Bravais-Peason correlation coefficients are used:

R = σ xy σ x σ y (11)
where σ xy is the covariance of variables x and y, and σ x and σ y are their standard deviations. A value of |R| close to 1 implies that x and y are perfectly linearly correlated. A value of 0 indicates that there is no correlation. The sign of R indicates the sign of the slope of the regression line.

The spark energy as a source of cyclic combustion variability is discarded since a very high amount of energy is deposited after each breakdown to pre-vent misfires or ignition delays. Preliminary experimental and numerical tests performed in the two unstable cases confirmed that variations of the energy transferred to the gas and of the spark life duration has negligible impact on the combustion process.

In the diluted case, the fraction of fuel and air is held constant for all cycles. As the mixture is stoichiometric, an incomplete combustion leads to the same proportion of unburnt fuel and air, but to a different amount of unburnt fuel mass. The in-cylinder fuel mass of cycle n is therefore partially correlated between subsequent cycles as shown in Fig. 14. In the lean case, there is no correlation between subsequent cycles.

In order to quantify the impact of thermodynamic variability on the ignition process, correlation coefficients between different local and global variables and the local laminar flame speed at ignition timing were computed. The explored variables were : temperature, equivalence ratio, fraction of dilution gases Y dil , dilutation rate τ dil and residual gas rate τ res . The three latter quantities were computed as described hereafter. The combustion chamber contains a total trapped mass (m tot ) which is the sum of the residual gas mass (m res ) and of the mass of air and fuel induced during the intake stroke. The residual gases are composed of variable proportions of fresh and burnt gases left over from the previous cycle. The trapped mass m tot can also be defined as the sum of the fresh gases (fuel + air) and of dilutant (m dil ):

m dil = m CO + m CO2 + m H2 + m H2O + m N2 - 3.76W N2 Y O2 W O2 m tot , (12)
where Y O2 is the oxygen mass fraction, W N2 and W O2 the molecular masses of nitrogen and oxygen. The last RHS term of Eq. ( 12) is the mass of N 2 coming from the fresh air introduced in the combustion chamber during intake. The mass m N2 in Eq. ( 12) includes the mass of N 2 dilutant introduced during intake in the diluted case. The dilution and the residual masses are normalized by the trapped mass to yield respectively the dilution rate (τ dil = m dil /m tot ) and the residual rate (τ res = m res /m tot ).

The local fraction of dilutant Y dil is defined:

Y dil = Y CO + Y CO2 + Y H2 + Y H2O + Y N2 - 3.76W N2 Y O2 W O2 , (13)
The mean in-cylinder temperature is not deeply analysed since its variation is negligible (about 2K) since it is related to the quasi-constant values of the incylinder total mass and pressure at spark timing. As a result, the variations of the dilution rate and of the trapped fuel mass are perfectly correlated, as shown in Figs. [START_REF] Aleiferis | Flame chemiluminescence studies of cyclic variations and a/f ratio of the reacting mixture in a si engine[END_REF]. Tables 5 and6 summarize the different correlations between the aforementioned quantities for the two unstable cases. The highest correlations are written in red, the lowest values appear in yellow and intermediate ones in orange. The mean and standard deviations (in percent) of each variable is also provided in the first row. For the diluted point, S L is strongly correlated with the local and global dilution gas fraction, which are also highly correlated with each other. Indeed, the additional dilutant (N 2 ) is perfectly premixed with the fresh intake gases. During intake and compression strokes it mixes with the residual gases so that at spark ignition timing, the mixture is almost homogeneous. A moderate correlation is found between S L and the temperature with a negative value, the temperature being partially correlated with the dilution rate. This result seems in contradiction with the fact that the laminar flame speed depends on the square of fresh gases temperature (using Eq. ( 5)). Furthermore, the analysis of Pera et al. [START_REF] Pera | Effects of residual burnt gas heterogeneity on early flame propagation and on cyclic variability in spark-ignited engines[END_REF] shows that the effect of temperature heterogeneity is dominant over the dilution effect. An explanation can be found when examining Fig. 16 showing the evolution of the computed laminar flame speed as a function of the temperature and of the dilution rate within the range of values encountered at ignition timing. The maximum variation of the laminar flame speed reaches 26 % of the mean value over the range of variation of the dilution rate, while it reaches only 7 % with the temperature. Therefore, the impact of the dilution is stronger than the one of the temperature in these conditions. The local fluctuation of temperature is generated by both the wall heat transfer and the mixing with hot residual gases, while the variation of the local dilution rate composed of N 2 dilutant from the intake (initially cold) and hot burnt gases from the previous cycle is rather related to the variation of the global composition. Consequently, an increase of the dilution gas fraction is followed by a decrease of the flame speed. We may conclude that in this case, the combustion progress is rather influenced by the variation of global quantities (the amount of available fuel mass is perfectly correlated with the burnt gas mass) than that of the local temperature. For the lean case, Table 6 shows high correlations between the laminar flame speed and temperature, dilution gas fraction and the local and global equivalence ratio, which are also well correlated to each other. For this case, the burnt gas rate is weak (5 %), so that an increase of the temperature or of the burnt gas fraction leads respectively to an increase and a decrease of the laminar flame speed. The increase of the temperature is also well correlated with the increase of the residual gas rate composed of burnt gases and excess air.

Correlations with different local and global quantities describing the flow motion were then examined. Fig. 17 The effect of the cycle-to-cycle dependency of combustion on thermodynamic and aerodynamic parameters was then investigated by studying their correlations with IMEP. First of all, single linear regression models were used. Fig. 18 presents the resulting histogram of the correlation coefficients between IMEP and different variables X i for the two unstable cases. Similar patterns filling the bars indicate variables which are correlated to each other. The correlations with IMEP of the following variables were examined (in the order of appearance in Fig. 18):

• ( w/ u) spk , the ratio between the resolved vertical and axial velocity components

• φ inhom spk
, a measure of the resolved equivalence ratio inhomogeneity defined in Eq. ( 14)

• φ, the mean in-cylinder equivalence ratio

• m f u , the in-cylinder unburnt fuel mass after the end of intake

• Y dil,spk , the resolved mass fraction of dilutant defined in Eq. ( 13)

• m air , the in-cylinder air mass

• τ dil , the dilution rate defined using Eq. ( 12)

• | U | spk , the resolved norm of velocity

• E Ty , the tumble energy defined in Eq. ( 16)

• IM EP n+1 , IMEP of the subsequent cycle

• E c , the resolved kinetic energy

• m tot , the total in-cylinder mass

• u sgs spk , the subgrid scale velocity (at the combustion filter size)

• τ res , the residual rate

• T inhom spk
, a measure of the resolved temperature inhomogeneity defined in Eq. ( 15)

The local equivalence ratio and temperature inhomogeneities in the spark plug vicinity Φ spk,inhom and T spk,inhom were evaluated following a monoparametric regression study as:

Φ inhom spkp = Φ spk -a Φ, a = 1.0785 (14)
T inhom spk = T spk -a ′ T , (15)
where a ′ = 1.6177 for unstab lean and a ′ = 3.3645 for unstab dil were found.

The tumble energy E Ty was determined following:

E T y (CA) = m tot [T y ω eng D(ST )/4] 2 , (16)
D is the distance between the piston and the top of the cylinder head and ST is the crank angle of spark timing.

In Fig. ) together seem to have a second order but still significant impact. It is noteworthy that | U | spkplg is well correlated with the axial velocity u which is generated by the tumble motion.The correlation of the IMEP with the subsequent cycle IM EP n+1 is weak. Other quantities are negligible : the total enclosed mass which varies only weakly, the residual gas rate and the local flow direction, described by the ratio between the vertical velocity w and the axial velocity u.

On the contrary, for the lean case, Fig. 18(b) shows that the main source of the IMEP variability is clearly the ratio w/ u with a coefficient of correlation |R| = 0.86. Then thermodynamic variations influence partially the combustion variability. For both cases similar correlations are found between the IMEP and the enclosed fuel mass (0.54) and the dilution gas rate (0.49).

Multidimensional data analysis

The previous analysis has highlighted the relationships between combustion characteristics and flow/thermodynamic parameters. However the single regression analysis reaches a limit when several parameters act simultaneously, as it is the case in the present results, making it difficult to identify the most important factors and to obtain relevant correlation coefficients. In order to achieve a more relevant overall picture of the most important parameter impacting CCV, a multivariate regression model was used. In the proposed approach, a relationship between the dependent response IM EP (n), where n is the cycle number, and different independent parameters or variables X i (n) is explicitely estimated.

As the type of regression (linear, gaussian, polynomial...) between the variables is unknown, and as in the present cases the standard deviation of any variable did never exceed 8 %, we assumed here a linear multiple regression.

To this purpose,the following steps were followed:

1. First the necessary condition for selecting the X i among all the possible parameters is that the X i vectors are non collinear, and that n > p + 1, where p is the number of variables. These variables may be defined based on a correlation analysis between every all available variables X i in order to determine whether quantites are collinear or not. This step could also have been carried out using Principal Component Analysis [START_REF] Rencher | Methods of Multivariate Analysis[END_REF]: it is a mathematical tool to reduce data space dimensions and compute linear combinations of the initial variables that are called principal components and which are uncorrelated to each other. However this method is beyond the scope of the present study as we need to manipulate physical quantities, and not complex linear combinations between them.

The matrix equation system IM EP

(n) = p i=1 β i X i (n) + ǫ(n)
is solved using a least-squares fit which provides the constant slopes β i and the residuals ǫ(n) of the equation. The X i (n) are normalised by their mean value to avoid that any variable has a larger impact due to its scale.

3. The degree of correlation between X i and the response IM EP is examined by calculating the coefficient of correlation between X i and IM EPp k =i β k X k , thus discarding the contribution of the other variables.

The degree of correlation between IM EP and the resulting model

i β i X i (n)+ ǫ(n) is examined. A correlation coefficient
R greater than 0.7 is usually considered as strong correlation.

5. If the contribution of variable X i is not significant in the resulting model (typically |R XY | < 0.4), then it is suppressed from the equation and step

(2) is repeated. If several sets of variables can be identified, steps (1) to [START_REF] Chen | Impact of ignition energy phasing and spark gap on combustion in a homogenous direct injection gasoline si engine near the egr limit[END_REF] are performed with each of them to find the most relevant one, providing the final, best fitting model.

Following steps (1) to (5), the following set of independent variables was identified :

• For the unstab dil case: m f u , | U | spk , u sgs spk , E Ty and T inhom spk .

• For the unstab lean case:

m f u , | U | spk , ( w/ u) spk , T inhom spk , φ inhom spk .
The modeled cylinder pressures are finally obtained as : For most variables, the correlation obtained from the analysis confirms the findings from single regression, but the level of correlation is generally higher. The noticeable exception is the correlation for temperature inhomogeneities at the spark plug : the multivariate correlation is much higher for unstab lean than for the single regression, the opposite trend appearing for unstab dil. Overall this parameter is found negligible for both cases in the multivariate analysis. In summary and to first order, the multivariate analysis indicates that for the diluted case, the main parameters correlated with IMEP variablity are in-cylinder fuel mass, local resolved and subgrid-scale velocity at spark and the tumble energy. In the lean case, IMEP variability is mainly correlated to variations of local resolved velocity components and amplitude, and to the in-cylinder fuel mass. Temperature inhomogeneities at the spark plug do have a weak impact.

IM EP model (n) = i=1,5 β i X i (n) + ǫ i with X i = (m f u , | U | spk ,

Qualitative classification of CCV

Finally, we compared our qualitative findings concerning the relative importance of different sources of CCV for the lean and diluted cases of the studied engine with the ones identified by Ozdor et al. [START_REF] Ozdor | Cyclic variability in spark-ignition engines. a literature survey[END_REF]. They reviewed different indicators of CCV from a substantial number of experimental studies. They involved two and four stroke, single and four cylinder SI engines, running on various fuels (iso-octane, propane, gasoline, toluene, methanol) either in premixed, carbureted or port fuel injection mode. Furthermore, different sparkignition systems and spark plug locations, different dilution rates by nitrogen (0 to 15%), several equivalence ratios (lean and stoichiometric), different combustion chamber geometries, piston speeds, loads and swirl ratios were included.

The result was a qualitative estimation of the phenomena influencing the CCV levels.

Following the results achieved in the present work, a similar classification could be obtained for the studied engine, but in our case based. The obtained results are summarized in Table 7, and compared with findings from Ozdor et al. [START_REF] Ozdor | Cyclic variability in spark-ignition engines. a literature survey[END_REF], distinguinshing the impact on the spark-ignition and main combustion phase respectively.

Concerning spark-ignition, the convection in the spark plug region is the variable having the highest impact on CCV, which confirms findings by Ozdor et al. For the diluted case, its importance is comparatively smaller and of comparable importance to that of the global dilution rate. No other main factor was identified by our LES, in constrast to the lean case and to work by Ozdor et al. For the latter, comparable secondary levels of importance are found for the the global equivalence ratio and dilution rates, and to a lesser extend for spatial mixture inhomogeneities at the spark plug. Generally speaking, the classification proposed by Ozdor et al. is quite similar to the one obtained for the lean case, but exhibit some discrepancies with the diluted case.

The effect of composition variation affecting the laminar flame speed and therefore the beginning of the combustion is second order in the lean case through the equivalence ratio, the temperature and the burnt gas fraction and first order in the diluted case through the burnt gas fraction. During the whole aries of the LES domain from system simulations of the complete engine set-up.

The proposed LES methodology based on the AVBP code was first shown to reproduce accurately the experimental findings for the stable reference operating point. It was then applied to perform multi-cycle LES comprising 30 full consecutive engine cycles for each of the two unstable cases. Detailed comparisons with experiments allowed validating the ability of the proposed methodology to accurately reproduce the phase-averaged mean and variability of in-cylinder pressure and combustion progress, as well as of other combustion characteristics and in particular IMEP.

The LES results were then post-processed with the aim of identifying causes for the observed CCV. A qualitative comparison between slow and fast burning cycles of the two operating points allowed illustrating how a combination of local velocity fields and geometrical features of the chamber as the spark plug cavity in the pentroof head concur to generate cyclic variability.

In order to yield a more quantitative estimation of the relative importance of different parameters on the occurrence of CCV, their individual impact was evaluated using single parameter regression models for their dependency on the cyclic varying IMEP. It was shown that despite the overall similarity between the diluted and lean operating point, the ensemble of most relevant parameters differ notably for the two cases. As the single regressions can hardly be used in a systematic approach to identify major causes of CCV, a multivariate regression was proposed and applied to both cases. It allowed identifying in a more systematic way the main parameters impacting IMEP variability. It also provided an explicit model linking the latter to first order to the variations of a limited set of independent, both local and global, parameters.

Finally a qualitative comparison of the relative importance of different local and global parameters as sources of CCV for the two studied cases was compared with findings from an extensive experimental survey of CCV sources by Ozdor et al. [START_REF] Ozdor | Cyclic variability in spark-ignition engines. a literature survey[END_REF]. Overall the latter findings could be confirmed by the presented analysis of LES, but also indicated that the detailed causes of CCV will probably depend on the type of engine and combustion mode. This approach may bring objective elements to engineers for better understanding CCV sources on a given engine design or specific operating conditions. Such information is essential for improving design rules and is inaccessible to experiments. In this sense the proposed exploitation of LES is today unique and opens a large application field to this CFD tool in the piston engine industry. Future work will concern comparative studies for gasoline direct injection engines, in which the complex interactions between the fuel spray and intake aerodynamics will require specific LES studies and could potentially require more sophisticated regression models (quadratic or non linear). CA deg after tign

Section 2

 2 starts by describing the experimental set-up and main characteristics of the SGEmac engine. Section 3 then details the numerical set-up of the LES of the SGEmac engine, and gives key elements of the proposed LES methodology. A first validation of the latter is provided in Section 4, which compares LES predictions of phase-averaged cylinder pressure and its variability for the reference low-CCV (or stable) operating point with experimental findings. Section 5 then presents the reproduction of experimental CCV findings for two high-CCV (or unstable) operating points of the SGEmac database achieved with the developed LES methodology. A first investigation of the CCV origins is proposed based on the visualisation of the flame propagation and by examining the degree of correlation between IMEP and the different combustion and flow characteristics. Finally, a multivariate regression model built from LES results, allows to quantify the relative importance of different local and global CCV sources and classify the degree of importance of several phenomena acting on CCV.

Fig. 1 (

 1 Fig. 1(b) illustrates the achieved discretisation of the computational domain.

Fig. 4

 4 Fig. 4 compares the crank angle evaluations of the cylinder pressure predicted by the 10 LES cycles with experimental findings resulting from 100 engine cycles. Fig. 4(a) first superimposes the individual cycles predicted by LES on the experimentally acquired data. Despite the limited number of simulated

Figs. 6

 6 Figs.6 and 7illustrate the cyclic variability predicted by the present LES for the unstab dil and unstab lean cases respectively. Top views of the combustion chamber are shown for 6 consecutive cycles taken at 35 CAD after spark-ignition, the flame being materialised by an iso-surface of the reaction rate coloured by the flow velocity. The time evolutions of the cylinder pressure of these cycles are also shown to distinguish slow burning cycles, leading to low cylinder pressure levels, from fast burning cycles exhibiting higher cylinder pressures.

Fig. 8 (

 8 Fig. 8 (a) superimposes cylinder pressure traces for individual LES cycles for the unstab dil case and an experimental probability density function (PDF) resulting from the observed probability of reaching a specific cylinder pressure at a given crank angle. Fig. 8 (c) shows the corresponding FMBe curves and experimental PDF. Qualitatively, the cycles predicted by LES seem to reproduce the high cyclic variability observed in the experiments. However it is also apparent that the LES has a tendency to predict overall slower cycles than in the experiments, no LES cycle corresponding to the fastest experimental cycles, and one LES cycle being slower than statistically observed. A quantitative analysis is proposed in Fig. 8 which compares the mean and statistical envelope of cylinder pressure predicted by LES with experiments. The evolution of the average

Fig. 11

 11 Fig. 11 left, compares results for the fast (C9) and the slow cycle (C17)

Fig. 17

 17 Fig. 17(c) to (f) correlates the resolved velocity magnitude and the axial

  Fig.17 (c) to (f) correlates the resolved velocity magnitude and the axial resolved velocity u in the vicinity of the spark at ignition timing with the tumble number. Surprisingly, no correlation is found between the local velocity amplitude and the tumble number. u and T y are only partially correlated in the lean case (R = 0.43), while no correlation is found in the diluted case. However we notice that the slowest cycle (C15) of the unstab lean case is characterized by both very low local velocity at spark and tumble number. As seen in Sec-

  18, lines corresponding to a coefficient of correlation |R| = 0.4 are plotted to indicate approximately the level above which correlations are significant. For the diluted case, Fig. 18(a) indicates that the thermodynamic quantities have the most important effect on the IMEP in terms of cyclic variations and that they are all correlated with each other. Aerodynamic quantitites characterizing the large scale turbulence (tumble energy E Ty and kinetic energy E c ) and the local resolved and subgrid scale velocity (| U | spkplg and u sgs spkplg

•

  u sgs spk , E Ty , T inhom spk ) in the diluted case and X i = (( w/ u) spk , | U | spk , m f u , φ inhom spk , T inhom spk ) in the lean case. The obtained best fitting models for the two unstable cases are shown in Figs. 19 and 20 respectively for the diluted and lean cases. The first bar on the left in addition gives the correlation coefficient between IMEP and the obtained best fitting model IMEP model , thus measuring the achieved degree of approximation by the multivariate model. The correlation coefficient between IMEP and its models is 0.93 for unstab lean, indicating a high accuracy. For the unstab dil case it is 0.8, showing a still sufficient yet slightly less accurate reproduction of the observations. A reason for this discrepancy is probably that the flame is more strongly perturbed by the turbulent flow during the whole combustion phase, so that the diluted case could require taking into account other relevant physical parameters to improve the model accuracy. This was not attempted in the present study. Two correlation coefficients between IMEP and the identified variables X k are shown in Figs. 19 and 20: The left bars filled with dots indicate correlation coefficients obtained between the X k and IMEP via the single regression model; • The right bars filled with lines indicate the correlation coefficients from the multivariate analysis between the independent variables X k and IM EPi =k β i X i (n), i.e. discarding the contributions of the other idependent variables.

Figure 1 :

 1 Figure 1: (a) Computational domain for the LES (rose frame), including the combustion chamber and a part of the intake and exhaust ducts; (b) View of the tetrahedral mesh during the intake stroke (-290 CAD).

Figures - 40

 40 Figures

Figure 4 :Figure 5 :

 45 Figure 4: Comparisons of the time evolutions of the cylinder pressure predicted by LES with experimental findings for the stab ref case.

Figure 6 :

 6 Figure 6: Visualisation of the flame surface using an iso-surface of the reaction rate (50 mol.m -3 .s -1 ) coloured by the velocity (blue: 0m/s, red: 20m/s) 35 CA degrees after ignition for 6 consecutive cycles of unstab dil.

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Visualisation of the flame surface using an iso-surface of the reaction rate (50 mol.m -3 .s -1 ) coloured by the velocity (blue: 0m/s, red: 20m/s) 35 CA degrees after ignition for 6 consecutive cycles of unstab lean.

Figure 10 :

 10 Figure 10: Matekunas diagrams for the three studied cases.

Figure 11 :

 11 Figure 11: Comparison of the flame propagation of a fast (C9) and a slow (C17) cycle of the unstab dil case. Left: vertical central cut showing the instantaneous reaction rate fields at four instants during the combustion phase. Right : crank angle evolutions of the fresh gas velocity components along the flame front (axial < u| c=0.001 > and vertical < w| c=0.001 >)following Eq. (10).

Figure 12 :

 12 Figure 12: Comparison of the flame propagation of a fast (C18) and a slow (C21) cycle of the unstab lean case. Left: vertical central cut showing the instantaneous reaction rate fields at four instants during the combustion phase. Right : crank angle evolutions of the fresh gas velocity components along the flame front (axial < u| c=0.001 > and vertical < w| c=0.001 >)following Eq. (10).

Figure 13 :

 13 Figure 13: Comparisons between experimental and LES correlations between CA2 and CA50.

Figure 14 :

 14 Figure 14: Scatter plot of the enclosed mass of fuel (bottom) of the subsequent cycle (n + 1) as a function of cycle (n).

  and of the dilution rate for the conditions encountered at ignition timing in the diluted case.

Figure 17 :

 17 Figure 17: Scatter plot of the tumble number as a function of the average in-cylinder resolved kinetic energy (top), of the resolved velocity magnitude (middle) and of the axial velocity (bottom) in the vicinity of the spark plug at ignition timing.

  details the cycle averaged mean and variability of different parameters characterizing the two unstable cases, The CA2 and CA50 values are defined as the crank angle at which the FMBe reaches respectively 2 % and 50 %. The CA2 is representative of the ignition phase while CA50 corresponds to the fully turbulent

	are reproduced within a 6 % error margin. The differences between LES and
	experiments in this case are related to a slightly too high overall combustion
	speed predicted by the simulations. LES predicts a mean CA2 which is 1.2 CAD
	in advance of the experiments, and up to 2.6 CAD in advance for the mean
	CA50. The level of variability predicted by the LES is generally satisfactory:
	propagation.
	Overall, the more accurate reproduction of experimental findings is obtained
	for the lean case: the mean trapped mass, IMEP and maximum cylinder pressure

Table 1 :

 1 Main characteristics of the SGEmac single cylinder engine. Crank Angle Degrees (CAD) are relative to combustion Top Dead Centre (TDC). Valve timings are given at a reference valve lift of 0.7 mm.

	Displacement		441 cm 3
	Bore		82 mm
	Stroke		83.5 mm
	Connecting rod length	144 mm
	Inlet Valve Open (IVO)	-348 CAD
	Inlet Valve Close (IVC)	-142 CAD
	Exhaust Valve Open (EVO)	140 CAD
	Exhaust Valve Close (EVC)	348 CAD
		stab ref unstab dil unstab lean
	F/A equivalence ratio	1.00	1.00	0.58
	Dilution by N 2 (% mass)	0	33	0
	IMEP (bars)	3.0	3.1	3.3
	Spark timing (CAD)	-20	-50	-52
	Engine speed (rpm)	1200	1200	1200

Table 2 :

 2 Characteristics of the three simulated fired operating points stab ref, unstab dil and unstab lean of the SGEmac database.

	Tables

Table 3 :

 3 Cycle averaged mean and variability of engine parameters for the stable stab ref case.

		unstab dil	unstab lean
		Experiment LES Experiment LES
	Number of cycles	200	28	200	28
	Mean trapped mass (mg)	259	259.2	277	280.4
	COV(trapped mass) (%)	0.1	0.1	0.3	0.1
	Mean IMEP (bar)	3.1	2.9	3.3	3.1
	COV(IMEP)(%)	7.6	9.1	8.2	8.1
	Mean p max (bar)	16.4	16.2	18.9	19.2
	COV(p max )(%)	12.5	12.9	14.5	13.8
	Mean CA2 (deg.)	-7.3	-8.0	-7.9	-9.1
	COV(CA2) (%)	5.2	6.3	5.9	6.4
	Mean CA50 (deg.)	21.3	20.5	19.0	16.6
	COV(CA50) (%)	7.6	8.5	7.7	8.3

Table 4 :

 4 Cycle averaged mean and variability of engine parameters for the unstable unstab dil and unstab lean cases.

	y / x	SL	T spk	Y dil,spk	τ dil	τres	Φ spk	Φ
	mean(x)	0.23m/s	549 K	0.39	39.1%	7.3%	1	1
	2σ (in %)	(8.6%)	(1.1%)	(2%)	(1.2%)	(1.5%)	0%	(0%)
	SL	1	-0.42	-0.93	-0.84	0.41	N.A.	N.A.
	T spk	-0.42	1	0.59	0.52	-0.51	N.A.	N.A.
	Y dil,spk	-0.93	0.59	1	-0.9	-0.48	N.A.	N.A.
	τ dil	-0.84	0.52	-0.9	1	-0.28	N.A.	N.A.
	τres	0.41	-0.51	-0.48	-0.28	1	N.A.	N.A.
	Φ spk	N.A.	N.A.	N.A.	N.A.	N.A.	N.A.	N.A.
	Φ	N.A.	N.A.	N.A.	N.A.	N.A.	N.A.	N.A.

Table 5 :

 5 Correlations coefficients R between several local and global thermodynamic variables at ignition timing in the diluted case.

	y / x	SL	T spk	Y dil,spk	τ dil	τres	Φ spk	Φ
	mean(x)	0.26m/s	515.5 K	0.05	5.1%	7.6%	0.56	0.56
	2σ (in %)	(3.9%)	(1.24%)	(13.69%)	(7%)	(3.2%)	(1%)	(0.7%)
	SL	1	0.77	-0.72	-0.61	0.32	0.71	0.6
	T spk	0.77	1	-0.4	-0.29	0.59	-0.1	0.033
	Y dil,spk	-0.72	-0.4	1	0.58	-0.29	-0.54	-0.39
	τ dil	-0.61	-0.29	0.58	1	-0.29	-0.63	-0.87
	τres	0.32	0.59	-0.29	-0.29	1	-0.16	-0.001
	Φ spk	0.71	-0.1	-0.54	-0.63	-0.16	1	0.76
	Φ	0.6	0.033	-0.39	-0.87	-0.001	0.76	1

Table 6 :

 6 Correlations coefficients R between several local and global thermodynamic variables at ignition timing in the lean case.

	Sources = variation of :	Effect on spark-ignition
			unstab lean	unstab dil	Ozdor et al. [1]
	Global A/F ratio	++	N.A.	++
		φ
	Global dilution	++	++	++
		τ dil
	Spatial inhomogeneities	+	-	+ to +++
	φ inhom spk	and T inhom spk
		Flow	+++	++	+++
	(w/u) spk , |U | spk
	Sources = variation of :	Effect on the main propagation phase
			unstab lean	unstab dil	Ozdor et al. [1]
	In-cylinder charge	++	++	+++
		m f u
	Global A/F ratio	++	N.A.	++
		φ
	Global dilution	++	++	++
		τ dil
		Flow	-	+	++
		E Ty

Table 7 :

 7 Qualitative impact of different sources of CCV for the studied lean and diluted case:

first those impacting the spark-ignition, second those impacting the main propagation phase are listed. The meaning of the symbols is as follows: +++: very strong correlation with (|R| > 0.7 for our LES findings), ++: moderate correlation ( |R| > 0.5 for our LES findings), +: weak correlation ( |R| ≈ 0.4 for our LES findings), -: no measurable impact and N.A. : not applicable.

http://pantar.cerfacs.fr/4-26334-The-AVBP-code.php

Acknowledgements

This work was financially supported by the European Union under project N. 233615 LESSCCV of the 7 th framework programme for research (FP7).

The project was granted access to the HPC resources of CINES under allocation N. 026139 from the GENCI (Grand Equipement National de Calcul Intensif) eDARI program. We also thanks Abdelaziz Faraj and Frédéric Nicolas from IFPEN for valuable discussions on the multivariate analysis.

propagation, the variation of the intake energy (fuel mass) has a strong impact in both cases. The tumble motion has moderate effect in the diluted case while it is negligible in the lean case. Concerning the impact on the main combustion phase, the lean and diluted cases do not exhibit any single major factor having a higher relative impact than others on CCV. The present LES underestimates the relative importance of factors as in particular the overall fuel mass, mixture inhomogeneities and flow energy as compared to Ozdor et al.

Globally, the LES results confirm major findings by Ozdor et al., despite some important differences do appear. These differences should however not be over-interpreted, as Ozdor's findings were the result of a compilation of a wide range of engines and modes of operation, while the present LES concerned only two operating points of a single engine type.

Conclusion

This paper presents comparisons between Large-Eddy Simulations and experimental findings on cyclic combustion variability in a single cylinder sparkignition engine fueled with a premixture of air and gaseous propane. Three operating points from the SGEmac database were studied : a stable reference operating point with low CCV levels, and two unstable points with high CCV, obtained from the reference case either by reducing the fuel/air equivalence ratio (lean case), or by adding an important external dilution by nitrogen (diluted case).

In the presented LES, spark-ignition was modeled using ISSIM-LES [START_REF] Colin | A spark ignition model for large eddy simulation based on a fsd transport equation (issim-les)[END_REF], which describes the electrical circuit of the spark and provides 3D modeling of the flame kernel growth during ignition, coupled to the ECFM-LES [START_REF] Richard | Towards large eddy simulation of combustion in spark ignition engines[END_REF] turbulent combustion model which describes the flame front propagation.

A simulation methodology was employed that imposes in/outflow bound-SI+5CAD SI+10CAD SI+15CAD SI+20CAD indicate the correlation coefficients from the multivariate analysis between the independent variables X k and IMEPi =k β i X i (n).


V N Zaytsev

Total P Biver

H Wackernagel

We B01 Geostatistical Simulations on Irregular Reservoir Models Using Methods of Nonlinear Geostatistics

Classical approaches to geostatistical simulations are not applicable directly on irregular reservoir models (such as Voronoi polygon and tetrahedron meshed models). One of the main difficulties is that the block marginal distributions are unique for every block due to volume support effect. We propose a methodology for geostatistical simulations which overcomes this difficulty in an analytical manner and provides a robust utilization of the small support petrophysical property distribution and the covariance model for irregular reservoir models. The proposed solution is based on the discrete Gaussian (DGM) model and operates directly on blocks of the target grid. This solution is also capable to improve the quality of the classical reservoir models, such as tartan meshes, by including the volume support effect into consideration and thus-providing geologically more realistic results. Applications to Voronoi polygon grid with local grid refinements and to a tartan-meshed offshore gas reservoir model are demonstrated.

Introduction

Geostatistical simulation on irregular grids is a challenging task, since the algorithms face a multiscale problem posed by the transition from the small scale core analysis and seismic data to the multiple uneven supports of the grid cells. Classical sequential Gaussian simulation (SGS) is not applicable on irregular support, since the data does not average linearly after the normal score transform. One approach to overcome this difficulty is using direct sequential simulation (DSS) algorithms, which avoid the transformation to the space of Gaussian variables and work directly in the space of the simulated variable [START_REF] Manchuk | Direct Geostatistical Simulation on Unstructured Grids[END_REF][START_REF] Soares | Direct sequential simulation and cosimulation[END_REF]). Another approach for simulating on irregular support with the methods of non-linear geostatistics was proposed by Brown et al. (2008). The authors demonstrate the application of discrete Gaussian model (DGM) for simulating the Cox process using multiple sample supports representing the samples as a tree graphical model and using an extended conditional independence assumption. The solution is based on a theoretical model proposed by [START_REF] Emery | On some consistency conditions for geostatistical change-of-support models[END_REF] for conditional simulations on regular support using DGM.

In this paper we provide an extension of the DGM-based geostatistical simulation algorithm for petroleum reservoirs, which uses the small support covariance derived from the data analysis and approximates with high accuracy the block-support marginal distributions. We first present two methods of transforming the problem of simulation on an irregular grid to the problem of simulating multivariate Gaussian random vectors. Then we demonstrate an application for simulating porosity on a tartan-meshed offshore gas field with significant volumetric differences between the cells.

Discrete Gaussian model for geostatistical simulations

Consider the problem of simulating on an irregular grid , 1 … the volumetric averages

| |
of the point-support random variable (RV) with a covariance function , ′ . Without additional assumptions, the marginal distribution of cannot be determined analytically from the distribution of and a Monte Carlo approach is applied, which in this context is referred to as mini-model simulation. In case of irregular grids, the mini-model approach is hardly applicable, since it requires multiple fine scale simulations for every block of the model. The proposed simulation model enables to derive the blocks marginal distributions in an analytical form by introducing additional theoretical assumptions. The model assumes that is a transformation of a standard Gaussian RV with a covariance function , ′ and that every block value is a block-dependent transformation of a standard Gaussian RV and , 1 …

The transform function can be fitted to the core sample values and is considered to be known. The DGM approach enables to shift from the problem of simulating the set of

∑ ⇒ ∑ (2)
Where the change of support coefficient for every block is defined by the equation

1 | | , (3)
And the covariance matrix for , … is defined by equations

1 | || | , , (4)
The multivariate Gaussian random vector , … can then be simulated using standard techniques, for instance SGS [START_REF] Chilès | Geostatistics: Modeling Spatial Uncertainty[END_REF]. Applying the transformation functions , 1 … gives an unconditional simulation of block average values , 1 … honoring the covariance function , ′ and approximating with high accuracy the block marginal distributions [START_REF] Chilès | Validity range of the discrete Gaussian change-of-support model and its variant[END_REF][START_REF] Chilès | Geostatistics: Modeling Spatial Uncertainty[END_REF][START_REF] Matheron | Change of Support for Diffusion-Type Random Functions[END_REF]. The simulated unconditional field can then be conditioned using block or point-support samples through simple or ordinary kriging and uniform conditioning (de Fouquet 1994). Another convenient way of applying DGM for geostatistical simulations consists in introducing an additional assumption, which leads to a weak DGM model -DGM 2

iii.

for any two blocks and , the joint distribution of

1 | | , (5)
, 1 1 , ′ (6)
DGM 2 enables establishing the simulation model in a simpler manner relative to DGM 1 by the cost of introduction a small bias into the model. Formally, using DGM 2 for simulation does not guarantee the correct reproduction of the covariance of the random vector , … due to assumption iii which is discussed by [START_REF] Chilès | Validity range of the discrete Gaussian change-of-support model and its variant[END_REF].

Difference between DGM 1 and DGM 2

In order to demonstrate the difference between DGM 1 and DGM 2, we provide statistical analysis of 100000 unconditional simulations on a Voronoi polygon 2D grid with 10 local grid refinements. The grid size is 20x20 km 2 with a vast diversity of block sizes demonstrated on Figure 1. In order to check the accuracy of DGM 1 and DGM 2, we compare the block marginal distributions implied by these models with the results of 50 000 Monte Carlo simulations on 2x2m 2 support. The difference between the observed block variance and the theoretical variance is also illustrated (Figure 2). It reveals that DGM 2 introduces a non-linear bias into the block distribution variance, which, however, is below 5% of the point support variance of .

Simulation of porosity on an offshore field tartan grid

Addressing the change of support problem in geostatistical algorithms not only provides a solution for populating fully unstructured grids (such as Voronoi and tetrahedron meshing) but also enables to ameliorate the simulation results for classical tartan grids. Consider a tartan mesh of an offshore gas field X with dimensions 140x200x2 km 3 (Figure 3). The classical approach consists of defining an average support for the model and the marginal distribution of porosity on this average support using Monte-Carlo simulations. The classical simulation proceeds in the IJK space, which leads to incorrect reproduction of block distributions and covariance. In order to treat the volume support effect in the simulation process, a proper simulation space is required in which the hypothesis of the stationary covariance model is applicable and the meshing structure of the reservoir is preserved. We adopt the GeoChron model [START_REF] Mallet | Elements of Mathematical Sedimentary Geology: the GeoChron Model[END_REF] corresponding to the target reservoir meshed model to perform DGM-based simulations (Figure 3).

Petroleum Geostatistics Biarritz, France, 7-11 September 2015

Conclusions

We propose a framework for geostatistical simulations on irregular reservoir grids using DGM. The proposed framework enables exact integration of the point-support covariance model on irregular support and approximating with high accuracy the block marginal distributions. It also enables to avoid the usage of the "average support" concept for simulations and to work on the target mesh directly. In addition, the proposed approach eliminates the artefacts imposed by the mesh, providing geologically more realistic full-size models.

Figure 1

 1 Figure 1 Voronoi-meshed 2D grid with 10 local grid refinements. Left) The block volumes and the well locations on the studied grid. Right) A simulation produced by DGM 2. The grid contains 3546 blocks; the smallest and the largest block sizes are 36x42 and 1035x1052 m respectively. The simulated variable has a standard lognormal distribution with logarithmic

Figure 2

 2 Figure 2 Quality check for DGM 1 and DGM 2. Left) The distribution for the largest block in the model. Right) The observed mismatch of the block variance versus block volume.

Figure 3

 3 Figure 3 Tartan meshed offshore gas field model. Left) Original model. Right) Corresponding GeoChron model.

Figure 4

 4 Figure 4 Result comparison. Left) Classical simulation in IJK on "average support". Right) Simulation in GeoChron model with DGM-based algorithm.