Large-Eddy Simulation of Diesel Spray Combustion with Exhaust Gas Recirculation.

— Simulation aux grandes e´chelles de la combustion d’un spray Diesel pour diffe´rents taux d’EGR — La simulation aux grandes e´chelles (LES Large Eddy Simulations) de la combustion est applique´e l’expe´rience du spray H re´alise´e le cadre de l’ECN ( Engine Combustion Network ). Le ADF-PCM, principe la re´solution de ﬂammelettes approche´es, a e´te´ adapte´ a` la de fac¸on a` simuler la combustion dans cette e´tude. Les effets chimiques complexes peuvent ainsi eˆtre pris en compte. La re´solution de la phase liquide est re´alise´e via une approche me´soscopique Eule´rienne Abstract — Large-Eddy Simulation of Diesel Spray Combustion with Exhaust Gas Recirculation — A Large-Eddy Simulation (LES) study of the transient combustion in the spray H experiment investigated in the frame of the Engine Combustion Network (ECN) is presented. Combustion is modeled using a LES formulation of the ADF-PCM approach, the principle of which is to tabulate approximated diffusion ﬂames based on the ﬂamelet equation to account for complex chemical effects. The liquid phase is resolved with an Eulerian mesoscopic approach coupled with the DITurBC model for the injection. The structure of the combustion resulting from the n-heptane liquid fuel jet is investigated and compared to the literature. A very good reproduction of experimental ﬁndings by the presented LES approach is reported for small EGR rates. Albeit the qualitative effect of increasing the EGR rate is captured, the quantitative quality of the LES predictions deteriorates with increasing EGR rate. One possible explanation for this poor reproduction of EGR effects might be related to the fact that the used semi-detailed scheme was not validated for high EGR rates.


INTRODUCTION
A detailed understanding of Diesel engine combustion is critical to improve engine efficiency while reducing pollutant emissions. To study Diesel engine combustion, constant volume chambers are often used [1,2]. They allow studying transient liquid jet combustion under conditions representative of Diesel engines, without the complexity inherent to such devices.
The simulation of constant volume chambers is still a challenge requiring to model spray formation, liquid evaporation, auto-ignition and flame stabilisation. Three major parameters are used to characterise the Diesel spray combustion: the Heat Release Rate (HRR), the auto-igniting delay and the Lift-Off Length (LOL) [1,3,4].
An  [9] combustion model to predict flame LOL and HRR. Novella et al.
[3] compared different chemical schemes assuming homogeneous combustion at the filter size. The chemistry was assessed on LOL, auto-ignition delay and HRR for different temperatures and oxygen concentrations.
Contrary to RANS, LES takes into account local, instantaneous, spatially filtered flow phenomena, resolving the largest flow scales and modeling only the effects of the smallest ones. It thus appears as having a good potential for addressing unsteady and highly stratified flow found in Diesel combustion. While the literature proposes several studies of LES of Diesel engines [10,11], very few has been published to date on LES of liquid spray combustion in constant volume chambers. Bekdemir et al. [12] used a FGM (Flamelet Generated Manifold) model [13] adapted to LES to predict LOL and auto-ignition delays and compare with experimental findings for Spray H.
In the present paper, the ADF-PCM (Approximated Diffusion Flame -Presumed Conditional Moment) turbulent combustion model initially developed for partially premixed and non-premixed combustion in the RANS context by Michel et al. [14][15][16] is adapted to LES.
It is combined with an Eulerian mesoscopic formalism describing the liquid jet and applied to the LES of the Spray H, experimentally studied by Sandia in the context of ECN [17][18][19]. It proposes a wide range of parametric variations and especially of EGR rates (Exhaust Gas Recirculation).
Section 1 presents the ADF-PCM model and its integration to the AVBP LES solver [20]. Section 2 then describes the numerical set-up and the different cases investigated. Finally, Section 3 compares the obtained results with experimental findings with emphasis on the reproduction of the impact of EGR rate on combustion characteristics.

THE ADF-PCM MODEL
The ADF-PCM model [14][15][16] tabulates auto-igniting approximated diffusion flames. These flames possess the advantage to be computed in a very short time, even for detailed chemical schemes. Indeed, unlike the computation of laminar diffusion flames which requires the resolution of the flamelet equation [21] for each species of the chemical scheme, the ADF approach resolves the flamelet equation only for the progress variable and extracts the chemical source term of this equation from a Homogeneous Reactor (HR) look-up table.

