3d Modeling of Mixing, Ignition and Combustion Phenomena in Highly Stratified Gasoline Engines

Résumé — Modélisation 3D du mélange, de l'allumage et de la combustion dans les moteurs à essence fortement stratifiés — Le présent article décrit les développements récents réalisés à l’ IFP sur la modélisation 3D de la combustion dans les moteurs à allumage commandé (SI, spark ignition ). Ces développements ont permis d’obtenir le modèle ECFM (extended coherent flame model) qui constitue une extension du modèle de flammes cohérentes CFM (coherent flame model) , particulièrement adaptée à la simulation de la combustion dans les moteurs à allumage commandé et à injection directe (DI-SI). L’idée principale de cette extension consiste à décrire localement la richesse, la composition (incluant les gaz résiduels) et la température des gaz frais. Une généralisation de ce modèle à des carburants multicomposants est ensuite proposée. L’effet des fluctuations à petites échelles est intégré dans le modèle ECFM à travers une combinaison d’un modèle de variance/dissipation scalaire et d’une méthode de densité de probabilité présumée. Le modèle d’allumage AKTIM, développé pour modéliser l’initiation de la combustion à la bougie, a été amélioré afin de prendre en compte le plissement du noyau de flamme


INTRODUCTION
The direct injection (DI) concept is a promising way to reduce fuel consumption and pollutant emissions of spark ignition (SI) gasoline engines.But in order to bring this new concept to a final industrialisation, some problems not relevant for the conventional port fuel injection (PFI) concepts have to be addressed by the engine designers.First, globally very lean mixtures are used at part load conditions.They are detrimental to the efficiency of the after-treatment catalyst system, usually laid out for stoichiometric mixtures.It is therefore important to reduce by design as much as possible the formation of pollutants in the combustion chamber.The second issue is the high sensitivity of consumption and pollutant formation on the chamber design and injection parameters.These two issues make the design of an efficient DI-SI engine a serious challenge for engineers.
In this context 3D computational fluid dynamics (CFD) codes are a very useful design tools, able to help reduce dramatically the development cost of such new engine concepts.It allows to explore and better understand the basic processes occurring inside the engine, and how they relate to global engine parameters.3D CFD also allows for a fast qualitative comparison of different promising designs.
Unlike PFI-SI engines mostly running under homogeneous conditions, operating conditions in DI-SI engines are characterised by a high charge stratification.This poses a new challenge for 3D CFD modeling, as the accurate description of this high stratification and its effects on spark ignition, combustion and pollutant formation not only needs improvement of the description of average flow quantities, but also the inclusion of the effects of fluctuations around these average quantities.This paper presents an overview of the developments made in this sense at IFP.We only present the last achievements on improving the description of combustion under DI-SI engine conditions, leaving out the work on improving liquid spray modeling.All these model developments were made in the context of the IFP 3D CFD engine simulation code KMB [1, 2] a multi-block structured in-house evolution of the KIVA-II code [3].
In the first part of this paper, recent developments of the extended coherent flame model (ECFM) are presented: -accounting for large scale stratification through an improved local evaluation of the fresh and burned gases state; -adaptations to include multi-component fuel effects on combustion; -accounting for small scale stratification through the formulation of a variance/scalar dissipation model and the inclusion of their effects on combustion; -improvements to the predictivity of the AKTIM spark ignition model by accounting for turbulence effects in early phases of combustion; -and finally, accounting for knock in a multi-component fuel context.
In the second part of the paper, these models are applied to the simulation of two engine configurations with the aim of validating their predictions, comparing them when possible to experimental results.A DI-SI optical access engine, based on a Renault PFI engine, is simulated in order to validate the multi-component, knock and scalar fluctuation models.A Mitsubishi GDI engine is simulated to assess the correct behaviour of the improved AKTIM model under high turbulence conditions, and to evaluate the effect of equivalence ratio fluctuations on combustion.

THE EXTENDED COHERENT FLAME MODEL ECFM
The coherent flame model (CFM) is a combustion model adapted to the flamelet regime.This modeling is particularly suited to the description of premixed flame combustion, which represents the main oxidation mechanism in spark ignition engines.It supposes that the chemical reaction of fuel oxidation occurs in a very thin layer.In premixed combustion, this layer separates burned and unburned gases and propagates toward the fresh mixture of fuel, oxygen and dilutant.
In what follows we present recent improvements made to the original CFM models as for example proposed in [4][5][6][7], yielding the ECFM model.

