Large Eddy Simulation of Spray Auto-ignition Under EGR Conditions

The paper describes the results of a computational study of the auto-ignition of a fuel spray under Exhaust Gas Recirculation (EGR) conditions, a technique used to reduce the production of NOx. Large Eddy Simulation (LES) is performed, and the stochastic field method is used for the solution of the joint sub-grid probability density function (pdf) of the chemical species and energy. The fuel spray is n-heptane, a diesel surrogate and its chemical kinetics are described by a reduced mechanism involving 22 species and 18 reaction steps. The method is applied to a constant volume combustion vessel able to reproduce EGR conditions by the ignition of a hot gas mixture previously introduced into the chamber. Once the prescribed conditions are reached the fuel is then injected. Different EGR conditions in terms of temperature and initial ambient chemical composition are simulated. The results are in good overall agreement with measurements both regarding the ignition delay times and the lift-off heights.


Introduction
Currently there are numerous economic and increasingly widespread ethical concerns regarding the burning of hydrocarbon fuels for energy generation and transport. Since over 90 % of energy generation and transport is fuelled by these means it is likely that this will remain the situation for the foreseeable future. Thus there is a clear need for more efficient and less polluting combustion systems. In the case of gasoline and diesel engines one method of achieving this is through the use of Exhaust Gas Recirculation (EGR), which can give improved combustion efficiency and reduces the emissions of pollutants. This technique is found to be effective in reducing oxides of nitrogen, but has as an observed drawback an increase in the production of particulate matter. Thus there is a clear need for further investigation, both experimentally and numerically. One difficulty in validating numerical simulation methods by comparisons with experimental data that often arises is through the lack of standardisation of experimental combustion facilities. Also the availability of measured data at high temperature and high-pressure conditions is sparse.
The Engine Combustion Network (ECN) [1], provides an experimental set-up capable of reproducing high temperature and pressure conditions with a high level of standardisation. The combustion vessel comprises a cubic constant volume chamber, and allows the measurement of the pressure in the chamber. A Schlieren setup is able to detect the onset of the cool flames, being sensitive to the refractive index variations related to them. However, a robust criterion to determine the exact beginning of the reactions is still missing [2]. The standard definition used by the ECN for ignition delay time sets a threshold at 50 % of the high-temperature chemiluminescence level. The simultaneous use of two pressure transducers in the constant volume vessel is found to be an useful tool for the characterisation of the spray ignition and the pressure rise produced by the combustion process. The lift-off length (LOL) is measured by averaging the maximum light intensity and the dependence upon ambient temperature was investigated. It was noticed that longer ignition delay times corresponds to higher LOL, major deviations in the penetration length were assumed to be caused by differences in the injector orifice diameter.
In the configuration considered in the present paper the fuel, n-heptane in the form of a spray, is representative of practical diesel fuels injected through high-pressure common rail systems. Computation methods capable of providing accurate predictions over a wide range of operating conditions would be a considerable aid to the design of improved combustion systems and it is towards this that the present paper is directed. In fact the lift-off height of a flame determines the temperatures to which the engine's wall may be exposed. Controlling auto ignition delay time allows fuel to be injected more efficiently in terms of timing, which results in lower pollution emissions and more efficient combustion. In addition sooting is strongly related to the efficiency of combustion.
There have been a wide variety of studies on auto ignition of single-phase flows reported in the open literature and a comprehensive review is provided by Mastorakos [3]. However internal combustion engines operates with liquid fuels and unlike single phase flows, very little information under controlled conditions is available. A review of the available experimental data on diesel sprays is provided by Soid et al. [4]. The combustion process is strongly characterised by small scale interactions between turbulence and chemistry. This has led some to undertake Direct Numerical Simulation (DNS) to study most closely and accurately the spray auto-ignition. In this regard Borghesi et al. [5] investigated (using a complex chemical mechanism) the dependence of the ignition delay time on mixture fraction distribution. DNS generally requires a prohibitively large computational effort, so it is not often suited to stimulation under engine-like conditions. The Reynolds Averaged Navier Stokes (RANS) approach has been widely adopted, because of its capability of reproducing the mean fields, often to a tolerable accuracy, at a relatively low computational cost. Numerous studies have been performed, e.g. [6][7][8] using RANS methods. Unlike DNS Large Eddy Simulation (LES) is largely free of Reynolds number dependencies and as a consequence computational costs are manageable, although somewhat larger than those of RANS. A comparison of the capabilities of RANS approaches with LES has been performed by Som et al. [9], emphasising the differences of the two methods while applied under engine-like conditions. LES lies midway between DNS an RANS in terms of computational cost and it provides a much more detailed description of flame structure and turbulent mixing processes that are crucial in diesel spray ignition and combustion. A comparison of different sub-grid-scale models and grid resolution is provided by Xue et al. [10]. The study is performed in the context of Eulerian-Lagrangian framework in diesel-engine like conditions, where it was concluded that LES offered advantages over traditionally adopted RANS models because of the greater details that it provided with computational costs remaining reasonable. Other LES studies include that of Bottone et al. [11] using the LES-Conditional Moment Closure (CMC) model and [44] where a two-phase LES-filtered density weighted pdf was applied. Amongst the findings of these studies is that the ignition delay times display a strong sensitivity to the chemical scheme selected.
The present work uses LES with turbulence-chemistry interactions being accounted for by the joint sub-grid probability density function (sgs-pdf ). Because of the large number of independent variables involved in the pdf equation the stochastic fields method is adopted for its solution. The stochastic field method provides a means of determining the joint sgs-pdf for the all the chemical species concentrations and enthalpy needed to describe chemical reaction. The dispersed phase is described using a Lagrangian particle method, where the spray is represented by an ensemble of stochastic particles. This approach was also successfully used in the simulation of a methanol spray flame [41]. The simulation of two-phase flows in diesel-like conditions presents numerous difficulties such as the modelling of spray vaporisation, breakup, coalescence and spray-wall interaction [12]. Furthermore the lack of data for the physical properties of fuels at high-temperature and high-pressure conditions introduces an important additional source of uncertainty in simulations.
Nevertheless, the oxidation process of the fuel spray and its evolution is crucial to the accurate simulation of the ignition and combustion process; comparative studies of different n-heptane reaction mechanisms are present in [13] indicating the primary importance of the appropriate choice of mechanism. On the other hand the aim of LES is to reduce the computational time to manageable proportions and thus the use of a detailed elementary chemical mechanism is not suitable for this purpose. A trade off between a detailed description of chemistry and computational cost is required.
In the first part of the paper the computational method used for the continuous and dispersed phase is described. A justification for the choice on the reaction mechanism adopted and the description of the method used to characterises the spray phase is provided. Finally ignition delay time and lift-off height for different cases are compared to the experimental data provided by Idicheria and Pickett [14]. Instantaneous temperature profiles through the flame provide an understanding of its evolution although a comparison with measured temperatures is not possible. The result confirms that the LES-pdf method is a valuable tool for the for the description of spray combustion under engine-like conditions.