Approximated Diffusion Flame
The ADF approach consists in resolving the flamelet equation [21] for a progress variable Y c : where Y flam c indicates the progress variable transported by the flamelet equation and Z denotes the mixture fraction. The progress variable Y c is defined as proposed by [22]: The scalar dissipation rate v is defined as: with D the mass diffusion coefficient. This quantity measures the rate of mixing and can be expressed, considering a counterflow diffusion flame with constant density, as a function of the mixture fraction Z and the strain rate a [23]: with erf À1 the inverse error function. The source term _ x Y flam c is extracted from a previously built look-up table based on Homogeneous Reactors (HR) at constant pressure. Each HR is initialised for different initial temperatures T u and mixture fractions Z. In the present case, the pressure is set to the experimental value of 42.25 bar. The set of initial conditions is chosen in order to map all the conditions that can be encountered during the simulation. This look-up table stores the species mass fraction Y i and the reaction rate _ x Y c as functions of discrete values of Z, T u and c, the normalized progress variable defined as: where the eq superscript denotes the thermodynamical equilibrium at constant pressure and enthalpy. The reaction rate of the flamelet (Eq. 1) is extracted from the HR look-up table following: with c flam the normalized progress variable reconstructed with Y flam c using Equation (5) and with T u ðZÞ the fresh gas temperature stratification considered as a linear function of the mixture fraction for the approximated diffusion flame. Figure 1 presents a schematic view of the ADF approach.
Equation (1)  The semi-detailed chemical scheme of Seiser et al. [24] containing 1 540 reactions among 160 species is chosen for the HR computation. It has been especially validated on counterflow diffusion flame configurations at 1 bar and on ignition delays up to 100 bar. This mechanism has been used in prior studies of the spray H experiment [4,5] and is recommended by the ECN group.

The ADF-PCM Look-Up Table
In order to take the subgrid scale mixture stratification into account, the species mass fractions and reaction rates of the ADF look-up table are finally integrated using a presumed Probability Density Function (PDF) which is assumed to be a b function. It is parameterised by the local mixture fraction Z, the filtered mixture frac-tionZ and its segregation factor S z . This PDF is denoted as PðZÞ for the sake of clarity. It results in the ADF-PCM look-up table containing all the species transported in the CFD code, as well as the reaction rate of the filtered progress variable: As the temperature stratification T u ðZÞ of the flamelet is considered linear, it is characterised by only two temperatures: the air temperature T a ðZ ¼ 0Þ and the fuel temperature T f ðZ ¼ Z s Þ. Z s denotes the mixture fraction saturation value set equal to 0.5. The mixture fraction segregation factor S z is defined as: where Z v is the mixture fraction variance. Finally, the filtered progress variablec is used instead of the time to Flamelet equation parametrise the look-up table. It is computed using Equation (5) withỸ ADFÀPCM CO andỸ ADFÀPCM CO 2 (Eq. 9). Finally, the look-up table is parameterised as a function ofc,Z, S z ,T a ,T f and a.

LES Transport Equations for the ADF-PCM Approach
The implementation of the ADF-PCM model in the compressible LES flow solver AVBP [20] requires to add some additional transport equations.

Filtered Mixture Fraction
The mixture fraction Z is defined as a fuel tracer, as already performed in several past studies about Diesel combustion modeling [9]. By definition, this tracer is convected and diffused exactly as the species it considers but it is not consumed during the combustion. It varies from 0 for pure air to 1 for pure fuel. The filtered mixture fraction transport equation writes: with D and D t respectively the mass diffusion coefficient and the turbulent mass diffusion coefficient. They are expressed as: where Sc is the Schmidt number equal to 0.75 and Sc t is the turbulent Schmidt number equal to 0.6. The turbulent viscosity m t is modeled via the Smagorinsky model [25] with a coefficient C s = 0.18. The molecular viscosity m is computed using a Sutherland law [25] assuming the viscosity to be independant of the gas composition and close to that of the air. C is the evaporation source term closed as in [26] assuming spherical droplets with uniform temperature. It writes: with n l the filtered droplet density and d the droplet diameter. The filtered Sherwood number Sh and the filtered mass Spalding number B M are also introduced.