Flame Surface Density Transport Equation
The basis of CFM models is to describe the rate of fuel consumption per unit volume by the product of the flame surface density Σ (i.e. the flame surface per unit volume) and the local speed at which it consumes the mixture.
The flame surface density is determined via a model transport equation [8]: (1) essential in the case of high turbulence levels.Finally, a model to predict knock in SI engines is briefly described.These developments are then validated on two engine configurations: an optical access engine, for which LIF (laser induced fluorescence) measurements are available and the gasoline direct injection Mitsubishi engine for global validations.
where µ and µ t are the laminar and the turbulent viscosities respectively, and Sc and Sc t are the laminar and the turbulent Schmidt's numbers.The consumption rate for the fuel present in fresh gases is then computed as: (2) where ρu | u is the Reynolds' averaged density conditioned in the fresh gases and s l is the laminar flame speed. is the Favre's average of the fuel mass fraction conditioned in the fresh gases 1,2 .The definition of this quantity will be specified later.
Apart from the classical unsteady, convection and diffusion terms, this equation contains different source and sink terms that are modelled as follows: -P 1 = αK t represents the flame surface production by turbulent stretch, where K t is the ITNFS stretch [9] and α = 1.6 is a model constant; quantifies the production by the mean flow dilatation; models the effects of the flame thermal expansion and curvature; -is a destruction term due to consumption, where β = 1 is a model constant.-P k is a source term applied during the ignition period by the spark plug (see Section 1.5).

Computing the Fresh and Burned Gases States
In its former formulation [7] the model supposed that the fresh mixture does not contain exhaust-gas recirculation (EGR), and that the fresh gases temperature is simply estimated using a polytropic compression law.Since in stratified gasoline engines the composition is not uniform (the F/A equivalence ratio varying from very rich to very poor mixtures), Duclos and Zolver [10] introduced the fuel and oxygen tracers which allow to compute locally the fresh gases composition even in the presence of EGR.In the same way, they defined an equation for the fresh gases enthalpy which allows to compute precisely the fresh gases temperature.This major evolution of the CFM model was called ECFM.As these authors focused mainly on the application of this model to a practical industrial case, we describe in this section the complete equations of species and tracers and how the fresh gases properties are computed.
ECFM is based on a conditional averaging technique that allows to compute more accurately the fresh gases state, which is used to evaluate the local laminar flame speed.As the reaction rate of CFM-type models is proportional to the laminar flame speed, we expect that the ECFM improved fresh gases description will allow to better account for large scale stratification effect on combustion.Of course small scale stratification is not accounted for in this way, and will be addressed in Section 1.4.

Average Species Equations
In the ECFM model, we solve a transport equation for the average quantities of chemical species O 2 , N 2 , CO 2 , CO, H 2 , H 2 O, O, H, N, OH.This equation is classically modelled as: (3) where is the combustion source term and is the average mass fraction of species x: (4)

V is the cell volume and m
-is the cell mass.
The fuel is divided in two parts: the fuel present in the fresh gases, , and the fuel present in the burned gases, : (5) (resp. ) is the mass of the fuel contained in fresh gases (resp.burned gases).A transport equation is used to compute : (6) where , specified later, is the source term quantifying the fuel evaporation in fresh gases.The computation of will be presented in Section 1.2.5.

The Conditional Averaging Technique
In all CFM models, the flame front is described as an infinitely thin interface that separates fresh and burned gases.
In order to compute correctly the laminar flame speed in the fresh mixture and the pollutant formation in the burned gases, we need to know precisely the species mass fractions in these two regions, and .These mass fractions are defined as: where m -u (resp.m -b ) is the unburned (resp.burned) mass in the cell.These conditioned mass fractions are related to the average mass fraction Y ~x in the cell through the mass balance equation: where is the local burned mass fraction, also called the progress variable, defined as: (10) The reader could remark that quantities and are different: -compares the quantity of the species x, contained in fresh gases, to the cell mass (see Eq. ( 4)); -is the "real" average mass fraction of species x in fresh gases (see Eq. ( 7)).The expression "mass fraction of species x present in fresh gases and conditioned in these gases" is usually used to distinguish the quantity .Of course, the same difference, between and , can be noticed.
We have presented above the transport Equation (3) used to compute .We show in the next two sections how and can be calculated.Finally in the third section, we present how is deduced from Equation (9).

Definition of the Progress Variable
The CFM model is based on the assumption that the flame can be seen as an infinitely thin interface separating fresh and burned gases.As there can be no accumulation of mass within this interface, the conservation of mass through the flame front allows to say that the local burned mass fraction is proportional to the fraction of fuel mass that has been oxidised since the beginning of combustion: (11) where Y ~TF is the mass fraction of fuel before the onset of combustion.In the case of a perfectly mixed charge, Y ~TF is a constant in space and time in which case the calculation of is straitghforward.On the contrary, for a highly stratified charge engine, Y ~TF varies in space and time due to the imperfect mixing of the charge.In order to calculate Y ~TF exactly for this type of applications, we introduce an equation for the fuel tracer .By definition, the tracer of a species x is convected and diffused exactly like the real species it represents, but unlike the species x, this pseudospecies is not consumed during combustion nor accounted for in thermodynamic balances.Therefore its evolution equation is the same as the one for x without the reaction term: (12) It can be seen from Equation (11) that the progress variable is only defined if Y ~TF > 0. If this condition is not satisfied, it indicates that the cell never contained fuel, in which case combustion could not occur.Therefore in practice, the combustion model is not applied to a cell if Y ~TF is inferior to a prescribed small value ε.

Fresh Gases State
The fresh gases state has to be defined accurately to yield a correct estimation of the laminar flame characteristics and of the auto-ignition delay time for a correct prediction of knock, as shown in Section 1.6.

