A New 0D Diesel HCCI Combustion Model Derived from a 3D CFD Approach with Detailed Tabulated Chemistry

A New 0D Diesel HCCI Combustion Model Derived from a 3D CFD Approach with Detailed Tabulated Chemistry — This paper presents a new 0D phenomenological approach to the numerical modelling of Diesel HCCI combustion. The model is obtained through the reduction of 260 Oil & Gas Science and Technology – Rev. IFP, Vol. 64 (2009), No. 3 TKI-PDF (Tabulated Kinetics for Ignition, coupled with presumed Probability Density Function) 3D CFD model developed at the IFP. Its formulation is based on physical considerations, to take into account the main phenomena and their mutual interactions that take place in the cylinder during the combustion process. Aspects relating to spray penetration, fuel evaporation, turbulence, mixture formation and chemical kinetics have been studied in detail. The original contribution of this work concerns the modelling of the formation and evolution of the equivalence ratio stratification around the spray, and of its connection to combustion kinetics. In order to achieve this, different tools commonly adopted in 3D modelling have been adapted to 0D modelling. Presumed PDF theory has been extended to a 0D formalism in order to characterize the mixture-fraction distribution. This approach has then been coupled with droplet-evaporation theory in order to have access to the thermodynamic conditions characterizing the mixture. The temporal evolution of the spray is computed in terms of volume and the entrained mass of gases, starting from conservation laws for mass, momentum and energy. An adapted κ − model is used to take into account the turbulence in the cylinder, which is very important, in an ICE (Internal Combustion Engine), especially during the mixing process. Further, combustion heatrelease is computed using an adapted detailed tabulated chemistry method inspired by the FPI (Flame Prolongation of ILDM (Intrinsic Low Dimensional Manifold)) theory. This look-up table allows the simulation of a large range of combustion regimes, since it takes into account the presence of EGR (Exhaust Gas Recirculation) in the mixture. The results of the 0D model are compared in an initial step to the 3D CFD results. Finally, the OD model is validated against a wide experimental database. NOMENCLATURE Latin abbreviations a [−] constant coefficient (a = 0.66) b [−] constant coefficient (b = 0.41) B [−] transfer number at steady-state Bm [−] mass transfer number Bt [−] thermal transfer number c [−] progress variable Ca [−] orifice area-contraction coefficient Cd [−] orifice discharge coefficient Cdiss [−] Z̃”2 dissipation-term coefficient Cev [−] evaporation-rate coefficient CK [−] K̃ dissipation-term coefficient Cp [J/kg/K] constant pressure specific heat Cv [J/kg/K] constant volume specific heat Cθ [−] spray opening-angle coefficient Cκ [−] κ̃ dissipation-term coefficient dh [m] injector hole diameter D [m] droplet diameter Ec [J] kinetic energy h [J/kg] enthalpy k [J/m/s/K] thermal conductivity K̃ [J/kg] mean specific kinetic energy L [m] liquid-phase fuel penetration Lv [J/kg] latent heat of vaporization Le [−] Lewis number m [kg] mass ma [kg] stoichiometric mass of air mcyl [kg] total in-cylinder gaseous mass M [kg/mol] molar mass n [mol] mole number nh [m] number of injector holes NK [−] K̃ dissipation-term exponent Nκ [−] κ̃ dissipation-term exponent p [Pa] pressure pcr [Pa] critical pressure pv [Pa] fuel vapor pressure P [−] probability function Pco [−] coupled probability function q [−] constant coefficient (q = 2.2) Qev [J/s] fuel evaporation thermal flux Qhu [J/s] fuel heat up thermal flux Qth [J/s] exchanged thermal flux Qtot [J/s] total thermal flux r [m] radius R [J/mol/K] universal gas constant s [−] mass stoichiometric coefficient S [m] spray penetration t [s] spray characteristic time-scale tin j [s] relative injection time T [K] temperature Tb [K] normal boiling temperature Tcr [K] critical temperature UFl [m/s] fuel injection velocity v [m3/mol] molar volume V [m3] cylinder volume VS [m3] spray volume x [m] spray characteristic length-scale x0 [m] removed cone height X [−] molar fraction Y [−] mass fraction Ylc [−] mass fraction linear combination Z [−] mixture fraction Z̃”2 [−] mixture fraction variance A Dulbecco et al. / A New 0D Diesel HCCI Combustion Model Derived from a 3D CFD Approach... 261 Greek abbreviations δ [−] segregation factor ̃ [J/kg/s] κ̃ dissipation rate θ [rad] spray opening angle κ̃ [J/kg] specific turbulent kinetic energy λ [m2/s] evaporation constant ν [mol] stoichiometric molar coefficient ρ [kg/m3] density σ [−] weighting coefficient τev [s] characteristic time of evaporation τt [s] turbulence characteristic time φ [−] Schwab-Zeldovitch variable Φ [−] equivalence ratio ψ [−] Dirac function ω [kg/s] reaction rate ΩS RK [−] acentric factor Superscripts S/I higher/lower Φ tabulated values 0 Φ reference value Subscripts A ambient gas c complementary mixture-fraction ds dissipation term eq chemical equilibrium conditions esp species F gaseous fuel Fin j injected fuel Fl liquid fuel g equivalent mixture i/o inlet/outlet m mass diffusion boundary layer max maximum value mix mixture O oxidizer pr production term r EGR s liquid-gas interface S spray S/I higher/lower Xr tabulated values t thermal diffusion boundary layer w withdrawn ambient gas β − PDF normalized β − PDF variables τ tracer (non reactive) conditions 0 reference value 1/3 1/3 law reference variable ∞ condition at infinity


TKI-PDF (Tabulated Kinetics for Ignition, coupled with presumed Probability Density Function) 3D
CFD model developed at the IFP. Its formulation is based on physical considerations, to take into account the main phenomena and their mutual interactions that take place in the cylinder during the combustion process. Aspects relating to spray penetration, fuel evaporation, turbulence, mixture formation and chemical kinetics have been studied in detail. The original contribution of this work concerns the modelling of the formation and evolution of the equivalence ratio stratification around the spray, and of its connection to combustion kinetics. In order to achieve this, different tools commonly adopted in 3D modelling have been adapted to 0D modelling. Presumed PDF theory has been extended to a 0D formalism in order to characterize the mixture-fraction distribution. This approach has then been coupled with droplet-evaporation theory in order to have access to the thermodynamic conditions characterizing the mixture. The temporal evolution of the spray is computed in terms of volume and the entrained mass of gases, starting from conservation laws for mass, momentum and energy. An adapted κ − model is used to take into account the turbulence in the cylinder, which is very important, in an ICE (Internal Combustion Engine), especially during the mixing process. Further, combustion heatrelease is computed using an adapted detailed tabulated chemistry method inspired by the FPI (Flame Prolongation of ILDM (Intrinsic Low Dimensional Manifold)) theory. This look-up