Mixture Fraction Variance
The determination of the mixture fraction segregation factor (Eq. 11) requires to transport an equation for the mixture fraction variance Z v : The variance source term due to evaporation _ S Z v is modeled following [27], with a coefficient a _ w set to 0.5: Equation (16) also introduces the filtered scalar dissipation rate computed in the CFD code as: The resolved partṽ RES of the total dissipation rate is directly computed from the filtered mixture fraction fields. In order to close the subgrid scale part, a classical equilibrium hypothesis between production and dissipation [28] of the mixture fraction segregation is retained, leading to a linear relaxation of with D the characteristic size of the filter chosen equal to the cell characteristic size. The strain rate a is then computed in the CFD code from the filtered scalar dissipation rate value:

Filtered Fresh Gases Temperature
In the ADF-PCM model, the approximated diffusion flames are computed for different linear temperature stratifications T u ðZÞ which are characterised by the temperatures of fuel and air boundary conditions. They are determined in the CFD code using transport equations for the unburnt gasesT u at Z ¼Z and the air tempera-tureT a at Z ¼ 0. The filtered fuel temperatureT f is then reconstructed by linear extrapolation as shown in Figure 2. The filtered fresh gas temperatureT u and the filtered fresh air temperatureT a are determined from their respective filtered enthalpy transport equations. The temperatures are deduced from the enthalpies using the mass fractions of the species composing either (fuel or air) stream, for which model transport equations are solved [29].
Concerning the fresh gas temperature, its filtered enthalpyH u is transported as: where _ S H u denotes the evaporation source term modeled as: The filtered enthalpy transfer terms by phase change K and filtered thermal conduction between liquid and gas U write [26]: In these equations, the liquid temperature T l and the sensible enthalpy of the liquid fuel h s;F at the reference temperature T l are introduced. The filtered Nusselt number Nu is also introduced. Finally the Prandtl number Pr is equal to 0.75 and the turbulent Prandtl number Pr t is equal to 0.6. The filtered fresh air enthalpy is transported as: with _ S H a the source term for the evaporation defined as: This transport equation represents a fictive unburnt state formed by air in a mixing between air and liquid fuel droplets. Conductive enthalpy exchanges due to temperature differences between liquid and gas are considered via the U term. The specificity of this equation is its conditioning at Z ¼ 0, where no evaporation occurs. As a result, this equation is similar to Equation (21) except that it does not take into account evaporation. A pre-combustion is used to obtain pressure and temperature conditions prior to the start of injection that are close to those found in Diesel engines. For the simulated case, the pressure at Start of Injection (SOI) is 42.25 bar, the density 14.8 kg.m À3 and the temperature 1 000 K. The liquid n-heptane fuel is injected at 1500 bar and 373 K during 6.8 ms, leading to a total injected mass of 17.8 mg. Details can be found on the ECN web site [17].

NUMERICAL SET-UP
The different initial compositions investigated are presented in Table 1. The presence of CO 2 and H 2 O is due to the pre-combustion phase, necessary to reach Diesel conditions. As the mass fractions of these species are small, they are supposed to have negligible influence over combustion. They are not considered in the simulations and replaced by N 2 . The experiment exhibits four reactive cases with different initial mass fractions of O 2 . Oxygen in the initial gases is gradually replaced by N 2 in order to mimic effects of dilution by EGR [19]. Four ADF-PCM look-up tables are built, each of them representing one EGR rate. The last case with no oxygen corresponds to an inert case dedicated to the determination of liquid and gas penetration.
Scheme of the extrapolation of the fuel temperatureT f from the fresh airT a and fresh gases temperatureT u .
The discretisation of the look-up tables may have a strong influence on the results. In the present simulation, we used the discretisation shown in Table 2. In order to limit the number of points while keeping a good accuracy, the mixture fraction, the progress variable and the strain rate are non-linearly discretised.
The LES of these five cases were achieved using the fully compressible flow solver AVBP for unstructured hybrid meshes [20] using a second order explicit Lax-Wendroff scheme [31]. Subgrid scale turbulence is modeled with the Smagorinsky model [25] using a constant coefficient C s = 0.18.
The description of the gaseous phase is based on spatially filtered Navier Stokes equations. It introduces transport equations for the momentum conservation as well as the species and the energy. The liquid spray is described using a Mesoscopic Eulerian Formalism (MEF) initially developed by Fevrier et al. [32] and adapted for piston engine conditions by Martinez et al. [33]. It requires introducing conservation equations for the droplets density, the volume fraction, the momentum, the sensible enthalpy and the uncorrelated energy resulting from the mesoscopic formalism. As MEF is only valid in regions of small liquid volume fractions, the DITurBC (Downstream Inflow Turbulent Boundary Condition) model is used. It consists in displacing the injection boundary condition downstream from the nozzle exit where the liquid volume fraction is small. Liquid and gas velocities, liquid and gas mass fractions and droplets distribution imposed on the displaced boundary conditions are determined using correlation detailed in [33]. Figure 3 displays a cut plane of the mesh used for these simulations. The mesh is composed of 22.1 million tetrahedral cells with a characteristic size of 60 lm close to the DITurBC inflow plane which gradually increases up to 600 lm at the end of the refined area. The mesh outside the cone is coarse in order to limit the overall mesh size. In this figure, the fuel mixture mass fraction is displayed. The fuel mass fraction goes up to around  Computational mesh in a cut plane through the injector in which the fuel mixture fraction at 2.5 ms is displayed. 0.38, which is still far from Z s chosen equal to 0.5 and below the maximum filtered mixture fraction tabulated equal to 0.42 (Tab. 2). The computation of 1 ms of physical time required 9 h on 120 processors on the CCRT Titane cluster. It reaches 23 h on 400 processors for a reactive case, mainly because of the additional transport equations (for ADF-PCM model and additional species) and the time step limitation due to the combustion.