Fresh Gases Composition
We assume that fresh gases are only composed of fuel, O 2 , N 2 , CO 2 and H 2 O.
The fuel mass fraction in the fresh gases is simply given by: (13) This equality simply states that the fresh gases mixture composition is the same as the average composition we would get if no combustion had occurred.This average composition is exactly the one given by the fuel tracer.By adding an equation for the oxygen tracer Y ~TO 2 , we can apply the same argument for the mass fraction of oxygen in the fresh gases : The mass fraction of N 2 is computed by summing up the average mass fractions of all species containing nitrogen (i.e.N 2 , N, NO) weighted by the mass of nitrogen composing them: (15) where M x is the molar weight of species x.In case of the presence of residual gases the rest of the fresh gases is made of CO 2 and H 2 O.A balance of the carbon and hydrogen atoms yields the excess of these quantities compared to those issuing from the combustion of the initial fuel.This allows to compute their mass fractions in the fresh gases as: where n c is the number of carbon atoms in the fuel, and: where n H is the number of hydrogen atoms in the fuel.

Fresh Gases Temperature
To accurately estimate the fresh gases temperature T ~u we solve a transport equation for the fresh gases enthalpy which integrates compression work, heat losses due to walls and liquid fuel injection: (18) where Dissip is the heat generated by turbulent dissipation (taken equal to the average term for each cell) and the wall heat losses.As shown in [11] this heat losses term allows to precisely describe the fresh gases temperature in the end gas near walls and consequently to predict knock at the right location and instant.The compression work term (ρ -/ ρ -u ) ∂P ~/ ∂t accounts for the enthalpy increase in the fresh gases which differs from the average enthalpy increase by a factor ρ -/ ρ -u .
The spray source term (heat of evaporation) should be evaluated using the local conditions in the fresh gases.As it would need an important computational effort to perform this in addition to the modeling of the average terms, we choose to use the average spray source term in the fresh gases enthalpy equation.To account for the fact that this will lead to an over-estimation of the heat of evaporation 1 we propose the following approximation: (19) where is the average fuel mass source term due to the spray evaporation.Knowing the fresh gases composition one easily computes T ~u.

Laminar Flame Characteristics
From the knowledge of the fresh gases state the laminar flame speed is computed using the experimental correlation of Metghalchi and Keck [12]: (20) (1) During combustion, the fresh gases temperature is lower than the mean cell temperature.Thus, the amount of fuel that can evaporate will be lower than that estimated by using the mean temperature. with: T ~u and P ~u represent the fresh gases temperature and pressure (taken equal to the average pressure P ~), φ is the average local equivalence ratio in the fresh gases and is the local molar fraction of EGR gases.Coefficients B m , B φ and φ m are only known for propane, C 3 H 8 , and iso-octane, i-C 8 H 18 (see Table 1).A linear interpolation is used in order to compute these constants for a fuel C x H y : (21) Note that any other composition, for which these coefficients are known, could be used to extend this simple approach.

TABLE 1
Parameters B m , B φ and φ m for Eq. ( 20), [12] Fuel As the correlations of Metghalchi and Keck range for φ going from 0.7 to 1.4, we have no data for the laminar flame speed of very lean or very rich mixtures.As a first answer to the problem, we chose to extrapolate the correlations linearly to zero for φ = 0 and φ = 3.
Note that in the case of high stratification, we employ a more precise expression for the average flame speed (see Section 1.4.2) which is still based on the correlation of Metghalchi and Keck [12] but calculated for different equivalence ratios accounting for the fluctuations of this quantity in the fresh gases.

Burned Gases State
The burned gases composition and temperature is used in the model to compute some flame characteristics as well as pollutant production.The fuel mass fraction in the burned gases is computed using a transport equation.The other species in the burned gases and the burned gases enthalpy are rebuilt as shown below from the average and fresh gases state.No transport equation is solved for them.