INTRODUCTION
Numerical simulation plays an increasingly important role in engine development processes; its most important advantage being due to the fact that, it is possible to perform a great number of numerical tests, at costs that are much lower than those associated with experiments. This becomes particularly important with future generations of ICE. In fact, in order to obtain best performance, engines are increasingly composed of complex system-layouts involving a great number of devices such as an EGR system, a CR (Common Rail) Direct Injection System, etc. Such systems need to be optimized [1][2][3][4].
In this context, a numerical engine simulator must be able to give an accurate estimate of the main values characterizing the performance of the engine, such as combustion heat release, combustion efficiency, and pollutant formation.
This paper presents a new 0D approach to the modelling of Diesel HCCI combustion derived from the reduction of a 3D CFD model. The new model is able to accurately predict the combustion process, while requiring low CPU resources. These two aspects outline a model profile adapted to globalsystem simulation applications. In the following, a global overview of the model is given first. The two main parts of the model are subsequently presented and detailed: the spray and the combustion modelling. Partial validations of these models have been performed using 3D results and experimental data. Finally, validation of the global combustion model is done by comparison with experimental results. The paper ends with a discussion of the model improvements and perspectives.

GLOBAL OVERVIEW OF THE MODEL
Conventional Diesel and HCCI combustion modes are governed by different physical mechanisms.
In conventional Diesel combustion, the mixture inside the cylinder is characterized by a high fuel mass fraction stratification. The first site of auto-ignition appears inside the spray, where chemical and thermodynamic conditions are the most favourable. The most important challenge concerning auto-ignition phase description is to capture the shortest ignition delay. Heat release, due to combustion at the first site, subsequently favors the multiplication of auto-ignition sites. This chain reaction leads to the sharp heat release process typical of the premixed Diesel combustion phase. The remaining fuel burns in a diffusion flame and the heat release is governed by the mixing process. In this phase, turbulence plays the most important role.
In HCCI combustion, the mixture inside the cylinder can be considered as homogeneous. Hence chemical kinetics plays the most important role in the combustion process. Because of the homogeneity of the mixture, all of the gas auto-ignites simultaneously. As a consequence, HCCI combustion is characterized by a stiff heat release process. For these reasons, a deep understanding of the spray formation mechanic is necessary [39][40][41][42][43][44].
In a general situation, the two combustion regimes are not so well distinguished. A modelling challenge, therefore, is to be able to manage all of the different scenarios. Moreover, the presence of EGR in the cylinder modifies the chemical kinetics, and phenomena, such as cold flames, could appear at the beginning of the fuel oxidation process.
The original modelling solution proposed in this work is derived from a reduction of the TKI-PDF 3D CFD model [33,34,45]. It has the potential to take into account all of the different phenomena, taking advantage of the common tools used in sub-grid 3D modelling. The first of these is the use of a presumed PDF [12,13,46] in order to describe the distribution of fuel inside the spray region. The second is the use of a look-up table in order to take chemical kinetics into account This technique has been inspired by the FPI theory [47]. The third tool is the use of an adapted κ − model to simulate the impact of turbulence on the mixing and combustion processes.
The synoptical diagram of the complete Diesel HCCI combustion model is presented in Figure 1. As shown, the c c Figure 1 Synoptical diagram of the complete Diesel HCCI combustion model. model is composed mainly of two parts, namely spray modelling and combustion modelling (these will be described in further detail in Sect. 2 and 3 respectively).
The purpose of spray modelling is to represent the formation and the temporal evolution of the spray volume (mixture zone), and to give constant access to the total spray mass and fuel mass-fraction distribution inside the spray. As shown in Figure 1, this model has only one input variable: the fuel injection-rate. The mixture zone has two mass sources: one is the evaporation from liquid fuel injected into the cylinder, and the other is the entrained fresh gas (in the most general case, it is a perfectly stirred mixture of air and EGR). The mixing process of gaseous fuel and fresh gas takes place in the spray volume and is controlled by turbulence. As shown in the diagram, the turbulence submodel has two inputs: one is the kinetic energy contribution associated with the injected mass of fuel and the other is the kinetic energy contribution associated with the entrained fresh gas. The mixture stratification within the spray volume is represented by a presumed PDF. The temporal evolution of the PDF shape is governed by the mass flow-rates concerning the mixture zone and the turbulence intensity. The output of the spray model is the complete spray description. It consists of quantifying the total gaseous mass contained in the mixing zone, and characterizing the fuel mass-fraction distribution inside it as well as the associated temperature field.
The goal of the combustion model is to describe the chemical composition of the mixture and the oxidation process as a function of the in-cylinder thermodynamic conditions. As shown in Figure 1, the only input required for the combustion model is the spray characterization given by the spray model. The outputs of the combustion model are the heat release and the temporal evolution of the species masses during combustion.

SPRAY MODELLING
As shown in Figure 1, the spray model contains four submodels (directional arrows indicate the different interactions): -The evaporation submodel (Sect. 2.1): this model computes the thermodynamic conditions at the liquid-gas interface during the evaporation process, and quantifies the liquid-fuel-mass evaporation-rate. It also gives qualitative information regarding the penetration length of the liquid fuel jet. -The gas entrainment submodel (Sect. 2.2): this model computes the entrainment of fresh air in the mixture zone. It also gives quantitative information regarding the penetration length and the volume growth of the spray. -The turbulence submodel (Sect. 2.3): the goal of this model is to compute the characteristic frequency of the turbulence, which is directly connected with the turbulent mixing process of gaseous fuel and fresh gas (air and eventually EGR). As well known in ICE combustion, molecular diffusion of species is negligible compared to the turbulent one: for this reason molecular diffusion will not be considered [48]. -The mixture submodel (Sect. 2.4): the role of this model is to quantify the total mass in the mixture zone and to characterize the temporal evolution of the fuel mass fraction distribution.