Mathematical Modelling
Reacting flow problems are described by mass, energy, species conservation and the Navier-Stokes equations. In LES a spatial filter is applied to the equations of motion with the consequence that large scale motions of scale larger than the filter width are described exactly while the contributions of sub-filter scale motions are unknown and must be modelled.

Filtered equations of motion
For a two-phase flow, with the contributions of the dispersed phase being regarded as point sources of mass, momentum and energy, the LES equations can be written: where ρ g is the density of the gaseous mixture, u j is the gas phase velocity,p is the pressure andē ij is the rate of strain tensor. The sub-grid scale stress tensor τ ij = ρ u i u j −ũ iũj is determined via the dynamically calibrated version of the Smagorinsky model proposed by Piomelli and Liu [15]. The filter width is taken as the cubic root of the local grid cell volume. The source terms appearing in the gas phase equations can be evaluated as: (i) where the summation is defined over the number of the droplets present in the cell volume under consideration andṁ (i) is the appropriate source term arising from the i-th droplet. The filtered forms of the conservation equations for specific mole number [kmols/kg] of the chemical species contain the filtered net formation rates of the chemical species through chemical reaction. The direct evaluation of these poses serious difficulties and to overcome these a joint sgs-pdf evolution equation formulation is adopted

Sub-grid joint Pdf
An exact equation describing the evolution of the joint sub-grid (or more strictly the filtered fine grained) pdf, P sgs can be derived by standard methods, e.g. [16] and for the present two-phase flow can be written: where σ is the Prandtl or Schmidt number which are presumed constant and equal, ie unity Lewis number. All of the terms on the left had side of equation (3) appear in closed form an no modelling is required. In contrast the two terms on the right hand side of the equation describe the sub-grid transport and micro-mixing. Since these terms describe processes occurring at the small scales that are not resolved by the LES the terms remain unknown and need to be modelled. In the present work the dynamic Smagorinsky sub-grid viscosity model is used for transport and the Linear Mean Square estimation closure (LMSE) [17] is applied for micro-mixing. Including these models, equation ( where σ sgs is assigned the value 0.7 whereω α (ψ) is, in the case of chemical species the net formation rate through chemical reaction. The termsṁ andṁ α represent the rate of addition of mass and mass of species α to the continuous phase per unit volume through droplet evaporation. The number of scalar quantities, N is equal to the number of chemical species considered plus one (enthalpy). The micro-mixing time scale is obtained from τ sgs

Eulerian stochastic field method
The equation describing the evolution of the pdf, equation (4), is solved using the Eulerian stochastic field method, [18,19].P sgs (ψ) is represented by an ensemble of N s stochastic fields for each of the N scalars namely ξ n α (x, t) for 1 ≤ n ≤ N s , 1 ≤ α ≤ N . In the present work the Itô formulation of the stochastic integral is adopted and the influence of the sgs fluctuations of the dispersed phase is neglected. Thus the stochastic fields evolve according to:ρ where represents the total diffusion coefficient and dW n i represent increments of a Wiener process, different for each field but independent of the spatial location x. This stochastic term has no influence on the first moments (or mean values) of ξ n α . The stochastic fields given by equation (5) form an equivalent stochastic system (both sets have the same one-point pdf, [20]) smooth on the scale of the filter width. The number of stochastic fields, N s employed was eight, whilst preliminary computations were performed using a single field. The influence of an increasingly larger number of stochastic fields would provide a better description of the pdf, although the results to be presented display good convergence at the selected threshold. Further studies focusing on a larger number of stochastic fields may be desirable in future.

PDF modelling of fuel sprays
The dispersed phase is described in terms of a set of macroscopic variables, the droplet radius, r, the droplet velocity v and the droplet temperature θ . In the present case the Weber number may be presumed small, (W e < 20). The equation for the evolution of the corresponding filtered joint pdfP (r, v, θ; x, t), [21] is: where a,Ṙ and˙ are the conditional particle acceleration and the conditional rates of change of the droplet radius and droplet temperature through evaporation and by heat transfer from the surrounding gas phase: they can be written in the general form is the expect value of Dψ k Dt conditioned upon = anywhere in the filter volume. These quantities are unknown and models are required. In order to first model and then solve equation (6) it is replaced with an equivalent system [20] of stochastic ordinary differential equations describing the trajectories of stochastic particles in the phase space {V, R, }. The LES solution provides only the filtered values of the gaseous phase properties at the particle position. The particle are treated as providing point sources and so the influences of unresolved sgs fluctuations of the carrier flow must be modelled. The particle acceleration is represented by a conventional drag-law evaluated in terms of the filtered gas phase properties plus a stochastic Markov model, [22,23] which is used to represent the influence of unresolved sgs velocity fluctuations on particle accelerations and hence dispersion: where V (i) is the velocity of the i th particle, where the subscript p represents carrier gas properties at the particle position, k sgs is the unresolved kinetic energy of the gas phase, C o is a model constant, dW t represents the increment of the Wiener process, g is the gravitational acceleration and τ t is a sub-grid time scale which affects the rate of interaction between the particle and turbulence dynamics, defined as: where the drag coefficient C D is obtained from the Yuen and Chen law, [24]: where Re is the Reynolds number based on the droplet diameter and the relative velocity of the droplet with respect to the gas phase. The sgs kinetic energy is obtained from: k sgs = 2 ν sgs S ij S ij 2/3 , an expression derived using equilibrium arguments.
In the present study, the droplet evaporation is represented using an equilibrium model, [25,26]. In cases where the evaporation rates are low and thus the evaporation is dominated by mass transfer effects, the droplet temperature may be assumed homogeneous and phase equilibrium conditions may be presumed to prevail at the droplet surface. Under these assumptions the rates of change of temperature T (i) and mass m (i) of a single droplet i can be expressed as [27]: where T g is the carrier gas temperature at the i th particle position, h fg is the latent heat of evaporation, Cp g and Cp are the gas and liquid heat capacities at constant pressure and P r g and Sc g are the gas-phase Prandtl and Schmidt numbers. The Nusselt and Sherwood numbers, Nu and Sh are given by the Ranz-Marshall convective correlation, [28,29], together with the modification proposed by Sirignano [30] to account for the effects of Stefan flow on heat and mass transfer. The quantity Sh st is a stochastic Sherwood number contribution representing the influence of sub-grid velocity fluctuations on the evaporation rate, given by: where μ g represents the viscosity of the continuous phase with H M being the convective correction term for mass transfer, [30]. The model constant, C v is assigned a value of unity [21,31]. The particle diameter is then derived from the particle mass. A similar term could be added to equation (9) but was not included in the present work and in the present case the influence of the stochastic Sherwood number is small. The properties of the liquid phase are evaluated using the 1/3 rd rule [26].

Chemical scheme
A detailed description of n-heptane combustion involves a very large number of chemical species and reaction steps. From a computational standpoint an essential aspect of combustion modelling using the sgs-pdf methodology (and most others approaches) is a reduction in the number species for which transport equations have to be solved to manageable proportions. For this reason a reduced mechanism involving 22 species and 18 steps is used. The reduced mechanism is based on a detailed high temperature mechanism containing 1008 steps and 168 chemical species [32], from which a skeletal mechanism with 185 steps and 45 chemical species is derived. By introducing "steady-state assumptions" and leaving unchanged the main reaction paths, [33] the reduced mechanism, Table 1, comprising 18 steps and 22 species is obtained. The majority of simulations involved the reduced mechanism although a small number of computations were conducted with the skeletal mechanism involving 185 steps and 45 chemical species. The cases computed with the 45 species were those with an initial temperature T = 1000K and different oxygen concentrations. For these cases a coarser mesh and one field were used for a preliminary study. The results obtained were closely similar to those resulting from similar computations carried out using the shorter mechanism although the computational times were substantially longer. For this reason the mechanism was discarded for further use. A study on the importance of chemical kinetics by Novella et al. [13] provides evidence, in the context of RANS, of the different performances of different mechanism at low temperature. However computational costs preclude the use of detailed and large reaction mechanisms in LES. It is worth noting at this point that, like the detailed Table 1 Reduced chemical reaction mechanism for n-heptane (18 steps) [33] Step

Results
The case considered is a constant volume combustion vessel, as presented in Figs. 1a and 2, studied experimentally by Idicheria and Pickett, [14,34]. The chamber has a cubic geometry with sides of 108mm. Initially it is filled with a flammable gas-air mixture, with a density, ρ = 14.8km/m 3 , that after ignition, burns to form hot products. After a short cooling period (1-2 s) to allow the mixture to fall to a prescribed temperatures, the spray is injected into the vessel. The chemical composition of the ambient is adjusted by controlling the initial mixture composition. This experimental setup is capable of reproducing EGR conditions and the measurements have been validated by comparison with results obtained in other experimental facilities available in the ECN network. The operating conditions for the vessel are summarised in Table 2.
The environment is at high-pressure around 4 MPa, and the fuel is injected through an injector with a very high pressure ( P inj ) in order to provide a fine atomisation of the droplets: the droplet diameters are of the order of 10μm. Because of the high atomisation pressure, the droplets distribution is approximately homogeneous and the small dimension of the droplets enhances the evaporation process. When, in specific zones, gas-phase equivalence ratios of around unity are reached together with sufficiently high temperatures, the ignition process of the mixture begins. Both from the experiment and the simulation it is (a) Experiment schematic [14].
(b) Computational domain. It is common to identify two parameters to characterise the ignition process: "lift-off height" and "ignition delay time". Lift-off is the distance between the base of the flame and the tip of the injector. In the initial transient part of the ignition, this value may fluctuate, but it eventually stabilises to an equilibrium position. This is reached when the spray jet that tends to convected the flame away from the injector is balanced by the tendency of the flame to propagate towards regions with higher mixture strength, which are located close to the injector.
The auto ignition delay time is the time that elapses between start of the injection and ignition. An accurate prediction of auto ignition delay time is often related to the reliability of the chemical mechanism. The reactions occurring in the "cool phase" (prior to ignition) are rate controlling for this process. In fact the presence of radicals, which are generated in this phase, is responsible for the ignition itself. If the chemical scheme does not include the relevant species or reactions then the auto ignition time may be inaccurately reproduced. Moreover chemical reaction mechanisms are often constrained to be used at specific pressure and temperature ranges and their applicability under EGR conditions is often critical and needs further investigation. In the present study different cases have been simulated corresponding to the different initial burnt gas composition and temperature conditions. These are summarised in Table 3; the oxygen concentration is reduced from atmospheric conditions (21 %) to 10 % and these together with the wide range of temperatures considered provides an important test of the proposed method under engine conditions.

Numerical details
The results to be presented below were obtained using the in-house block-structured, parallel, boundary conforming coordinate LES code, BOFFIN-LES [35]. The code utilises a Fig. 2 Schematic of common rail fuel injector used in the experiment [14] pressure based, low speed variable density formulation and has been applied to an extensive range of flows; further details of method and solution algorithms can be found in, for example, [36,37] and [38]. The stochastic fields equation, equation (5) is discretised using the Euler-Maruyama scheme [39], which is a variant of the commonly used Euler scheme. The Weiner process is represented by time step increments η i √ δt where η i is a [−1, 1] dichotomic random vector. Random numbers η 2i−1 for 1 ≤ i ≤ N/2 are selected from a dichotomic distribution and the remaining numbers are determined from η 2i = −η 2i−1 for 1 ≤ i ≤ N/2. Providing an even number, N of fields is selected the procedure ensures that the mean and variances of the random vector are zero and unity regardless of the number of fields. The stochastic particle equations were also approximated using a similar Euler-Maruyama method. The computational domain, Fig. 1b, which is a cube of sides 108mm, is

Mechanisms comparison
The choice of the chemical mechanism must satisfy multiple constraints. In the present case a low computational cost and accurate representation of the auto ignition delay time and liftoff height is desired. Two different schemes are compared in order to provide a clear idea of the effects that this choice will have on the results. For this reason a well-stirred reactor test case was set-up to evaluate the behaviour of the two mechanisms in isolation. The auto ignition delay times were computed for premixed gaseous n-heptane and air at different equivalence ratios spanning from 0.8 to 1.5. The initial temperatures were selected to cover the same range as used for the LES simulation summarised in Table 3. The comparison between the 22 species and the 44 species mechanisms are shown in Fig. 6. The variation of ignition delay times with equivalence ratio, φ, is shown in Fig. 3b and d, with the former comparing the two mechanisms at different initial temperatures whilst the latter compares them at an initial temperature of T = 1300K. The temperature profiles as a function of time for the two mechanisms, for initial temperatures of 1000K and 1300K and φ = 1.0, are displayed in Fig. 3a and c respectively. The results indicate that the difference between the skeletal and the reduced mechanisms are strongly related to temperature. At the higher initial temperatures the ignition delay times for φ < 1.2 are closely similar with the differences increasing progressively as the mixture strength is increased, Fig. 3d. As the initial temperatures are reduced the differences between the ignition delay times given by the two mechanisms increases progressively, Fig. 3b. The 44 species skeletal mechanism requires very large computational effort when combined with LES and for this reason the 22 species reduced mechanism was selected for the majority of computations. While the reduced mechanism appears to be reliable at high temperatures it has to be admitted that the results for temperatures below 1000K are likely affected by this choice. Thus any discrepancies that arise between simulations and measurements at low temperatures are probably attributable to the choice of kinetic mechanism.

Spray characterisation
Data regarding the droplet distribution at the injector exit is not available so the spray properties at the inlet must be estimated. It is known that the injector used in the experiment is a Bosch solenoid activated pressure injector often used in diesel common rail systems. These injectors operate at high-pressure conditions producing droplets of nearly homogeneous sizes. The spray velocity can be estimated from the n-heptane flow rate, the injector geometry and pressure drop given. To characterise the droplet distribution a Rosin-Rammler expression is used: (12) where Q is the cumulative distribution function for droplet diameter. This needs two parameters, χ and q to be fully specified. The first parameters is related to the Sauter Mean Diameter (SMD) while the second provides a measure of the homogeneity of the distribution. In the initial stage of injection and after interaction with the ambient mixture fuel ligaments break up to form droplets. As the ambient temperature is elevated the droplets evaporate quickly and as a consequence virtually no secondary breakup takes place. The number of 'droplets' in the domain is increased by the injection and reduced by the evaporation and the number of 'droplets' is constant at around 5 × 10 5 . The SMD obtained from the injector specification can be compared to the SMD value computed using various methods [40]. The most relevant to the present configuration is that proposed by Elkbot [42]:  (13) where ρ g and ρ l are the densities of the gas and fuel, ν g is the kinematic viscosity of the gas, σ is the surface tension and where P inj is the injection pressure. Equation (13) is a corrected version of an earlier expression. The SMD obtained from equation (13) is found to be slightly lower than the one given in the technical specification. Simulations performed to investigate the influences in the initial SMD shows that negligible changes result from the two estimates. The evaporation is sufficiently high for the gas phase mixture to rapidly reach its critical value. A small change to the initial droplet SMD, which is already small, has little influence on the on the breakup dynamics and the mixture strength is determined by the mass flow rate rather than by the SMD.
The other parameter to be defined for the spray is the angle of spread and this determines the 'width' of the spray. This can be estimated qualitatively from the available visualisation of the experiment. This estimated value was also compared with the correlation proposed by Siebers [43], a method well suited to high pressure, engine-like conditions such as those studied. The correlation is given by: where the correction factor, which is dependent upon the nozzle geometry is c inj = 0.255. The angle depends on the ambient density in which the droplets are injected and the density of the liquid itself. The resulting computed spray angle is θ spray ≈ 10 • , a value close to that estimated from the visualisation data. Droplets are injected into the vessel at a constant speed and with a varying injection angle. An homogeneous angle distribution is imposed with values varying from zero to θ spray . The parameters used for the Rosin-Rammler distribution, based on the above analysis are: χ = 15μm and q = 4. The cumulative probability distribution function obtained from these parameters is presented in Fig. 4 and a snapshot of the droplet distribution is shown in Fig. 5, where it is possible to see how the droplets are distributed in the vessel and their interaction, depending on their size, with the turbulence of the main flow. The droplet sizes decrease as they travel away from the injector because of evaporation. The evolution of droplet temperature and size is illustrated in Fig. 5a and b. The figures indicate the dynamics of the evaporation: after being injected at low temperature (black), the droplets quickly reach the n-heptane boiling temperature in the external part of the spray and it is evident that the boiling point is reached very close to the injection point. The evaporation rate is

Lift-off height
Typical results illustrating the lift-off height are shown in Fig. 6. Different oxygen concentrations are examined, and corresponding behaviour is observed. However, in all cases the lift-off height stabilises after a transient period. The lift-off is determined as the smallest distance between the injector tip and the nearest cell to it in which combustion occurs, with the same temperature threshold being employed as in the definition of auto-ignition. It is Here the reaction rates are lower because of the lower concentrations of oxygen available, although similar levels of fuel vapour are present. The oxygen is consumed and the combustion process is thus generally slower and occurs at lower temperature. Since the chamber is closed, no fresh air can be entrained in the flame to supply it with oxygen. Using the temperature profiles, the specific lift-off distance from the injector to the region where ignition is first observed can be evaluated. The experimental results are obtained by measuring the chemiluminescence of species such as the OH because of the difficulty of measuring the temperature in the vessel.
The results of varying the oxygen concentration at an initial temperature of 1000K are summarised in Table 4. The calculated lift-off height are in reasonable agreement with the experimental results with the differences, the exception of Case III, being less than 3 %.The  lift-off height has also been simulated for Case I for initial temperatures ranging from 850K to 1300K. The results obtained are summarised in Table 5. As it is evident, there is not the same consistent trend that is observed with varying the oxygen concentration while maintaining the temperature constant. However, the results are in reasonable agreement although the discrepancies are significantly larger for high initial temperatures, although at most they correspond to a distance of 4mm. The comparison with measurements is also influenced by the different techniques used to identify lift-off with chemiluminescence being used in the experiments and temperature in the simulations. Moreover, the location where the flame stabilises is not always sharply defined as the flame location fluctuates and this together with flame 'wrinkling' introduces some uncertainty in the evaluation of lift-off height.

Auto-ignition delay time
In order to calculate the ignition delay, the average temperature across the chamber must be evaluated. The experimental data associates the increase in pressure with the heat released by burning. However, the most accurate way to capture the ignition would be to measure the temperature field directly, although this is clearly not possible in the present configuration. In the simulations, the maximum temperature is used as basis for the determination of ignition delay. The definition is somewhat subjective as a threshold temperature has to be defined. However the very rapid increase in temperature following ignition allows some flexibility on the temperature limit. In this study a threshold of 5 % of the initial ambient temperature is used. Results involving different oxygen concentrations are summarised in Table 6 and Fig. 7a. The simulated ignition delay times show a discrepancy of around 25 % with respect to the experimental data. The observed trends are, however, correctly reproduced and the accuracy improves as the oxygen concentration is reduced. As expected the ignition delay time increases as the oxygen content of the initial mixture is reduced. Table 7 summaries the ignition delay times for different initial ambient temperature ranging from 850K to 1300K. The changes in temperature are examined at constant oxygen  concentration as in Case I. The increase of temperature reduces the ignition delay times as expected, the fuel evaporates much faster thus allowing the ignition probabilities to increase significantly. It is important to note that the simulations generally represent the measured data more closely as the initial temperature is increased. The case with an initial temperature of 1200K is an exception, where the discrepancy is unusually high. The measured value for this case is peculiar because it does not follow the usual correlation between temperature and ignition delay time, being only 0.01ms higher than the T = 1300K case. To help understand better the reason for the differences between the experiments and simulations, ignition delay time against LOL is plotted in Fig. 8. The size of the marker is proportional to the temperature in Fig. 8a and to the oxygen concentration in Fig. 8b. It can be observed that the trends are correctly reproduced for all cases where the temperature is greater than 1000K. In Fig. 8b, where the marker size is proportional to the temperature, the simulations consistently following the measured trend. For the lower temperature (smaller marker) the shift in the horizontal axes, associated with the discrepancies in the auto-ignition delay time, is not correlated with a systematic discrepancy in the lift-off height. Otherwise the simulations are consistent with the experimental data. The discrepancies at low temperature are most probably due to limitations in the reduced reaction mechanism, which, as mentioned earlier, is applicable to high temperature conditions (T 1100K). As far as we are aware there are no alternative reduced mechanisms available covering the range of temperatures encompassed in the experiment. (b) T = 1000 concentration proportional to markers sizes.

Flame propagation
Having considered the two main parameters that provide a quantitative description of auto ignition, it is also worth considering the evolution of the flame through the combustion chamber. This provides information on how the flame is expected to develop and how it interacts with the walls of the enclosing chamber. The instantaneous temperature profiles presented in Fig. 9 show the sequential change in temperature within the domain at different times from injection. The images show that the flame development in the x-z plane through the injector; the injection axis is parallel to the z-axis. The instantaneous temperature plots exhibits the highest temperatures and indicates the flame structure. Initially the cold spray enters the vessel in the form of a jet. The droplets begin to evaporate and consequently reduce the local gas temperature as they absorb heat. Further downstream of the injector the temperatures in the spray begin to rise again until the necessary temperature and mixture strength allows the first auto ignition kernel to appear. Interestingly, the flame initially propagates "backwards" towards the injection point. This is due to the higher values of the equivalence ratio closer to the injection point. However, after a transient period, the flame stabilises at the lift-off height. The stabilising process is shown in Fig. 10, where the lift-off height against time for Case I at an ambient temperature of 1000K is displayed.
The simulated kernel occurs at a larger distance from the injector compared to the experimental measurements. However, the solution rapidly converges to the measured lift-off height after the development of the flame structure. Results for the stabilised location of the lifted flame are in excellent agreement with the measured ones. The stabilisation of the flame at its lift-off height is due to the balance between two different tendencies. On one hand, the spray tends to convect the flame away because of the dynamics and low temperatures of the droplets. On the other hand, the flame propagates towards regions with higher fuel vapour concentrations which lie closer to the injection point. Some of the differences in the dynamics of the stabilisation process of the flame may be associated with the initial  velocity field; in the simulation the chamber was initialised with zero velocity everywhere. This is unlikely to be the case in the experiment for which no initial velocity data is available.

Conclusions
Large Eddy Simulations have been performed on the auto-ignition of an n-heptane spray injected into a constant volume combustion vessel under EGR conditions. The chemical reaction was described by a reduced mechanism involving 22 species and 18 reaction steps. Sub-grid-scale turbulence-chemistry interactions were accounted for by an sgs-pdf transport equation approach, with the equation being solved by the stochastic fields method. The results produced are generally in reasonable agreement with experiment both in terms of liftoff height and auto ignition delay time and their dependence on the O 2 concentration in the initial gaseous mixture. In general the accuracy of the simulations increases with increasing initial temperature, with the largest discrepancies occurring for an initial temperature T = 850K. It was noted that the dependency on the initial temperature of the gaseous mixture was a determining factor in the accuracy of the simulated ignition-delay time and lift-off height. The discrepancies that do occur, however, are most likely to be related to the chemical reaction mechanism used, which is essentially applicable to high temperature regimes, ie T 1100K. As an alternative to the 22 species reduced mechanism a skeletal mechanism (on which the reduced mechanism is based) involving 44 species and 185 reaction was also used in LES simulations but this did not bring any significant improvement. Another source of uncertainty lies with the spray evaporation and dispersion models at high pressures but further measured data is required to investigate this. Although the simulations reproduces the experiments data to a reasonable accuracy the method could be further improved. Future studies may include the atomisation of the droplets by simulating the flow in the injector itself. Modelling effort would therefore be needed to describe the breakup process. In comparison with other LES studies available [44], the present work uses a smaller mesh, and a reduced chemical mechanism with less reaction steps and species which implies a lower computational cost but provides results at least comparable accuracy. The LESpdf method appears to be a reliable tool for reproducing turbulent spray combustion and appears to be applicable to high-pressure EGR conditions.