Burned Gases Composition
In the case of very high local F/A equivalence ratios, there is not enough oxygen in the fresh gases to completely oxidise the fuel; therefore a part of the fuel going through the flame is not oxidised and ends up in the burned gases (where it can To account for this we introduce an equation for the fuel mass fraction located in the burned gases: (22) The source term represents the fuel going through the flame which is not.The source term due to the spray in Equation ( 22) is modelled as: (23) where is the average spray source term for the fuel.The equivalent source term for the fuel in the fresh gases (see Eq. ( 6)) is: (24) These source terms guaranty that the total source term of fuel mass is indeed .For all other species contained in burned gases, their mass fraction in these gases is given by Equation ( 9), which can be rewritten: (25) We note that for species O, OH, NO, and N, their mass fraction in the fresh gases is zero.
After updating the average state and by the ECFM results, an equilibrium computation is performed to correct the burned gases state and temperature.The set of equilibrium reactions used is that of the CO-CO 2 system taken from [13].This system is solved for based on the burned gases composition and temperature and allows to correct the estimation of the species mass fractions given by the turbulent combustion model.Finally, these new estimates of are reported to the average mass fractions Y ~x using Equation (9).

Burned Gases Temperature
Like for the species mass fractions in the burned gases, knowing the average and the fresh gases enthalpies allows us to directly compute the burned gases enthalpy: (26) Like for , the chemical equilibrium computation leads to a new estimate of h ~b which is reported to the average enthalpy h ~using Equation (26).
For small values of c ~Equation (26) can cause numerical problems.Therefore for c ~≤ ε we use instead the estimation: (27) where the enthalpy difference ∆H is computed as: (28) with hf x the heat of formation of species x and coefficients α 1 , α 2 and α 3 defined by the fuel composition and local F/A equivalence ratio.

Extension to Multi-Component Fuels
The fuel used in engines is not pure but made up of many components.Since these components have different properties, they evaporate at different speeds and consequently the composition of the fuel vapour within the combustion chamber is not uniform.On the other hand, in fuel droplets the lightest components tend to evaporate faster.So it is probable that droplets impinging on walls will be mainly made up of heavy components.These arguments show that we can expect a composition stratification in the chamber, which can be expected to have a direct influence on: -the combustion, since the laminar flame speed depends on the local concentration of the fuel vapour; -the knock tendency, as the auto-ignition tendency can differ from fuel to fuel; -the pollutants formation, since the composition of the vapour fuel could be characterised by a strong stratification.This is why we adapted the ECFM to multi-component fuels [14].As this paper is uniquely devoted to combustion modeling, we will not present here the description of the multi-component injection and film models developed at IFP, nor the definition of a multi-component fuel representing gasoline.We will only describe how the combustion of a multi-component fuel is modeled.
The multi-component combustion model is based on the definition of a "mean" local fuel which represents all the components of the gasoline.This mean fuel will be the only fuel species used during the combustion process.Once this process is accounted for, a reconstruction of the components of the fuel is performed in order to update the quantity of each component.These procedures are performed for each cell and at each time step.

Determination of the Mean Fuel
We define a fictitious mean fuel which is a purely local quantity: its composition is computed in each cell.The chemical formula of this fuel can vary from one cell to another.
Each fuel contained in fresh gases is represented by its molar concentration conditioned in these gases N -TF i and its formula C x i H y i .The molar concentration of the mean fuel F in fresh gases is simply N - . The number of carbon atoms x and hydrogen atoms y per molecule of mean fuel (defining the formula C x H y of the mean fuel) are simply obtained by summing the total number of carbon and hydrogen atoms for all the fuels: The computation of the molar weight M F and the heat of formation hf F of the mean fuel is performed using the same procedure: (30) where M F i (resp.hf F i ) is the molar weight (resp.heat of formation) of the component F i .

Computation of the Laminar Flame Speed
To model the influence of a multi-component fuel on the local flame propagation, we could ideally use correlations like the one used for iso-octane, which would take into account the proportion of the different components.Unfortunately, this type of correlation does not exist to our knowledge.This is why we chose a simple interpolation technique for the laminar flame speed, based on the local mean fuel composition.As in the mono-component case we use the correlation of Metghalchi and Keck [12].Constants B m , B φ and φ m of Equation ( 20) are computed using the linear interpolation (21).

Reconstruction of the Fuel Components
Once the reaction rate has been computed by the ECFM model, Equation (2), it is necessary to rebuild the components of the gasoline in order to follow the evolution of these components.Let us use the exponent n (resp.n+1) for all quantities computed at time step n (resp.n+1).If we note the concentration of the mean fuel consumed between the times n and n+1 by dN -F , we can write: (31) where: (32) The new value of the molar concentration of the component F i in the fresh gases is: (33) Following a similar procedure, it is possible to compute the molar concentration of each component F i in burned gases.

Modeling the Effects of Equivalence Ratio Fluctuations on Combustion
In the context of Reynolds averaged Navier-Stokes (RANS) calculations, the fresh gases properties only represent statistical average values of these quantities.In the case of a weak stratification, these average values correctly represent the real values observed locally during an engine cycle.In other words, the probability density function (PDF) of the F/A equivalence ratio, temperature, etc., is narrow.But when the stratification becomes stronger, which is the case of stratified charge engines, the local value observed during a cycle is not correctly represented by the average value anymore.A major consequence on combustion is that using the average equivalence ratio to calculate the laminar flame speed used in the ECFM model can lead to erroneous values.
In this case, it is necessary to estimate for each physical property its statistical distribution around its average value, which is commonly called the fluctuation of this property.As a first approach, we choose to neglect the temperature fluctuations in the fresh gases.On the contrary, the F/A equivalence ratio fluctuations are described using a twoequation model for the fuel mass fraction variance and dissipation which we describe in detail in another paper (see [15]).We present here how the F/A equivalence ratio fluctuation obtained from this model is used to compute a correct estimation of the average flame speed.

Variance and Scalar Dissipation Equations
The instantaneous mass fraction Y TF of the fuel tracer can be written as: where Y ~TF is the Favre's average of Y TF and its fluctuation.Note that in the multi-component case, Y ~TF represents the mass fraction of the sum of the different components.The scalar fluctuation variance υ is defined as: In [15] we present the procedure that allows to derive the unclosed equation of υ: where is the source term generated by the vaporisation of the injected fuel: is the evaporating source term in the transport equation of the fuel tracer (see Eq. ( 12)).As it can be For the turbulent transport (I) and the production (II) terms, a classical gradient assumption is adopted [15].The vaporisation source term (IV) in the variance equation contains four contributions (see Eq. ( 37)).The first two ones are unclosed correlations.To model them, we follow Demoulin and Borghi [16] assuming that the fuel mass fraction around an evaporating droplet k is equal to the mass fraction of the fuel at the saturating pressure vapour Y * k : (38) The above summations are performed over all spray particles contained in the current cell volume V with the rate of vaporisation on the surface of droplet number k.
The term (III) is the scalar dissipation χ which represents the main decreasing term in the variance equation (i.e.Eq. ( 36)).To close χ two different approaches can be chosen: 1) Algebraic closure [17][18][19] which assumes that χ is inversely proportional to the turbulent time scale: where k is the turbulent kinetic energy, ε its dissipation and C a constant of the order of 2. This closure is not well adapted to two-phase flows (flows occurring in gasoline direct injection engines) where the mixing scales during the vaporisation of the fuel are related to both droplet scales and turbulent scales.On the other hand, Beguier et al. [20] achieved an experimental study which showed that C is not always constant.2) Establishing a transport equation for χ.This kind of approach was adopted in several previous works, like [21][22][23][24], in order to ensure a more accurate description of χ.
We choose the second strategy.The reader could find all details about the derivation of the χ equation in [15].