Evaporation Submodel
The evaporation submodel has two main purposes: firstly, to compute the thermodynamic states of the fluids at the liquid-gas interface, and secondly to quantify the mass evaporation rate of the injected liquid fuel. The evaporation of droplets in a spray simultaneously involves heat and mass transfer processes. The overall rate of evaporation depends on many factors, such as pressure, temperature and transport properties of fluids, drop geometry and relative velocity of drops respective to the ambient gas. In the case of a droplet immersed in a gas at a higher temperature, part of the heat flux increases the droplet temperature, Q hu , and the rest evaporates it. At steady-state, the droplet temperature reaches its wet-bulb value and all the heat is used for the evaporation process. When steady state conditions are established, the drop diameter, D, diminishes with time according to the well-established D 2 law: where λ is the evaporation constant. The thermodynamic interfacial conditions are determined using the approach proposed in [49]. It describes the evaporation process, in steady-state conditions, of a single spherical liquid drop of pure fuel at rest in a gaseous mixture with a given composition, Figure 2. In such conditions, it is possible to define a mass transfer number, B m , and a thermal transfer number, B t , characterizing the mass and thermal boundary layers respectively, as follows: where Y F s and Y F ∞ are the fuel mass fractions of the mixture at the interface and ambient gas conditions respectively. T s and T ∞ are the temperature of the mixture at the interface and ambient gas conditions respectively, and L v s is the latent vaporization heat of the liquid at the interface temperature, computed using the Watson relation [50]. C p g is the equivalent gas mixture heat-capacity at constant pressure computed using the 1/3 law proposed by Sparrow and Gregg [49]. According to this law, average properties of the mixture are evaluated at reference temperature and composition defined as: Consequently, the reference ambient-gas mass-fraction, Y A 1/3 , is given by the relation: At infinite distance from the drop, fuel concentration is assumed to be zero. C p g is then obtained as: where C p A and C p F are the heat-capacities of the ambient gas and of the gaseous fuel respectively. According to [49], steady-state conditions hold: The system is solved iteratively varying the value of T s in order to verify Equation (6). With a maximum absolute error of 1e −5 , the number of iterations necessary to satisfy Equation (6) varies around 18. The link between the two transfer numbers is detailed as follows. At a given fuel saturation temperature at the liquid-gas interface, T s , the fuel vapor pressure, p v s , is estimated from the Clausius-Clapeyron equation as: in which: where T b is the normal boiling temperature of the fuel, and p cr and T cr represent its critical pressure and temperature respectively. Fuel mass fraction at the interface can now be expressed as: where p ∞ is the ambient pressure and, M A and M F are molar masses of ambient gas and fuel respectively.
According to [49], the fuel mass flow rate,ṁ F , can be computed as:ṁ and the mass of fuel, m F , constituting a drop of given diameter is: where k g is the mixture thermal conductivity, computed using the 1/3 law, Le is the Lewis number (here supposed equal to unity) and ρ F l is the density of the liquid fuel estimated using the Hankinson-Brobst-Thomson (HBT) relation [50]. HBT relation holds to: in which v s is the molar volume of the liquid at the liquid-gas interface temperature, that can be computed as: where: and Ω S RK is the acentric factor of the fuel. Moreover, the constant coefficients (that are valid for all hydrocarbons) assume the following values: HBT correlation is valid in a temperature domain defined as: By differentiating Equation (11) and combining it with (1) and (10), it is possible to express the evaporation constant as: Figure 3 shows the dependence of the evaporation constant on the ambient pressure and temperature conditions for n-heptane. As shown, there is good agreement between the values computed by the model and experimental data [51]. Figure 4 shows the influence of ambient thermodynamic conditions on the mixture temperature at the liquid-gas     Dependence of the temperature at the liquid-gas interface on the ambient conditions. interface. The computed thermochemical properties of the mixture at liquid-gas interface, such as Y F s and T s , will be used in the mixture characterization model and in the combustion model. These models will be presented in Sections 2.4 and 3.3 respectively. In particular, T s is a very important variable in the spray modelling, as it permits an estimate of the local temperature of the mixture. In what concerns the computation of the global fuel-mass evaporation-rate, several aspects must be highlighted. Equation (10) refers to single-drop evaporation. The real spray scenario is very different as droplet size distribution, ligaments and the undefined geometry of the liquid phase make it much more complicated. Consequently, because of its CPU cost, it can not be solved in a global-system simulation approach.
For that reason, the global liquid-mass evaporation-rate has been estimated using a simpler approach, based on a characteristic time of evaporation, τ ev : where C ev is a dimensionless coefficient and m F l is the liquid mass of fuel. In Section 4.1 the results of this model are compared with those of the 3D CFD code. The original method used to estimate the evaporation characteristic-time is an approach based on the theory proposed in [52,53]. This theory permits an estimate to be made of the liquid penetration length in a pressurized vessel as a function of the injector geometry, and of the compositions of the fluids and their thermodynamic states. In this approach, the spatial, x + , and temporal, t + , scales are defined as: where C a is the orifice coefficient of area-contraction, d h is the injector hole diameter, ρ A is the ambient gas density, a is a dimensionless coefficient (authors recommend a = 0.66), θ is the opening angle of the spray (its definition will be given in Sect. 2.2), and U F l is the fuel velocity in the injector nozzle. The theoretical expression for the maximum liquid penetration, L, is given as: where b is a dimensionless coefficient (authors recommend b = 0.41) and B is the transfer number computed using Equation (6). τ ev is computed by inverting the equation for the jet penetration length as a function of time and imposing a penetration length equal to L [13]. That means that τ ev represents the time necessary to fully evaporate a droplet whose lifetime corresponds to a penetration equal to L, once the velocity law is known. The expression of τ ev is: in which:L Figure 5 Spray volume description. Although the presence of two distinct volumes, the fresh-air and the spray volumes, the model is a one-zone approach. Hence, only one energy balanceequation is computed for all the in-cylinder gases.

Gas Entrainment Submodel
This submodel is based on the well known theory developed in [54] for estimating spray penetration and dispersion. The spray region is supposed to be a perfect truncated cone, Figure 5. Its volume, V S , is completely defined by the height (spray penetration, S ), and the angle (spray angle, θ).
These two variables depend on the nozzle geometry, the injection pressure, and the injected and ambient fluid thermodynamic states. They are determined from the following semi-empirical correlations: and: where q is a dimensionless coefficient (authors recommend q = 2.2), t in j is the relative time coordinate, with origin at the beginning of the injection, and C θ is a constant coefficient. Hence: with x 0 defined in Figure 5. Through experimentation of injection in a constant volume vessel containing a gas at a given pressure and temperature, these correlations have been widely validated. Results are reported in various works by different authors [20,42,54,55]. Vapor phase penetration. n-heptane, T ∞ = 800 K, ρ A = 30 kg/m 3 .
Figures 6 and 7 compare the computed spray penetrations with experimental data obtained in a high pressure vessel [42]. The results are in good agreement with experiments. As can be seen, the model is sensitive to the variation of parameters such as injection pressure, nozzle geometry and the thermodynamic state of the ambient gas.
In ICE conditions, in-cylinder thermodynamic conditions vary in time because of the piston motion, thermal exchanges and combustion process. In order to take into account the time history of spray penetration, the differential form of Equation (27) has been implemented in the model. Consequently, the spray volume is computed as: The mass of entrained gas in the spray volume is deduced from geometrical considerations: (29) where dm A /dt is the entrained mass-flow-rate and ρ A is the density of the gas still outside the spray region.