Spray Formation
The validation of the LES prediction of the spray is performed by comparing liquid and gas penetrations with experimental findings for the non-reactive case (0% of O 2 ). This is achieved by following in time the evolution of the smallest axial distance from the injector outlet at which a specific variable reaches a threshold value. The threshold values used to post-process the LES are those proposed by the ECN group: a liquid volume fraction a l ¼ 0:0015 for liquid penetration and a mixture fraction Z ¼ 0:001 for gas penetration. Figure 4 shows the temporal evolution of liquid and gas penetrations. Numerical and experimental liquid penetrations reach very fast a nearly constant value which is overestimated by the LES compared to the experimental one. Even so, the LES predicted liquid penetration is lower than the minimal LOL value observed which should limit the possible influence of liquid on combustion. The gas phase penetration is accurately reproduced with a small under-estimation after 1.5 ms. Figure 5 presents the evolution of the chamber pressure variation (on the left) and the HRR (on the right). The latter is deduced from the former using integral thermodynamic relations. Before t = 0.380 ms, almost no heat release is visible, as the spray is forming, creating by evaporation and mixing a premixture in which chemical reaction can start. After this induction phase, the autoignition of the so formed pre-mixture eventually leads after t = 0.380 ms to a non-zero combustion heat release, and to an increase of pressure. The HRR curve exhibits a strong peak until t = 0.500 ms, as a result of the burning of the premixture formed in the induction phase. It is consumed fast in what is generally called the premixed phase of Diesel combustion. After t = 0.500 ms the HRR reaches an almost constant value, leading to a linear increase of pressure with time. This phase corresponds to the non-premixed phase of Diesel combustion. This HRR profile is typical of Diesel combustion [18]. The temporal description of the combustion can be spatially observed in Figure 6 showing the temperature fields for different times after SOI. The snapshot at 0.380 ms shows the first auto-ignition spots, characterised by a temperature rise. They are located in regions of lean mixture fractions (Z % 0:05Þ in the downstream part of the jet, which is consistent with experimental observations [18]. The combustion then rapidly propagates and reaches the leading edge of the spray 0.550 ms after the SOI. The snapshots at t = 1.0 ms and 2.5 ms correspond to the non-premixed phase characterised by an almost constant HRR. The flame is anchored at a quasi-fixed axial position called the LOL. It is defined as the minimum axial location for which the temperature rise reaches half of the maximal temperature rise in the domain [12]. This temperature T lift writes:

Reference Case Without EGR
where T init is the temperature of the ambient gas before injection and T max is the maximal temperature in the domain. The time evolution of the LOL is presented in Figure 7. It confirms that the combustion starts downstream in the spray and rapidly propagates upstream towards the injector until it stabilises. The experimentally observed stable LOL is accurately reproduced by the present LES. Temporal evolution of the liquid and gas penetrations for simulation and experiment.

Predicting the Impact of Increasing EGR Rate
Figure 8 compares experimental and simulated pressure rise and HRR against time for the four studied EGR rates. The HRR peak during the premixed phase decreases with increasing EGR rate, as a direct consequence of the slower chemistry resulting from O 2 depletion [1]. Increasing the EGR rate is found to have a negligible effect on the constant HRR level reached in the non-premixed phase, which remains the same for all studied EGR rates. This confirms the view presented in [19] that in the stabilised diffusion combustion phase in which the flame is anchored at the stabilized LOL, lowering the O 2 mass fraction in the initial fresh mixture is compensated by an increased LOL that increases the mixing time prior to combustion. As a result, both effects cancel each other. For all the EGR cases except the highest one, the ADF-PCM model accurately predicts the magnitude of the HRR premixed peak as well as the HRR plateau. Nevertheless the auto-ignition delays and the time occurrence of the HRR premixed peak are increasingly under-estimated with increasing EGR rates. It results in a correct reproduction of the experimental pressure curve for low EGR rates which deteriorates for high EGR rate. The auto-ignition delay, defined as the time to reach a mean pressure rise of 1% of the mean pressure rise at 2.5 ms, is presented in Figure 9 as a function of the EGR rate. The value predicted by the present LES increases with EGR, reproducing the experimentally observed trend. However the LES increasingly underpredicts experimental values as the EGR rate is increased. Temperature fields at four different times after SOI for the configuration with 0% EGR. The grey lines represent equivalence ratio isolines at 0.5, 1 and 2. Temporal evolution of the LOL for the LES of the 0% EGR case. The dashed line represents the experimental stabilised LOL value. Figure 10 compares numerical and experimental LOL as a function of the EGR rate. The ADF-PCM model qualitatively reproduces the effect of the EGR to increase the LOL. Predictions show good agreement for cases without or with moderate EGR rates. For the highest EGR rates however, the LOL is under-estimated.
LES successfully predicts combustion for cases with no or medium EGR rates but its accuracy decreases with increasing EGR, especially for cases with EGR rates superior to 43%. This may be related to the chemical scheme which has not been validated for oxygen mass fraction under 0.17, which in the present case would correspond to 19% of EGR. Another explanation concerns the mesh resolution. With higher EGR, the auto-ignition occurs farther from the injector nozzle, in coarser regions. Finally, the flamelet hypothesis might be blamed as the addition of EGR decreases reaction rates magnitude and therefore impacts the Damkoler number.  Auto-ignition delay as a function of the EGR rate. Evolution of the LOL as a function of the EGR rate.

CONCLUSION
A LES formulation of the ADF-PCM turbulent combustion model is based on auto-igniting strained diffusion flamelets. It was coupled to an Eulerian mesoscopic formalism to be able to study liquid spray combustion. The models were implemented into the AVBP flow solver and applied to the LES of the spray H experiment from Sandia national laboratories [17]. The chemistry of the liquid n-heptane that is injected into a constant volume chamber under Diesel like conditions was tabulated using the Seiser [24] semi-detailed scheme containing 1 540 reactions among 160 species. In the MEF for the liquid spray, the zone close to the injector was modeled using the DITurBC approach [33].
For the reference case with no EGR, the reproduction of experimental findings by the presented LES was very accurate, with a very good reproduction of the experimentally observed mean chamber pressure evolution, HRR, auto-ignition delay and LOL. The impact of increasing the EGR rate was qualitatively reproduced, with an increase of the auto-ignition delay and of the LOL, while preserving an identical HRR once the flame has reached its stable position. However the quantitative prediction deteriorates as EGR rate is increased, with an increasing under-estimation of the auto-ignition delay and LOL with increasing EGR. This is particularly true for the highest studied EGR rate of 62%, for which the LES predicts a too fast combustion compared to the experimental findings. The Seiser chemical scheme might be blamed for this poor reproduction of the impact of increasing EGR rate as it has not been validated for low O 2 concentrations (i.e. high EGR rates). Other possible explanations concern the mesh refinement which is coarser far from the nozzle where high EGR cases autoignites and the validity of the flamelet approach, as EGR decreases the reaction rates magnitude. This will be the subject of future work.

ACKNOWLEDGMENTS
This work was granted access to the HPC resources of CCRT under the allocation 2012-026139 made by GENCI (Grand Equipement National de Calcul Intensif).