Estimation of the Average Laminar Flame Speed
Supposing a local stratification means that a PDF (i.e.P(Y TF ) since , see Eq. ( 13)) can be defined for the fuel mass fraction in the fresh gases.In this case, the statistical average laminar flame speed is given by: (40) To perform this integration, the probability density function for the fuel mass fraction can be related to the average and variance through the use of a truncated Gaussian model PDF [14]: where the coefficients α, β and γ are computed by imposing Y ~TF (resp.u) as the average (resp.the variance) of P(Y).Note that this equation contains two Diracs δ at Y = 0 and Y = 1.
The integration of Equation ( 40) could be performed explicitly, but would lead to an important CPU overhead since it is necessary to perform it for each cell of the grid.We chose instead to perform a five point interpolation which proved accurate enough on tests calculations: (42)

The AKTIM Spark Ignition Model
In the calculations performed in [10] with the ECFM model a simple phenomenological model was used to describe the initiation of combustion by the spark plug.This model performed accurately in some simple homogeneous engine test cases, but it clearly showed a lack of predictivity when engine parameter variations on spark ignition were performed.The need to include phenomena like charge stratification, available electrical energy, heat losses to the spark plug and the influence of turbulence on the early flame kernel thus became clear and lead to the development of a new spark ignition model called AKTIM (arc and kernel tracking ignition model).This model has been described in [25] and we only give here a reminder of its general features and describe the most recent development realised at IFP.

General Description
AKTIM is based on four sub-models which describe realistically the different parts of the spark plug initiation: -The first sub-model describes the secondary side of the electrical inductive system.-The spark is represented by a set of particles regularly placed along the spark path, as shown in Figure 1.This spark is elongated in time by the mean flow field.-The flame kernels are described by Lagrangian marker particles.Each kernel can be seen as the initial flame development of one particular cycle of the engine.It is convected by the mean and turbulent flows and it receives energy from the electrical circuit and loses energy by contact with the spark plug.-Finally, a fourth sub-model describes the spark plug itself by a set of fixed discrete particles which are independent of the chamber mesh.This is why no remeshing for the plug is needed.
3D view of the spark plug, electrical arc and flame kernel centres.

Recent Improvements
In the model presented in [25] each flame kernel is seen as a sphere containing burned gases.The surface of this sphere consumes fresh gases exactly like in the ECFM formalism: where S i K is the sphere surface and the effective laminar flame speed on the surface (see [25] for a complete description of these quantities).This description of the fresh gases consumption implicitly assumes that the flame surface is not wrinkled by turbulence, which is a good approximation at the beginning of the flame development or when the turbulence intensity is low.
Yet, as shown in Section 2.2.2, in cases of high turbulence intensity on a high load operating point of a stratified charge engine, the first version of the model was shown to underestimate the consumption of fresh gases.We thus further improved AKTIM to enable it to represent the action of turbulence on the flame kernel, by introducing a turbulent wrinkling coefficient Ξ of the flame surface.
The laminar surface S i K in Equation ( 43) is replaced by the turbulent surface Ξ S i K , leading to an enhancement of the reaction rate as Ξ ≥ 1.Note that Ξ S i K is exactly the Lagrangian equivalent to the turbulent flame surface V Σ in the Eulerian description of ECFM.This is why we can define an evolution equation for Ξ by integrating the equation of Σ along the kernel radius.The final equation for Ξ reads: (44) where P 1 is the flame surface production term calculated like in the ECFM model, Equation (1), and D k is the destruction term equivalent to D in the ECFM model.The pressure term P 2 and the curvature term P 3 which appear in the Σ equation are respectively included in the description of the kernel growth through the evolution in time of the burned gases density inside the flame kernel and through the evolution of the ratio between the surface and the volume of the kernel.
In the laminar version of AKTIM presented in [25], when an arbitrary prescribed mass m 0 = 4πr 2 0 d (where d is the inter-electrode distance) of fresh gases has been consumed by the flame kernel, its surface is transferred to the equation of Σ.We now add another criterion for transition: if the wrinkling factor Ξ exceeds a prescribed value Ξ max of order 1 before the mass m 0 has been consumed by the kernel, the transition is operated.If this prescribed value is not reached during the kernel growth, the transition is still imposed when the total prescribed mass m 0 of fresh gases has been consumed.This new criterion guaranties that the flame shape remains approximately spherical, which is an hypothesis used to calculate the kernel characteristics.The total source term P k in the Σ equation becomes: where w k is the weighting factor of kernel k (inversely proportional to the total number of flame kernels used in the simulation) and V the volume of the cell receiving the surface.
With this development, the description of the initiation of combustion by AKTIM and of the following combustion process by ECFM are now nearly equivalent.As will be seen in the engine application results it allows the initiation process to be weakly dependant on the transition criteria.