Turbulence Submodel
In ICE, turbulence plays a very important role, increasing mass, momentum and energy transfer-rates. The dynamics of turbulent motion can be described as follows. At first, kinetic energy associated with the mean flow is transferred to the large scales of the fluctuating motion (integral scale). After that, kinetic energy is transferred isentropically from large to small scales of turbulence (Kolmogorov scale). Finally, the kinetic energy is dissipated by viscous forces and is transformed into thermal energy. The aim of the turbulence submodel is to obtain an estimate of the turbulence characteristic time, τ t , which is very important in the mixing process modelling. The proposed approach is based on the well-known κ − turbulence-model theory. According to this theory, the kinetic energy is decomposed into two contributions: -the mean kinetic energy, associated to the big structures of the mixture flow, -the turbulent kinetic energy, associated to the local mixture velocity-fluctuations respect to the mean value. The procedure to estimate τ t is detailed in the following. The turbulence submodel constitutes of two state variables: the spray global kinetic-energy, E c , and the mean specific kinetic-energy,K (1) .

Spray Global Kinetic-Energy
The spray global kinetic-energy rate can be most generally expressed as the sum of a production and a dissipation term: The production term in an ICE application can be detailed as follows: The spray energy, E c S , represents 98% of the total kinetic energy, according to [20]. For that reason, with good approximation, it is possible to affirm that: (1) The symbol tilde means that the mean kinetic energy value has been computed using the Favre mean operator.
in which: (33) where m F in j is the mass of injected fuel, C d is the discharge coefficient of the orifice and n h is the number of holes of the nozzle.
The spray global kinetic-energy dissipation-term, mainly due to flow viscous forces, has been modeled by using empirical closure terms [13].
Equation (30) can now be rewritten as: where m S is the total mass of the spray, C K and N K are empirical dissipation-term constant-coefficients associated to the mean specific kinetic-energy,K, and C κ and N κ are empirical dissipation-term constant-coefficients associated to the mean specific turbulent kinetic-energy,κ.

Mean Specific Kinetic-Energy
Similarly to κ − theory, by spatially filtering the 3D momentum balance-equation, it is possible to obtain an expression for the temporal evolution of the mean specific kinetic-energy [13]: where, as can be seen, the last right-hand-side term represents the mean specific kinetic-energy dissipation-term, already present in Equation (34). The filter size used to obtain Equation (35) is the spray volume, V S .

Characteristic time of turbulence
It is now possible to compute the mean specific turbulent kinetic-energy value as:κ The term indicating the mean specific turbulent-kineticenergy dissipation-rate,˜ , is supposed to have an expression similar to the one used to represent the dissipation rate of the mean specific kinetic-energy: as shown in Equation (34).
According to the κ − turbulent model theory, the characteristic time associated to turbulence is defined as: The proposed 0D turbulence model, derived from the 3D approach, has the advantage of describing well the mechanism involving the transfer of kinetic energy in the frequency domain from large to small scales of turbulence due to the shear forces. Results concerning the turbulence submodel will be discussed in Section 4.1.

Mixture Submodel
One of the main issues in Diesel HCCI combustion modelling is that auto-ignition and combustion processes take place in stratified-mixture conditions. In order to define a model able to predict both these phenomena, the modelling of the mixture stratification is of primary importance. To improve this aspect in the model, one solution is to use a presumed PDF approach, commonly used in 3D codes, which describes the mixture stratification using a statistical tool. A PDF is then used to estimate several local variables, such as mixture composition, temperature or reaction rates. Moreover, using the PDF approach implies the addition of only one more state variable and, as a results, a minimum impact on the differential equation system size and consequently on the CPU time cost.

Mixture Fraction Variable
According to the theory developed to study the laminar diffusion flame [48], combining the transport equations of fuel mass fraction, Y F , and oxidizer mass fraction, Y O , it is possible to define a combustion-independent scalar (Schwab-Zeldovitch variable), ϕ, as: where s, the mass stoichiometric coefficient, is defined as: Here, ν i represents the molar coefficient of the stoichiometric combustion reaction and M i is the species molar weight.
The normalization of the variable ϕ with the mass fractions relative to the fuel and oxidizer streams gives the mixture fraction variable, Z, whose values through the diffusive layer vary between 0, in pure oxidizer, and 1, in pure fuel. Z can be expressed as: where Y F 0 is the fuel mass fraction in the fuel feeding stream and Y O 0 is the oxydizer mass fraction in the oxydizer stream. Φ is the equivalence ratio of the nonpremixed flame defined as: It is possible to demonstrate that: where Y F τ is the fuel-tracer mass-fraction, that is the local fuel mass-fraction in the case of a non-reactive mixture.

Presumed PDF
Different kinds of PDF are available in the literature. In order to simulate the mixing process taking place in diffusion flames, the PDF described by a β-function has been judged the best adapted [56]. A β-PDF function is a normalized statistical tool representing the probability, P(X * ), of a generic variable X, defined in its own domain of existence (0 ≤ X ≤ 1).
A β-PDF is completely described by two parameters. The first one is the mean value of the distribution,X, defined as: and the second is the segregation factor of the distribution, δ, defined as: where X 2 is the value of the variance of the distribution computed as: and X 2 max is the maximum value of the variance defined as: The mathematical expression of P(X) is: with: and the coefficients α and β are defined as: Figure 8 shows the typical evolution of the β-PDF. Representation of β − PDF evolution. At the beginning representing a perfectly unmixed mixture (a); at an intermediary stage still presenting two peaks of mixture fraction (b); at an intermediary stage presenting one peak of mixture fraction (c); and at the end of the mixing process a perfectly stirred mixture (d).
The studied configuration is representative of the mixing process taking place in a closed constant volume initially containing pure fuel and oxidizer. The value ofZ is constant in time once the mass proportion of the two gases is fixed.Z can be computed as:Z At t = 0, the oxidizer and fuel are perfectly unmixed (δ = 1). In this case, β-PDF presents two peaks: one representing the pure oxidizer (at Z = 0), and the other representing the pure fuel (at Z = 1). At times t 1 and t 2 , two intermediary configurations are presented. The mixing process is going on so consequently δ decreases; β-PDF is a continuous curve. At t = t end , the mixture is perfectly stirred, and in this case, β-PDF presents one peak centered atZ.

Application of the Submodel to ICE
As shown in the previous paragraph, β-PDF is a tool potentially adapted to reproduce the mixture stratification evolution during a mixing process. The last step in completing the mixture model is to determine the equations giving the temporal evolution of these two parameters defining the PDF (and so of the mixture stratification). For the sake of clarity, it is important to remember that the mixture submodel details the gas mixture composition inside the spray volume. Compared with the previous example, referred to two perfectly-unmixed gases in a closed constant-volume, the ICE application is much more complicated: -both volume and total mass grow in time, -the upper limit of the mixture-fraction existence-domain, Z s , that corresponds to the value Y F s computed by the evaporation submodel, evolves in time as a function of the in-cylinder thermodynamic conditions, -the fuel/ambient gas ratio, and so likewise the mean value of the mixture fraction, varies in time.

