A PDF method for HCCI combustion modeling: CPU time optimization through a restricted initial distribution

– Probability Density Function (PDF) is often selected to couple chemistry with turbulence for complex reactive ﬂows since complex reactions can be treated without modeling assumptions. This paper describes an investigation into the use of the particles approximation of this transport equation approach applied to Homogeneous Charge Compression Ignition (HCCI) combustion. The model used here is an IEM (Interaction by Exchange with the Mean) model to describe the micromixing. Therefore, the ﬂuid within the combustion chamber is represented by a number of computational particles. Each particle evolves function of the rate of change due to the chemical reaction term and the mixing term. The chemical reaction term is calculated using a reduced mechanism of n-heptane oxidation with 25 species and 25 reactions developed previously. The parametric study with a variation of the number of particles from 50 up to 10 4 has been investigated for three initial distributions. The numerical experiments have shown that the hat distribution is not appropriate and the normal and lognormal distributions give the same trends. As expected when the number of particles increases for homogenous mixture (i.e. high turbulence intensity), the in-cylinder pressure evolution tends towards the homogeneous curve. For both homogeneous and inhomogeneous (i.e. low turbulence intensity) cases, we have found that 200 particles are suﬃcient to model correctly the system, with a CPU time of a few minutes when a restriction of initial distribution is adopted.


Introduction
The necessity to reduce the emissions of CO 2 from automotive engines leads since many years to the generalization of diesel engines, whose one of the advantage compared to spark ignition engines is the lower emission level of this greenhouse effect gas.Unfortunately, diesel engines produce greater NO x emissions and particulate matter (PM) which leads to expensive, difficult and complex exhaust gas after treatment systems.That's why research is heading towards alternative combustion modes which allow a drastic reduction of engine-out emissions.One way to achieve these objectives can be through Homogeneous Charge Compression Ignition (HCCI) combustion.HCCI operation is based on the combustion of a homogeneous mixture of air, fuel and burned gas (internal or external EGR) with reduced combustion temperatures.Furthermore, theoretically the HCCI process eludes locally lean high temperature regions and rich low temperature regions compared to the combustion process for conventional diesel engines, thereby reducing NO x and particulate matter formation.The ignition of this combustion mode is controlled by chemical kinetics [1][2][3][4][5] without flame propagation.The ignition of the mixture occurs in several points in the combustion chamber leading to a sudden combustion.
Experimental and theoretical published investigations [1][2][3][4][5][6][7][8][9][10][11] show that the HCCI combustion process is very rapid, leading to a very rapid rate of heat release due to the near constant volume combustion.The power density is then limited by combustion noise and high peak pressures.Beside this, as observed experimentally by some authors [12,13] the very homogeneous charge is obtained only when the injection of fuel is done far from the engine, whereas for the other usual types of injections (port fuel or direct) inhomogeneities of temperature and of fuel concentration are observed in the combustion chamber.These inhomogeneities can substantially reduce peak pressure, which is an advantage in term of combustion noise.Consequently, the mixture is more often heterogeneous and understanding inhomogeneous autoignition may lead to the development of strategies with several combustion modes to meet the ultra low emissions engine strategies.
To understand this inhomogeneous combustion, modelling tools can provide at low cost more accurate analysis of the in-cylinder processes for the optimization of the combustion chamber design.A direct numerical simulation or 3D modelling of this kind of combustion would require solving a set of coupling equations with a prohibitive computational time.Single zone or multi-zone approaches are not suitable to consider the fluctuations in the zones and the chemical source term is calculated using the mean temperature and gas composition in each zone.One way to reduce computational time is to assume that species and temperature are random variables with a Probability Density Function (PDF) that does not vary in the combustion chamber.PDF approach is often selected to couple chemistry with turbulence for complex reactive flows, since complex reactions can be treated without modeling assumptions, while molecular mixing has to be modeled.Different mixing models can be found in the literature, the most popular is the IEM model (Interaction by Exchange with the Mean) [14,15].This model is also known as the LMSE model (Linear Mean Square Estimation) [14,16].Other models are available in the literature as the MC model (Modified Curl model) [17,18], the EMST model (Euclidian Minimum Spanning Tree) [19] and the PMSR model (Pairwise Mixing Stirred Reactor) [20].This paper describes an investigation into the use of the particles approximation of this transport equation approach applied to HCCI diesel combustion.The goal of this paper is to apply this approach to model Homogeneous Charge Compression Ignition (HCCI) combustion with a minimum number of particles, with a view to generalizing this approach to different diesel combustion modes.This paper describes the complete model, focusing on the optimization of the computation time.The model used here is an IEM (Interaction by Exchange with the Mean) model to describe the micromixing.Therefore, the fluid within the combustion chamber is represented by a number of computational particles.Each particle evolves function of the rate of change due to the chemical reaction term and the mixing term.The chemical reaction term is calculated using a reduced mechanism of n-heptane oxidation with 25 species and 25 elementary chemical reactions developed previously [21].From the mean initial conditions (inlet engine conditions: mass and temperature) obtained from engine test bench, random initial mass and initial temperature of each specie of each particle have been calculated by using three initial distributions: normal, lognormal and hat distributions.The parametric study with a variation of the number of particles has been investigated.To limit the random initial fluctuations, the initial distribution has been restricted through a fixed parameter, estimated from 3D simulation.This approach allows an important reduction of particles number (i.e.CPU time) to reproduce with a good agreements the measured in-cylinder pressure profiles.
The model formulation adopted is described in the next section.Section 3 is divided into three sub-sections.In Sect.3.1, the engine geometry and the operating conditions are specified.Section 3.2 presents the simulation procedure, where the optimization of the CPU time is discussed.The comparisons between the in-cylinder pressure profiles predicted by the model and measured on the engine are presented and discussed in Sect.3.3.The conclusions are reported in Sect. 4.