The Knock Model
In order to predict knock in SI engines we use a simple model for knock prediction, compatible with the ECFM model.This model, described in [11], is based on a correlation for the auto-ignition delay θ which has been fitted using the data specified in [26].This correlation is a function of the local equivalence ratio, pressure, fresh gases temperature, EGR rate and octane number.This delay is then integrated in time using a fictitious species Y ~p which is initialised at zero: (46) where F(θ) is a quadratic function of θ ensuring that for constant fresh gases conditions, Y ~p equals Y ~TF after a time θ.The choice of a quadratic function was made to ensure that the auto-ignition delay calculation does not depend on the time at which the calculation starts.Ignition is defined to be reached when Y ~p = Y ~TF .As shown in [11] this simple model accurately predicts the auto-ignition delay in engine Flame kernel centers Spark plug applications, even though θ evolves in time-due to the evolution of pressure and fresh gases temperature-and space in the case of stratified charge engines where the equivalence ratio is not uniform.When auto-ignition starts, the fuel present in the fresh gases is consumed at a constant rate and a source of flame surface density is imposed in the Σ equation, corresponding to the development of a laminar flame: , (c -is defined as the volume fraction of burned gases).
In the case of multi-component fuels, the octane number used in the correlation has to be calculated locally as a function of the local fuel composition.We decided to use the following expression: where IO F i is the octane number of the component F i .

ENGINE APPLICATIONS
The ECFM model has already been validated on different engine applications, for instance [27] and [28].The aim of this section is to present first applications of the new models presented in this paper, keeping in mind that more complete validations will be presented in a near future.A first engine is used to test the multi-component, the knock and the scalar fluctuations models.A second one tests the development of AKTIM and the coupling between the scalar fluctuations model and the combustion model ECFM.

The GDI Optical Access Engine
To validate the multi-component combustion, the knock and the scalar fluctuation models, a comparison with experiments, achieved on a gasoline direct injection optical access engine [29], has been performed.This engine is a mono-cylinder engine based on a Renault four-cylinder PFI gasoline engine.Table 2 summarises the geometrical characteristics of the engine and the operating conditions where CAD stands for crank angle degrees and BTDC for before top dead centre.

Validation of Multi-Component Combustion and Knocking Models
The experiments we simulated were carried out on the engine using a fuel composed (in mass) of 35% of n-nonane, 45% of iso-octane and 20% of n-pentane.Using Equation (47) and Table 3 yields a low octane number for the mixture of , indicating that auto-ignition might occur in the case studied, which is confirmed by the experimental results.TABLE 3 Octane number of the used fuels, [26]

Fuel
Octane number n-nonane -10 iso-octane 100 n-pentane 62 This case is very interesting to calculate since it allows to validate the multi-component combustion and knock models at the same time.Figure 2 compares the cylinder pressure evolution predicted by the ECFM model with experimental findings.The pressure curve, obtained without the knock model, clearly shows a bad combustion which leads to an under-estimation of the pressure levels.Activating the knock model in the calculation leads to auto-ignition and clearly improves the predictions.On the other hand, we can say that with this kind of multi-component fuel, the engine runs almost under controlled auto-ignition (CAI) mode.The multi-component knock model coupled with ECFM accurately predicts the fast heat release generated by the onset of auto-ignition near the top dead centre observed experimentally: the burned fuel mass fraction is well predicted.