Distribution Mean Value Equation
The mean value is computed as the ratio of the total fuelevaporated mass, m F , and the total mass of gas in the spray volume, m S :

Distribution Variance Equation
By definition, the variance of the distribution of Z in the spray is: Differentiating this expression (2) [13]: The use of the β-function to represent the PDF, as already seen, brings several positive aspects, but also imposes some limitations on the mixture-distribution description. For example, Figure 9 shows a situation that cannot be represented by a β-function. To reduce the impact of these limitations on modelling, in the following the upper limit of the distribution domain will be considered as constant and equal to the maximum value reached during the evaporation process, Z s max .
In Equation (54), it is possible to distinguish three contributions. The first term, I, corresponds to the mixing-process contribution (term of dissipation of the variance). It is only thought to depend on turbulence and it is proportional to the turbulence characteristic-frequency via the C diss coefficient.
(2) Complete proof of Equation (54) is given in Appendix. Representation of β − PDF limitation: mixture fraction distribution at the instant t = t 1 (a), mixture fraction distribution at the instant t = t 2 = t 1 +Δt (b). During Δt temperature increases and liquid fuel evaporates at Z = Z s,2 .
The second term, II, corresponds to the entrained ambientgas contribution. The last term, III, corresponds to the fuelevaporation contribution. In this term, it has been done on the assumption that liquid fuel evaporates at saturated thermodynamic conditions in pure ambient gas.

β-PDF Parameters
Once the mean value and the variance of the mass fraction distribution have been computed, they have to be normalized in order to be adapted for use with a presumed PDF. The proposed normalization takes into account the fact that Z s max can be smaller than unity. Mean value and variance are respectively normalized as follows: and: The evaporation of the liquid fuel in the cylinder, especially during combustion, takes place under transcritical or supercritical conditions [57,58]. The evaporation approach proposed in this work is not adapted to describe the evaporation process under these conditions, for the following reasons: -the perfect gas hypothesis is no longer valid, -Clapeyron's law is not valid, -the ambient-gas dissolution in liquid is no longer negligible, -the concept of a liquid-gas interface vanishes.
According to [57,59], using a subcritical evaporation model instead of a tran/super-critical one leads to an underestimation of the value of the maximum fuel mass-fraction. For these reasons the maximum value of fuel mass-fraction has been arbitrarily imposed as equal to unity. This assumption will be partially justified in the following by showing the good agreement between the evolutions of the fuel mixture-fraction distributions computed by the 0D and 3D models, Section 4.1.

Complex Chemistry Approach Motivation
The second part of the Diesel HCCI combustion model, Figure 1, concerns the modelling of the combustion process. Due to the limitations on pollutant emissions, a deeper understanding on the pollutant formation mechanisms must be achieved. The introduction of non-conventional combustion regimes is seen as a possible key point. This is the case for PCCI (Premixed Controlled Compression Ignition), LTC (Low Temperature Combustion) and HCCI combustions, Figure 10 [60]. As shown, the tendency is to obtain a combustion process characterized by lower temperatures and a higher level of local air-excess ratio, in order to get off the region relative to NO x and soot emissions. Nowadays, to reduce the combustion temperatures, it is usual to dilute the fuel-air mixture with EGR because of its high specific heat-capacity and immediate availability. The EGR dilution effect has a direct impact on the chemical kinetics of combustion. In particular, it slows the formation process of active radicals of combustion and reduces their local concentration. It is common, by using high EGR rates, to split the heat release process into two stages: the first, called cold flame, associated to the fuel decomposition in intermediary species; and the second, called main flame, associated to the oxidation process.
To predict all of these scenarios, a complex chemistry approach is required. Figure 11 shows the temperature evolutions relative to two examples of combustion processes Results from SENKIN for EGR rate influence on auto-ignition of a perfectly stirred mixture of n-heptane/air/EGR in a constant volume reactor. Initial states: T 0 = 750 K, p 0 = 25 bar, Φ = 0.6, and X r = 0% and 50% for solid and dashed lines respectively. performed with SENKIN from CHEMKIN [61]. Both configurations are relative to the combustion of a homogeneous fuel/air mixture in an adiabatic constant-volume reactor. Initial conditions do not change between the two configurations, except for the EGR dilution rate. Nevertheless, in the configuration with a higher EGR rate, it is possible to clearly distinguish the two-step auto-ignition process and the reduction of the final combustion temperature.