Model formulation 2.1 General formulation
The Partially Stirred Reactor (PaSR) model has been widely used to study the interaction of turbulence and chemistry.In particle methods this model has provided a sophisticated way for chemical and mixing schemes in the field of combustion [22,23].This model assumes statistical homogeneity of the mixture in the reactor.It accounts for mixing and include large coupled reaction mechanisms.
A stochastic reactor model on the basis of the PaSR is used in this study to simulate the HCCI combustion process.The HCCI engine treated here involves a very early injection of fuel during the cycle.This injection occurs into an air and burned gas (EGR: Exhaust Gas Circulation) mixture.This mixture describes compression and expansion cycle in the combustion chamber (compression, auto-ignition and expansion).The stochastic reactor model is derived from the PDF transport equations assuming statistical homogeneity and can be written as: with the initial condition: where: -F is the mass density function; ψ stands for scalar variables such as chemical mass species (here 25 species) and temperature T (i.e.ψ = ψ(ψ 1 , . . ., ψ S , ψ S+1 ) = (m 1 , . . ., m S , T ); denotes the mean according to F ; -C is a measure for the intensity scalar mixing; the second term of the left side hand side of Equation (1) denotes the effects of chemical reactions, volume change and convective heat losses; the right hand side of Equation ( 1) describes the mixing of scalars by turbulence diffusion.
The mixing term on the right hand side of Equation ( 1) represents the interaction by exchange with the mean (IEM) model.This model is a deterministic model whose main features are simplicity and low computational cost.This model works on the principle that the scalar value approaches the mean value over the entire volume with a characteristic time τ t and assumes statistical homogeneity of the mixture.
The different terms of Equation ( 1) are now described in more detail.The intensity of scalar mixing is introduced through the parameter C in the mixing term given by: where: Cϕ is a constant fixed equal to 2 [15] (some authors [15,24,25] have shown that this constant value can vary from 0 to 2) and τ t is the turbulent time scale.
The term H i in Equation ( 1) describes the change of the MDF due to chemical reactions, change in volume and heat losses.With the assumption that the heat transfer mechanism is due only to forced convection and assuming ideal gas behavior, this term is written as: where: -M i is the molar mass of species i; -M is the molar mass of the mixture; m is the total mass in the cylinder; ρ is the density of the mixture; -A is the heat transfer area function of time; -V is the volume of the combustion chamber function of time; c v is the specific heat capacity at constant volume of the mixture; -R is the universal gas constant; -T w is the wall temperature; -T is the temperature of the mixture function of time; h i is the specific enthalpy of species i; h g is the heat transfer coefficient function of: the mixture temperature, the motored temperature, the volume of the combustion chamber and the in-cylinder pressure.
The molar production rate ωi of the ith species introduced in Equation ( 4) is calculated from: and the rate of progress of jth reaction is given by: where [X] denotes the molar concentration.The numbers k j,f and k j,r are the forward and reverse rate constants for the jth reaction calculated by Arrhenius law, the details of this step with the reduced mechanism used here are given in [21].
The volume of combustion chamber is function of time and its variation versus crank angle degree (θ) is given by: where, V 0 is the clearance volume, B is the cylinder bore, l is the the connecting rod length and a is the ratio of l to the crank radius.The mean density of the mixture is calculated as: and the in-cylinder pressure is given by the ideal gas law: The convective heat transfer coefficient h g is calculated from a modified Woschni correlation [26], which is more suitable to HCCI process than the traditional Woschni correlation [27].

Stochastic process
Equation (1) can be rewritten as: (11) the stochastic process ψ(t) corresponding to the solution F (t) of Equation ( 11) (i.e.Eq. ( 1)) is determined by the system of ordinary differential equations: where, E (ψ i (t)) denotes mathematical expectation, the initial state ψ(0) is distributed according to F 0 (Eq.( 2)).Note that according to Equation (11): In order to create some trajectories of this stochastic process the expectation is approximated as: where, N is the number of stochastic particles.This stochastic particle set approximates the mean density function.Therefore, each particle has a temperature and fluid composition with 25 species, and is thereby related to the fluid particle.The fluid within the combustion chamber is then represented by N computational particles.Each particle evolves according to a set of ordinary differential equations (Eq.( 12)).

Random initial conditions
The IEM model used here is known to have some deficiencies, mainly statistical homogeneity assumption and the parameter N of particle approximation (see Eq. ( 14)).The main advantage of the IEM model is its simplicity.However, as discussed in this paper, the IEM model gives reasonable results with an acceptable CPU time when the mixture in the combustion chamber presents slight heterogeneities of temperature and of concentrations, meaning that the error level of statistical homogeneity assumption is small.On another hand, the number of particles N determines the error level between the approximated solution and the exact solution.The computation time is directly function of this number and its optimization is related to the minimization of particles number.In our first attempt, we have found that a particles number higher than 50 000 is required to reach an accurate in-cylinder pressure profile with high computation time.One way to reduce this factor is to limit the random initial fluctuations through a restricted initial distribution.This approach naturally leads to a reduction of particles number.Additionally, this restriction allows to keep a physical meaning: for example excluding particle temperature equal to 250 K.As presented below, all the tested distributions (hat, normal and lognormal) have been restricted through a parameter ε (with (ε ≺ 1)) representing the initial heterogeneities coefficient.When this parameter is equal to zero the mixture is completely homogeneous.Contrarily, when ε increases the mixture becomes heterogeneous, the level of heterogeneities is directly linked to the value of the initial heterogeneities coefficient.
Based on mean initial conditions (mass and temperature corresponding to the engine operating point), for each stochastic particle, the random initial mass and the random initial temperature are generated with a given distribution.Besides, each particle contains 25 species and each specie has a mass (i.e.m s (0) > 0), generated also by the distribution.Note that at initial conditions, the mixture includes only air and species coming from EGR (Exhaust Gas Recirculation).Therefore, each initial stochastic particle has a temperature and a total mass containing four species with their own mass: N 2 , O 2 , H 2 O and CO 2 .The fuel injection occurs after the beginning of the compression phase, therefore the fuel mass is introduced in the simulation at t > 0. Consequently, at t = 0 the mass of fuel and the mass of the other species included in the reaction mechanism are equal to zero.Note that during fuel injection, the fuel mass is distributed for each stochastic particle according to the distribution.Three initial distributions have been tested in this study: normal, lognormal and hat distributions.As a reminder, for a given distribution noted L and mean initial conditions values noted E χ ≥ 0 (where χ ∈ S ∪ {T }), then for each state variable χ ∈ S ∪ {T }, N random initial conditions are generated.
The hat distribution C (μ, ε) defined by: can be limited by its parameter ε.
In the cases of normal and lognormal distributions, this limitation is directly linked to the variance value of each distribution.The well known normal distribution defined with the parameters μ, σ 2 is equivalent to solving the following system: where the error level is defined by δ − > 0 and δ + > 0 (fixed in this study to: ) and X → N or μ, σ 2 .Using the property that the normal distribution is symmetric, σ 2 is then solution of: where F is the cumulative distribution function of the normal distribution.
The lognormal distribution defined with the parameters μ X , σ 2 X is equivalent to solving the following system: with an error level δ − > 0 and δ + > 0 (fixed as above: F X is the cumulative distribution function of the lognormal distribution.

Experimental description
The experiments have been done with a Renault single cylinder turbocharged engine.The specifications of this engine are given in Table 1.The engine was running with premium diesel fuel in near homogeneous mode, the premium fuel has a cetane number of 55 very close to the cetane number of n-heptane (around 56) used to calculate the reaction term in the model.The amount of EGR was regulated by a valve in the pipe connecting the exhaust and inlet system and was evaluated by measuring the CO 2 concentration in the inlet and exhaust gases.The cylinder pressure was measured with a pressure transducer mounted in the cylinder head.Three injections timing have been investigated: 70, 80 and 90 CAD BTDC.All engine tests have been done on a warmed up engine conditions, the wall temperature required for the simulations has been estimated around 470 K.The engine was equipped with a standard combustion bowl (diameter equal to 5.8 cm).

Simulation procedure
The dimension of the system (12) is high; a splitting technique is then used to decouple the effects of reaction and mixing (more details about the splitting technique are given in [22][23][24]).Note that the term H i in Equation ( 12) has been solved using CVode software package [28], which is very efficient deterministic approach for stiff problems.All the simulations have been done on 13 bi-processor computers, 6 Go Ram, 3 GHz.

Initial distribution effects
The effects of initial distribution presented here have been obtained with a restricted distribution, the particles number set to N = 2000 and to N = 10 4 with different initial heterogeneities coefficient values (but with ε ≤ 20% low heterogeneities).On another hand, the same trends were found for N = 2000 and N = 10 4 ; therefore only results obtained at N = 2000 particles are plotted in Figure 1.In Figure 1, mean temperature history versus CAD is plotted for the three distributions: normal, lognormal and hat with the same turbulence time scale (τ t ) equal to 0.1 s (low turbulence intensity) and the same ε = 10%.This figure shows that the temperature history obtained with the hat distribution is different from the solution obtained with the normal and lognormal distributions.These two distributions give the same solution.The same trends have been observed for several engine operating conditions and model parameter (ε and τ t ).These  differences can be explained by the shape of these three distributions.As expected Figures 2 to 4 show that the hat distribution spreads more around the mean temperature (equal to 350 K).This explains why the hat distribution gives a less stiff solution as observed on Figure 1 in comparison with the other distributions.Note that when the distribution is not restricted a minimum of 50 000 particles number is required to reach the same solution.

Particles number effects: CPU time optimization
In all cases increasing the number of particles, N , decreases the numerical error.However, an increase in N also increases the computational expense in terms of CPU time.For the three tested distributions, the CPU times  was very high (around one week) when a distribution restriction is not adopted, the number of particles required in this case was higher than 50 000.When the distribution is restricted the number of particles was reduced drastically from 50 000 down to 200 as shown in Figure 5.This figure reproduces the temperature profiles obtained at various particles numbers (50, 100, 200, 500 and 2000) with normal distribution (zoom from 6 CAD BTDC to TDC during the main ignition) when the engine speed is equal to 1500 rpm and the initial heterogeneities coefficient equal to 10% (ε = 10%).We note that above 200 particles the same trends are observed and the error level reported in Figure 6 between the solution obtained with 2000 particles, 500 particles and 200 particles is very small, with a stabilized error level above 200 particles for both normal and lognormal distributions.The same trends were found for these two distributions at several engine operating conditions and for all species included in  3.40526e+01 3.10197e+01 2.04166e+02 T the reaction mechanism used in this study (as an example Fig. 7 shows the CO 2 mass evolution with the same conditions of Fig. 1).Besides, for all the cases, the hat distribution requires a higher number of particles to reach the same error level (with a higher computational time) as shown in Figure 1.The variance values reported in Table 2 illustrate that the normal and lognormal distributions have very close variance levels; as expected the hat distribution yields a higher variance.
Therefore, the normal and lognormal distributions are more appropriate than the hat distribution for HCCI combustion.Figure 5 illustrates that a choice of N = 200 stochastic particles is sufficient to achieve an acceptable accuracy in the temperature profiles (i.e.pressure profiles).The CPU time required was about 7 min instead of one week when the initial distribution is not restricted.Amit Bhave et al. [29] have found that a number of stochastic particles equal to 125 is sufficient to achieve an acceptable accuracy in the pressure profiles (i.e.temperature) during HCCI process when the reactive term is not taken into account.These authors have used the same number of particles when this reactive term is introduced.Our results at varying engine conditions point out that 200 stochastic particles are required to reach an acceptable error level, when the initial distribution is restricted.

Turbulent time scale and restriction coefficient effects
In this modeling approach, the turbulent time scale (see Eq. ( 3)) and the initial heterogeneities coefficient ε control the shape of the solution.The parameter ε initiates the initial heterogeneities in the combustion chamber and the turbulent time scale control the dissipation process through the mixing term of the IEM model.In all cases presented in this paper, ε has been estimated from 3D simulations (Star-CD kinetics) [30], using the same reduced reaction mechanism and the same engine operating conditions.However the turbulence time scale (i.e.parameter C of Eq. ( 3)) has been calibrated in order to feet experimental in-cylinder pressure profiles.Figure 8 shows the effects of the parameter ε on temperature histories predicted by the model when the turbulent time scale τ t is set to 0.1 s (low turbulence intensity for this engine operation) during the main ignition at 1500 rpm.These results have been obtained at the same engine conditions as above with normal distribution and with ε = 10% obtained from 3D simulations, this parameter is deduced from fuel concentration gradients amplitude Reference [30].These plots point out at first that the solutions obtained with a particle number set to 200 are very close to the solutions obtained with a particle number set to 2000 for both values of ε (10% and 20%).Secondly, the in-cylinder temperature profile is linked to the value of ε.Besides, as expected when ε increases from 10% (low initial heterogeneities level) to 20% (higher initial heterogeneities level), the in-cylinder temperature evolution is smoother and the main ignition delay is advanced due to the increase of heterogeneities in the mixture.These heterogeneities control the combustion process (i.e.peak pressure and heat release rate) as was shown in [30].
Figure 9 illustrates the effects of the parameter ε at the same operating conditions; we observe that for a given measured pressure the restriction parameter has an important effect on the solution, its value must be estimated accurately.The effects of the turbulent time scale (i.e. the mixing term) when ε is equal to 10% are plotted on Figure 10.This figure reproduces the comparison of the measured and the predicted in-cylinder pressure profiles at the same engine operating conditions for different time scale values.We observe that the time scale has a slight effects on the solution compared to the effects of ε (Fig. 9).On another hand, this figure shows that the in-cylinder pressure evolution obtained with the hat distribution differs from the experimental profile.

Comparisons between model and experiments
Figure 11 presents the in-cylinder pressure profiles predicted by the model and measured on the engine.The engine has been operated with an inlet pressure equal to 2.04 bar, an inlet temperature equal to 380 K, one    2000 particles, predicts a too early ignition (around 5 CAD) and requires a higher particle number (i.e. higher CPU time around 6 h) to fit experimental data.
The same trends as above are observed in Figure 12 with the same injection timing (80 CAD BTDC) and engine speed but with different inlet conditions: inlet pressure equal to 2.56 bar, inlet temperature equal to 400 K, injected fuel mass equal to 19 mg.cycle −1 and EGR rate equal to 64%.For this operating point ε and τ t have been set respectively to 20% and 0.2 s.Here again, the comparison between calculated and measured pressure profiles presents good agreement for the normal distribution.
This comparison has shown that the two parameters ε and τ t are very important to obtain reliable results and that they depend on the engine operating conditions.These parameters need to be modeled in order to make the model autonomous.However, model against experimental data shows overall satisfactory agreement.

Conclusions
Homogeneous Charge Compression Ignition (HCCI) combustion has been modeled using the particles approximation of the PDF transport equation.The model used to describe the micromixing is the well known IEM model (Interaction by Exchange with the Mean).In particle methods this model has provided a sophisticated way for chemical and mixing schemes to account for the inhomogeneities in the combustion chamber using stochastic approach.The random initial mass and initial temperatures of each specie of each particle have been calculated by using three initial distributions: normal, lognormal and hat distributions.The chemical reaction term is calculated using a reduced mechanism of n-heptane oxidation.
The optimization of the CPU time has been done by reducing the number of stochastic particles through a restricted distribution.This step allows reducing drastically the computation time from a week down to around 7 min.In all cases, we have found that a number of particles equal to 200 is a good compromise between the accuracy of the solution and the computation time for normal and lognormal distributions.The hat distribution requires higher stochastic particles with higher CPU time in comparison with the normal and the lognormal distributions.
The in-cylinder pressure predicted by the model was validated against the experimental results.This step has shown overall satisfactory agreement.The differences observed between the measured and calculated pressure profiles presented in this paper can also be linked to possible errors of the in-cylinder pressure measurements due to calibration of the transducer.In this study, we have found that the IEM model gives good results when the mixture presents small heterogeneities.
However, the model needs to be improved by the modeling of the initial heterogeneities coefficient and of the turbulence time scale.As shown in this paper, the solution is strongly τ t linked to the initial conditions; here the parameter ε and the turbulent time scale τ t .Those parameters initiate the heterogeneities in the combustion chamber through ε and the dissipation of these heterogeneities through τ t (i.e.mixing term).In engine applications these parameters are directly function of the operating point and have to be modeled in order to include the physical phenomena they stand for turbulent kinetic energy and its dissipation rate.This step is under progress.
single injection at 80 CAD (BTDC) with an injected fuel mass equal to 16.4 mg.cycle −1 , engine speed equal to 2000 rpm and an EGR rate equal to 54.3%.The simulations were done with the same inlet conditions except the EGR composition, calculated from the model over five consecutive cycles.The particle number was set to 200, ε = 10% and τ t has been calibrated to obtain the best fit with the experimental data, equal to 0.1 s.In this figure, the numerical results obtained with normal and hat distributions are shown (note that the lognormal distribution gives the same trends that the normal distribution).The predicted pressure curve with normal distribution is close to the measured one.The hat distribution including

Table 2 .
Variance of different state variables versus distributions: normal, lognormal and hat.