Validation of the Scalar Fluctuation Model
The engine also served in [15] to validate the scalar fluctuation model.Hereafter we only present an example of the obtained results to illustrate the model performance in predicting charge stratification and composition fluctuations.The case presented corresponds to an injection timing of 23.6 CAD BTDC, leading to a highly stratified charge.For this case the fuel is pure iso-octane and no combustion is simulated, as our goal is to test the mixture preparation before combustion.
Using the LIF visualisation technique, the experiments provide a 2D view of the instantaneous equivalence ratio in the vertical plane defined in Figure 3.The injector and the spark are located on this plane: the symbol "×" (resp."↓"), on the experimental fields of Figure 5, corresponds to the spark (resp.injector) position.The images are the result of an averaging of 300 engine cycles.Four visualisation instants were acquired: 11, 9, 7 and 5 CAD BTDC.Experimental measurements allow to compute the average equivalence ratio and its standard deviation using the classical ensemble averages.We present in [15] how we compute the corresponding quantities in our calculations starting from the average mass fraction Y ~F of the fuel and the variance υ.We also show that due to the discrete acquisition performed on the visualisation plane by the CCD camera in the experiments, we need to filter the variance υ obtained in the calculations in order to compare with the experiments.This filtering is realised thanks to an original procedure also presented in [15].
The comparisons between computations and experiments were performed in two different ways.We first compared the F/A equivalence ratio fluctuation in the spark zone (note Figure 3 A general view of the GDI optical access engine mesh, with the location of the LIF visualisation plane. that a particular spark is used, see reference [29]).For the experimental measurements, this zone is defined as a disk of 2 mm radius, centred on the electrodes gap and located in the vertical measurement plane.In the computations this zone is represented by a sphere of the same radius.Then we compared the fields of equivalence ratio fluctuation in the experimental measurement plane at the four instants of acquisition.
Figure 4 shows the evolution of the average and the fluctuation of the equivalence ratio at the spark location.The high level of fluctuations observed during the spray vaporisation and the following rapid decrease due to turbulent mixing are accurately reproduced by the scalar fluctuation model.

The Mitsubishi GDI Engine
The aim of the simulations of the Mitsubishi GDI engine is to check the good behaviour of the ignition model AKTIM and of the combustion model ECFM when the variance/scalar dissipation model is used.We are not aiming here at reproducing precisely the experimental results; for this purpose a more complete validation with experiment will be the subject of future work.

Engine Description
The engine used in our calculations is the same as the one presented in [10]: it is the first 1.8 l GDI Mitsubishi engine that has been used in a passenger car.Only half of the chamber is meshed as there exists a symmetry plane between the four exhaust and intake valves, see Figure 6.Table 4 shows the two simulated engine test cases.
A first operating point at high load is used to validate the turbulent wrinkling description in the AKTIM model (see Section 1.5).This point is based on an experimental working point with nearly homogeneous conditions (injection in the intake stroke) in an over-charged adaptation of the engine with a late spark timing imposed to avoid knocking.This point was chosen because it presents high levels of turbulence and is also a difficult test case for previous ignition models.The second case parameters are representative of a late injection strategy.This test is aimed at reproducing high levels of equivalence ratio fluctuations during the

Validation of the Ignition Model
We present in Figure 7 the pressure evolution for the different calculations (see Table 5) of case A using the AKTIM model for ignition with the ECFM model.When the wrinkling factor Ξ is not used in the AKTIM model (cases A1 and A2), it is clear that the initiation of combustion is too slow, leading to an underestimated pressure.Besides when the critical mass m 0 (or radius r 0 ) is changed, the pressure profile changes strongly.These results indicate that unlike in the test cases presented in [25], the ignition model is no more independent of m 0 for this case.
The reason for this is that for large values of m 0 (case A2), the flame kernel lives long enough to be wrinkled by the strong turbulence of case A. But as we suppose the flame to be laminar, we underestimate the combustion intensity.On the contrary, for small values of m 0 (case A1), the kernel rapidly transfers its surface to the ECFM model which correctly describes the flame wrinkling by turbulence (through the Σ equation).
Pressure evolution for case A on the GDI Mitsubishi engine.As can be seen in Figure 8, the factor Ξ increases regularly from 1 (which corresponds to the laminar flame) to 1.2 for case A3 after approximately 1 CAD, and to 2 for case A4 after 3.2 CAD.Also displayed in this figure is the production term P 1 and the destruction term D k normalised by the curvature strain which represents the increase of flame surface due to curvature effects.It can be seen that initially the turbulent strain P 1 is negligible compared to the curvature strain E l , but as the flame radius increases, the curvature term tends to zero (the flame becomes planar) and the turbulent term becomes dominant.For cases A3 and A4, the turbulent wrinkling Ξ is too low to generate a high level of surface destruction D k .
Introducing the factor Ξ in the AKTIM model intensifies the combustion during the flame kernel life, which leads to a faster transition to the 3D combustion model and to a higher source term P k (factor Ξ in Eq. ( 45)).As can be seen in Figure 7, cases A3 and A4 exhibit comparable pressure curves, unlike cases A1 and A2, indicating that adding a turbulent description to the flame kernels allows to use longer initiation periods.