Tabulation Methodology
Commercial Diesel fuel is not a pure hydrocarbon but a mixture of many different ones; moreover, a univocal definition of Diesel fuel does not exist. Nevertheless the mixture, in order to be commercialized as Diesel fuel, must respect rigorous specifications.
A Diesel fuel-surrogate is often adopted in combustion modelling. It is usually a pure fuel that presents the same auto-ignition properties as the real fuel. For this paper nheptane has been chosen as the Diesel fuel substitute.
In literature several complex chemical kinetics schemes exist for n-heptane. Here, the complete scheme proposed in [62] involving 544 species and 2446 reactions has been chosen. In a 0D model dedicated to a global-system simulation context, in which the CPU computation-cost is a first order evaluation parameter, it is an unreasonable choice to solve the complete chemical scheme. Hence a tabulated chemistry method, inspired by 3D approaches [45], has been set up.
In 1992, [63,64] have demonstrated the existence of attractive trajectories, in the chemical phase space. Based on this, they proposed a theory to simplify the chemical kinetics description, namely the ILDM method. In 2000, [47] proposed to extend this theory by prolonging trajectories with premixed laminar-flame computations: the FPI method. In 2004, [65] replaced premixed laminar flames with auto-ignition flames, and obtained good results on autoignition prediction.
To establish the look-up table, it is necessary to define the main species involved in the chemical scheme. The choice of these species is crucial, because their evolutions have to be sufficient to represent the whole chemical scheme. Here, 8 major species have been identified from mass and energy considerations: C 7 H 16 , N 2 , O 2 , H 2 O, CO 2 , CO, H 2 and H [20].
The thermochemical state of the gaseous mixture of fuel, air and EGR can be completely described with the following four parameters: p 0 : pressure, -T 0 : temperature, -Φ: equivalence ratio, -X r : molar EGR rate in ambient gas. It is defined as: where X O 2 and X F are the oxygen and fuel molarfractions respectively. They constitute an orthogonal base in the phase space, and they have been chosen as inputs of the table. In order to describe the progress of the combustion process, a fifth parameter representing the progress variable, c, is necessary.
The progress variable must satisfy the following properties: c varies in a closed interval between 0 and 1: 0 in the fresh gases and 1 in the burnt gases, c must be monotonic with the progress of combustion, c must be representative of all main reactions.
The generic definition of c used in a FPI approach is: with the mass fraction linear combination, Y lc , expressed as: where t eq stands for the time necessary to reach the chemical equilibrium, σ esp is a coefficient depending on the species and N is the number of species. According to [13], in the model, Y has been defined as: The reaction rate associated with a given species is expressed as: These correspond to an initial condition map of 12960 operating points. The combustion process computation using SENKIN is then performed at each operating point, considering a perfectly stirred constant-volume reactor and adopting the complete chemical kinetics scheme. From each simulation the evolution of the progress variable is determined and discretized in 44 points (fifth table input). Finally, the table is filled up with the reaction rates of each tabulated species. Figure 12 shows the discretization of the temporal progress variable evolution.
Among the eight major species retained, only four have been tabulated: O 2 , CO, CO 2 and H. The others are determined using C, O, H and N atomic balance equations. For example the carbon atomic balance equation is: This tabulation method for n-heptane has been validated and widely discussed [13]. Dependence of the main ignition-delay on the temperature for a perfectly stirred mixture of n-heptane/air in a constant pressure reactor. SENKIN results are compared with experimental data. Initial state: p 0 = 42 bar, Φ = 0.5.
For the sake of clarity, and in order to show the potential of using a complex chemistry approach, some results are reported below. Figure 13 shows the evolution of the main combustion ignition-delay for a given mixture at a given initial pressure as a function of initial temperature. Values computed using SENKIN are plotted against experimental data [66]. The ignition delays are determined from the maximum temperature gradient. It is interesting to note the NTC (Negative Temperature Coefficient) zone in the interval of temperature from 850 K to 950 K. In that interval, an augmentation of the initial temperature corresponds to an augmentation of the ignition delay. Figure 14 shows the evolutions of four major species (C 7 H 16 , O 2 , CO 2 and CO) during the combustion process. Results computed using the tabulated chemistry and the complete chemical scheme (SENKIN) are compared. It concerns an homogeneous lean-mixture with 50% of EGR. It is worth pointing out two aspects: -the species evolution is well predicted if the species is a tabulated one (e.g. CO 2 ). This is not the case for nheptane as it is a reconstructed species that includes all carbonate species other than CO 2 and CO, -potentially, the tabulated-chemistry method allows for the description of any species evolution during combustion. Concerning the study of pollutant emissions, this is a very interesting aspect (e.g. it allows for the determination of thermochemical conditions that favor pollutant formation). Evolution of several species during combustion. Comparison between complete chemical scheme (top) and tabulated chemical scheme (bottom). The difference between reconstructed (tabulation method) and computed (complex chemistry method) evolutions concerning the n-heptane is put in evidence. Curves refer to the combustion process of a perfectly stirred mixture of n-heptane/air/EGR in a constant volume reactor. Initial state: T 0 = 750 K, p 0 = 25 bar, Φ = 0.6, and X r = 50%.

Matching of Spray and Combustion Models
The last step in obtaining the global Diesel HCCI combustion model, Figure 1, consists of matching the spray model with the combustion model.
First of all, it is important to define the relation between the fuel mass-fraction, Z, used to characterize the mixture stratification in β-PDF and the equivalence ratio, φ, representing an input variable of the look-up table: with: where M r and M mix are the molar masses of the EGR and of the mixture respectively, and m a is the mass of air necessary to burn the fuel in stoichiometric conditions. As shown in Equation (63), Z is a function of two look-up table inputs. Therefore particular attention must be paid in interpolating the table values. In the presented ICE application, the initial values of pressure and temperature, p 0 and T 0 respectively, correspond to the in-cylinder pressure and temperature conditions associated with the non-reactive engine-cycle simulation: The most general form of the equation used to compute the mean reaction rate in a stratified mixture associated to a tabulated species is: with: and: where P co is the coupled PDF. Equation (66) is very difficult to solve because of the complex interactions between the different variables.
For that reason, in order to compute the mean reaction rate, several assumptions have to be made: p τ and X r are constant in the domain, -T τ depends only on the fuel mass fraction, Z, c is considered as homogeneous in the domain during combustion: This assumption, certainly the strongest one, probably has an impact on the combustion process. However, it is worth noting that it does not impact the computation of the start of combustion [13]. According to these assumptions, the coupled PDF reduces to a simple single-variable PDF: Now, the mean reaction-rate can be expressed as: The equation defining the local initial temperature of the mixture, T τ , as a function of the mixture fraction is: where Z s and T s are the values of the mixture fraction and the temperature at the liquid-gas interface under evaporating conditions respectively, andT τ is the value of the mean tracer temperature, the definition of which is given below.
The mean progress variable is computed as: The equilibrium data are tabulated as a function of p τ , T τ , Φ and X r , too. To complete the look-up table set of inputs, the specification of the values of p τ and T τ remains. These variables evolve in time because of the thermal losses and the cylinder volume variation. In order to compute p τ , an energy balance-equation, independent of combustion, has been set up: where m cyl is the total mass contained in the cylinder, C v is the constant volume specific heat, Q th represents the thermal loss flux, Q ev represents the thermal flux to evaporate the liquid fuel, h is the enthalpy, n is the mole number and V is the cylinder volume. The subscript τ indicates the variables associated with the thermochemical conditions without combustion (tracer conditions), and the subscript i/o refers to inlet/outlet variables. Enthalpies and specific heats depend on the thermodynamic state of the system: Finally,T τ is computed using the ideal-gas state equation: where ρ mix τ and R are respectively the density of the mixture in tracer conditions and the universal gas constant. Thermal losses have been estimated by using the well established theory proposed in [67]. The reaction rates computed using the look-up table are exact only if the set of inputs refers to a table definition node. For all the other conditions, reaction rates are computed using extrapolation and linear interpolation techniques. It is worth investigating both these aspects. In the model, because of the highly non-linear behavior of the chemical kinetics, all kinds of extrapolation have been avoided; that is, in reading the look-up table, the species reaction rate is set to zero if initial conditions do not belong to the definition domain of the table. Concerning interpolation, equivalence ratio and molar fractions (table coordinates) do not linearly depend on fuel mass fraction (PDF coordinate); consequently an adapted interpolation method is necessary. Interpolation method.
The method, presented in Figure 15, is detailed below. The subscripts/superscripts S , I and 0 are relative to the higher and lower tabulated values (with respect to the reference value), and the reference value respectively. If subscripts, they are associated to the EGR rate, otherwise, if superscripts, they are associated to the equivalence ratio. The interpolation method, used to find the reaction rate associated to Φ 0 and X r 0 , is resumed in the following steps: -determining of the upper and lower tabulated values of X r with respect to X r 0 , -once the value of the variable X r is fixed as X r S /I , using Equation (63), reaction rates ω Φ, X r S /I can be expressed as ω Z S /I , X r S /I , -determining the upper and lower values of Z S /I with respect to Z 0 S /I , -linear interpolation in Z S /I , -to express the reaction rates ω Z 0 S /I , X r S /I as a function of ω Z 0 S /I , Y r S /I , where Y r indicates the EGR mass fraction in the ambient gas, -linear interpolation in Y r , -to express the resulting reaction rate ω Φ 0 , Y r 0 as a function of ω Φ 0 , X r 0 .

MODEL VALIDATION
Validation of the global Diesel HCCI combustion model has been carried out by comparing the computed in-cylinder pressures with engine test-bench experimental data (3) . The engine used is the four cylinder turbo-charged G9T-NADI TM (Narrow Angle Direct Injection). More details concerning the engine can be found in [2].
The model has been validated on a set of eight steadystate operating points that are representative of the engine operating domain, Table 1. The operating point No. 1 has been chosen as the reference operating point. In order to take into account the cylinder-to-cylinder variations, both the minimum and maximum mean cylinder pressures obtained from experimentation are plotted. Two aspects concerning the simulations are noteworthy: -the coefficients of the 0D model are tuned only so as to obtain best results at the reference operating point; they are not modified for the other points, -the initial values of the simulations, such as pressure, temperature and gas mass-fractions, are obtained directly from the engine test-bench. Consequently, their values can not be adjusted to optimize the simulation results. This aspect is an additional challenge for the model. The reference point has been studied in detail in order to have a deeper understanding of the combustion process.

Analysis of the 0D and the 3D Model Results
In order to validate the different aspects correlated to the Diesel HCCI combustion model in an ICE application, such as liquid-fuel evaporation, spray formation, mixing process and chemical kinetics, the reference engine-operating-point, Table 1, was simulated using both the Diesel HCCI 0D combustion model and the TKI-PDF 3D CFD model. The results were then compared, Figures 16-23. (3) The presented experimental data represent the mean value of hundreds of acquired cycle pressure signals.   Figure 16 shows the evolution of the computed incylinder pressures plotted against the experimental pressurecurve envelope. As shown, the simulations are in good agreement with the experiments. Figure 17 compares the evolution of the spray volumes (4) . As shown, near 5 • crankshaft angle, the spraygrowth rate, as obtained from the 3D model, decreases. This is for two major reasons, both depending on geometry: -the first, due to the fact that the spray can hardly penetrate the interstitial volumes (e.g. the volume between the piston squish-area and the cylinder head, when the piston is close to the TDC (Top Dead Center), (4) Concerning the 3D CFD model, spray volume has been defined as the total volume of the cells containing gaseous fuel.

Figure 18
Equivalence ratio field in the cylinder at 5 • crankshaft angle. A cross view on a transverse plane close to the TDC (left), and a cross view on a longitudinal plane containg the injection axis (right). -the second, associated with the toroidal coherent turbulent-structure, which is generated by the impact of the spray on the piston bowl. Figure 18 shows two views of the in-cylinder equivalence ratio field at 5 • crankshaft angle. Figure 19 compares the two total liquid-fuel masses in the cylinder. As shown, the crankshaft-angle intervals necessary to complete liquid-fuel evaporation are very close for both models. The masses of liquid fuel in the cylinder are comparable even though the 0D model underestimates the liquid-fuel evaporation rate. Figure 20 compares the total masses of gases contained in the sprays. The two curves match each other well. The slight variation is mainly due to the geometrical reasons mentioned above. Figure 21 compares the specific turbulent kineticenergies associated to the sprays computed by the 0D and the 3D models. The comparison of the two curves shows  Spray specific turbulent kinetic-energy at the operating point No. 1. Comparison of 0D with 3D modelling results. that the OD model well represents the first-order magnitude variations of the specific turbulent kinetic-energy. Figure 22 shows the evolution of the fuel mixture-fraction PDF inside the sprays for different crankshaft-angle values. The crankshaft angle interval considered extends from the SOI (Start of injection) proximity to the EOC (End Of Combustion). As can be seen all along the interval, the 0D spray mixture model reproduces the air/fuel ratio stratification in the spray well.
Finally, Figure 23 compares the two heat release curves. The two heat-releases are in good agreement, and show that the 0D Diesel HCCI model is able to predict auto-ignition delays (see also Fig. 16) and the chemical reaction kinetics. The sharp shape of the curves shows that the heat release process is more representative of conventional Diesel combustion.

Comparison with Experimental Data
The major challenge for the Diesel HCCI combustion model is to be able to adapt itself to all of the different engine operating-points. An engine operating point is completely defined once initial and boundary conditions are fixed. It follows without saying that model behavior must be sensitive to initial and boundary conditions. Figures 24-26 show the in-cylinder pressures relative to eight different engine operating points. Computed results are plotted against experimental data. Figure 24 shows the model's sensitivity to a variation of the IMEP value (points Nos. 4, 5 and 6). Computed pressure evolutions agree with the experimental curves. The most important aspects in combustion prediction, such as SOC Sensitivity of the 0D model to IMEP variations. Comparison of 0D modelling with experimental data. Curves of pressure refer to the operating points Nos. 4 (top), 5 (middle), 6 (bottom).
(Start Of Combustion), pressure peak value, and pressure peak position, are well reproduced. All three combustion processes are representative of conventional Diesel combustion. Figure 25 shows the model's sensitivity to a variation of the SOI value while the other parameters stay constant (points Nos. 7, 8 and 3). The computed pressures follow the experimental results well. In particular, the main ignition- delay growth and the reduction of the pressure peaks, with a reduction in the SOI, are well reproduced. All three operating points concern high EGR rates that, as already seen, favor cold-flame phenomena. Thus, the combustion processes are representative of HCCI combustion. This aspect is clearly put in evidence by the presence of two peaks in Sensitivity of the 0D model to EGR rate variations. Comparison of 0D modelling with experimental data. Curves of pressure refer to the operating points Nos. 1 (top), 2 (middle), 3 (bottom). the pressure curves: the first associated with the cold-flame heat-release, and the second associated with the main heatrelease.
Finally, Figure 26 shows the model's sensitivity to a variation of the EGR rate value (points No. 1, 2 and 3). This Correlation graphs relative to computed results and experimental data, for the engine operating points listed in Table 1. The circles represent the values referred to the mean experimental curves, cylinder to cyinder variations are pointed out by using error bars. Curves refer to the maximum in-cylinder pressure values (top) and the corresponding crank angle values (bottom). series of operating points considered highlights the potential of the model to predict the combustion process. In fact, with the growth of the EGR rate, the combustion process changes from conventional Diesel to HCCI.
In order to have a more detailed comparison of the 0D Diesel HCCI combustion-model results with experimental data, some characteristic variables of the combustion process are discussed in the following. Figures 27-29 show some correlation graphs between computed results and experimental data (5) . Figure 27 shows the correlation graphs relative to the maximum values of the in-cylinder pressure obtained for the different engine operating points and the corresponding crank angle values. As can be seen there is a good correlation between computations and experimental data.  Correlation graphs relative to computed results and experimental data, for the engine operating points listed in Table 1. The circles represent the values referred to the mean experimental curves, cylinder to cyinder variations are pointed out by using error bars. Curves refer to the crank angle values at which the 5% (top), 10% (middle) and 50% (bottom) of the global injected energy is released. Figure 28 shows the correlation graphs relative to the combustion progress for the different engine operating points. The crank angle values at which the 5%, 10% and Correlation graphs relative to computed results and experimental data, for the engine operating points listed in Table 1. The circles represent the values referred to the mean experimental curves, cylinder to cyinder variations are pointed out by using error bars. Curve refer to the final values of the burned-mass fractions of injected-fuel.
50% of the global injected-energy are released have been pointed out. Correlation graphs referred to crank angle values relative to the end of combustion (e.g. corresponding to the 90% of the global injected-energy) are here omitted as these values, in reason of the very particular adopted engine-control strategies, are not often reached, Figure 29. As shown, the model well predict the mixture auto-ignition angular position (top), which is the most important variable to catch in conventional Diesel and HCCI combustion modelling. The evolution of the combustion process (middle and bottom) is quite well predicted by the model, at least up to the release of half of the injected energy. Figure 29 refers to the final values of the burned-mass fractions of injected-fuel. As can be seen, the combustion processes are often incomplete and the experimental data present wide intervals of variation. These two aspects, if occurring at the same time, are symptomatic of unstable combustion processes. For that reason computed results are considered as satisfactory.

CONCLUSION
This paper presents an original approach to Diesel HCCI combustion modelling. The main contribution of the model concerns the detailed description of the spray using a 0D formalism. In order to carry this out, tools hitherto commonly adopted in 3D CFD models (such as PDF and κ − turbulence model) have been modified and integrated in an innovative 0D model. As a consequence, it is possible to have access to a statistical local description of the air/fuel mixture. It consists of an estimate of the local temperature and local composition values. The spray model has then been coupled to a complex-chemistry tabulation method that computes reaction rates associated with the local thermochemical-conditions. The spray model has been widely discussed and validated with results obtained with the TKI-PDF 3D CFD model developed at the IFP. Finally the global Diesel HCCI combustion model has been validated with experimental data on a wide range of engine operating conditions. This approach is seen to be very promising, mainly for the following reasons: -it is essentially based on physical considerations (this allows a direct comparison of the computed results with experimental data), -it gives a detailed description of the dynamics of spray formation and the combustion process, and can therefore be a useful tool in engine control strategy development, -the detailed chemistry tabulation method gives potential access to all species playing a role in the combustion process; this aspect is very interesting in expanding the study to include pollutant formation mechanisms, -the model is able to predict both conventional Diesel and HCCI combustion, -the model is sensitive to both boundary and initial conditions, -the CPU time required to perform a complete enginecycle simulation is of the order of 10 2 of real time.
Below are some short-term actions envisaged to further improve the existing model: -extending the spray model to the multi-injection strategy management, -investigating the impact on the oxidation process of a combustion progress-variable variance in the reactive mixture, -defining a new complex-chemistry tabulation more orientated towards prediction of pollutant emissions. Main changes shall concern the surrogate fuel choice and the tabulation method. In particular, a species mass-fraction tabulation approach, instead of a species reaction-rate one, seems to be more promising in ICE modelling. Further investigations of these topics will be undertaken by the authors in their future works.

APPENDIX: β − PDF VARIANCE FORMULA PROOF
Mathematical proof of Equation (54) holds as follows. The fuel mass-fraction distribution variance is defined as: (77) hence, expliciting the PDF of the fuel mass-fraction: Differentiating Equation (78): Developping the term (a) of Equation (79): As shown, the term (a) does not contribute to the variance variation.
Concerning the term (b) of Equation (79), it represents the variation of the PDF shape in time. In order to correctly describe the PDF evolution, the different contributions to the variance variation will be determined. More precisely, the impacts of the mixture-zone inputs, namely the entrained ambient-gas and the evaporated fuel, and of turbulence will be presented in detail. Henceforth, global PDF will be subdivided into three contributions: -P A , the probability to have a local mixture-fraction in the spray at Z = Z A = 0; that is, the composition of pure ambient-gas (a perfectly stirred mixture of air and eventually EGR) entrained in the spray, -P s , the probability to have a local mixture-fraction in the spray at Z = Z s ; that is, the mixture composition at liquidgas interface during evaporation process, -P c , the probability to have a local mixture-fraction in the spray different from Z A and Z s . Henceforth, it will be called the complementary mixture-fraction probability.
P A and P s at each instant of time have a well-known composition. Consequently, in the composition space, they can be represented using a Dirac function, ψ.
On the other hand, turbulence acts on the overall mixturefraction distribution and contributes to dissipate the distribution variance. Physically, its action represents the mixture homogenization caused by turbulent-mixing process.
In the following, the different term contributions to the variance variation will be detailed.

Entrained-mass Contribution
At a given time, t 0 , the probability to have in the spray a local mixture-composition at Z = Z A is: after a short interval of time, δt, the spray entrains a certain amount of mass, δm A . Hence: Computing the derivative of P (Z A ) as the limit of the incremental ratio gives: (83) Hence, its contribution to the variance is:

Evaporation-process Contribution
Adopting the same way of proceeding, the evaporationprocess impact on the PDF holds: hence, its contribution to variance is: Fuel evaporation process at liquid-gas interface at saturated conditions implies to introduce in the spray an homogeneous mixture of pure ambient-air and pure fuel in proportions established by the value of Z s . As the necessary ambientair mass is already contained in the spray, its transfer from one value of mixture fraction to another represents a further contribution to the variance variation. In the presented approach, it has been done the hypothesis that liquid fuel evaporates in pure ambient-air. Locally, by definition: where m A is the mass of pure ambient-air in the spray. Hence, for a given fuel evaporation-rate (evaporating at saturated conditions), dm F /dt, the total mass at saturated composition introduced in the spray is: Consequently, the corresponding withdrawn ambient-air mass, dm A w , is: It follows that the impact on the PDF is: and the impact on the variance variation is: as Z A = 0.

Complementary Mixture-fraction Contribution
The used β-PDF is a normalized statistical-tool. Consequently, the following relation must always be verified: Taking advantage of Equation (92), at every time, t, it is possible to know the probability to have in the spray a local composition Z = Z c :