Validation of the Fluctuations Model with Combustion
We present in Table 6 the calculations of case B. For case B1 no variance/dissipation model is used and therefore the laminar flame speed used in the ECFM model is defined by Equation (20).On the contrary, for cases B2 and B3, the variance is used to compute the laminar flame speed as in Equation (42).For case B2 we use the fluctuations model briefly presented in Section 1.4.1.For case B3, the same model is used but the scalar dissipation is artificially increased.This way, the variance destruction by scalar dissipation will be higher for case B3 than for case B2.This is what we can observe in Figure 9a where the average and the fluctuation of the equivalence ratio at the spark plug are shown for cases B2 and B3.Note that the average equivalence ratio is the same for the three calculations before combustion onset as the variance model does not influence the average quantities during this period.In Figure 9b we observe that the pressure increase is much faster for case B1 than for cases B2 and B3.This result indicates that defining the flame speed at the average equivalence ratio conditions can induce an important error in the average flame speed estimation in the case of high fluctuations levels.Comparing cases B2 and B3 clearly shows that as the fluctuation level decreases, the average laminar flame speed increases due to the smaller importance given in Equation (42) to the very rich or very poor mixtures which have very low laminar flame speeds.These curves also show that when the fluctuation level tends to zero, the ECFM model with the variance/scalar dissipation description tends toward the ECFM model of case B1.

CONCLUSIONS
This paper is a summary of a number of recent model developments performed at IFP and concerning the adaptation of a CFM turbulent combustion model to the simulation of DI-SI engines, yielding the ECFM model.In order to address the high engine stratification in such engines, and to include more realistic fuels, the following model developments were assessed in the first part of this paper: -Improving the description of large scale stratification on combustion: the local laminar flame speed is a crucial component of CFM type models, strongly depending on local fresh gases equivalence ratio, temperature and EGR fraction.We presented how the ECFM model allows to accurately compute these quantities thanks to the conditional averaging technique, to the use of species tracers and of a fresh gases enthalpy equation.-Modeling the combustion of multi-component fuels: we presented an adaptation of ECFM to be able to use a more realistic fuel which could contain many components having different properties in order to better describe the local flame speed, octane number and the possible formation of HC. -Accounting for equivalence ratio fluctuations around the average: a scalar fluctuation model was presented that allows to account for small scale stratification on combustion through a simplified integration of the laminar flame speed over fuel mass fraction using a presumed truncated Gaussian PDF for its distribution.-Improving the predictivity of the AKTIM spark ignition model by including turbulence effects on the early flame kernel growth.
-A simple knock model coupled to the ECFM model and based on a simple ignition delay correlation was presented in the context of multi-component fuels.
The second part of the paper was devoted to a presentation of test simulations performed using two different engine geometries.
The simulation of a DI-SI optical access gasoline engine was used to illustrate the validation of the multi-component, knock and scalar fluctuation sub-models of ECFM.A first series of simulations allowed to show that ECFM accurately reproduced the experimental evolution of the cylinder pressure in the case of auto-ignition with a multi-component fuel.A second series served to validate the scalar fluctuation model in the case of a mono-component fuel without combustion.It was shown that the model correctly reproduces the experimental evolution of equivalence ratio at the spark location, as well as the qualitative and quantitative spatial average and fluctuating fuel distribution in the chamber.
The simulation of a homogeneous, high turbulence case of a Mitsubishi GDI based engine allowed to check the correct behaviour of the AKTIM spark ignition model on an overcharged configuration.Finally the simulation of a partial load operating condition (high stratification obtained by late injection) of the same engine showed that the proposed scalar fluctuations model reduces the rate of combustion when the stratification is high.
Future work will be devoted to more extensive validations of the models presented in this paper, making extensive use of experimental measurements performed at IFP on industrial type DI-SI engines.In particular, the focus will be on the prediction of cycle-to-cycle instability by modeling the effects of equivalence ratio fluctuations on spark ignition with the AKTIM model.Another main research path will be the improvement of the prediction of pollutants such as carbon monoxide and soot.Furthermore, these models have already been used in different types of engines from GDI engines to formula one engines.
(36)  contains several terms to be closed.

Figure 2
Figure 2Experimental averaged and computed pressures in the combustion chamber of the GDI optical access engine.

Figure 4
Figure 4 Evolution of the average and the fluctuation of the F/A equivalence ratio at the spark location.

Figure 5
Figure5shows the fields of the equivalence ratio fluctuation for experiments (on the left) and calculation (on the right).Only the part of the computed fields located inside the frame drawn on these fields have to be considered.This frame represents the observation window through which the

Figure 6 1
Figure 6 1.8lMitsubishi GDI engine: fuel/air equivalence ratio at 27 CAD BTDC for case B.
A3: Ξ max = 1.2 A4: Ξ max = 2 Figure 8Evolution (minimum and maximum values over all kernels) of the flame kernel wrinkling factor Ξ, the production term P 1 and destruction term D k for calculations in case A with the AKTIM model including the turbulent description.
E l min.P 1 /E l max.D k /E l min.D k /E l max.Ξ max.= 1.2 (case A3) E l min.P 1 /E l max.D k /E l min.D k /E l max.Ξ max.= 2 (case A4)

Figure 9
Figure 9Results for Case B: a) Evolution of average and fluctuation of the equivalence ratio at the spark plug; b) Evolution of the cylinder pressure.

TABLE 4
Definition of the two cases on the GDI Mitsubishi engine