Complex chemistry DNS of n-heptane spray autoignition at high pressure and intermediate temperature conditions

Combustion and Flame 160 (2013) 1254–1275

Contents lists available at SciVerse ScienceDirect

Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e

Complex chemistry DNS of n-heptane spray autoignition at high pressure and intermediate temperature conditions Giulio Borghesi ⇑, Epaminondas Mastorakos, R. Stewart Cant Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, UK

a r t i c l e

i n f o

Article history: Received 26 November 2012 Received in revised form 5 January 2013 Accepted 8 February 2013 Available online 12 March 2013 Keywords: Diesel engines Autoignition Complex chemistry Spray DNS

a b s t r a c t Direct Numerical Simulations (DNS) of turbulent n-heptane sprays autoigniting at high pressure (P = 24 bar) and intermediate air temperature (Tair = 1000 K) have been performed to investigate the physical mechanisms present under conditions where low-temperature chemistry is expected to be important. The initial turbulence in the carrier gas, the global equivalence ratio in the spray region, and the initial droplet size distribution of the spray were varied. Results show that spray ignition exhibits a spotty nature, with several kernels developing independently in those regions where the mixture fraction is close to its most reactive value nMR (as determined from homogeneous reactor calculations) and the scalar dissipation rate is low. Turbulence reduces the ignition delay time as it promotes mixing between air and the fuel vapor, eventually resulting in lower values of scalar dissipation. High values of the global equivalence ratio are responsible for a larger number of ignition kernels, due to the higher probability of finding regions where n = nMR. Spray polydispersity results in the occurrence of ignition over a wider range of mixture fraction values. This is a consequence of the inhomogeneities in the mixing field that characterize these sprays, where poorly mixed rich spots are seen to alternate with leaner ones which are well-mixed. The DNS simulations presented in this work have also been used to assess the applicability of the Conditional Moment Closure (CMC) method to the simulation of spray combustion. CMC is found to be a valid method for capturing spray autoignition, although care should be taken in the modelling of the unclosed terms appearing in the CMC equations. Ó 2013 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction HCCI combustion has been attracting significant attention in recent years due to its potential for reduced NOx and particulate matter emissions whilst maintaining the high thermal efficiency typical of compression–ignition engines [1,2]. At the low temperature and high pressure conditions corresponding to HCCI, combustion occurs through volumetric autoignition. Controlling the ignition delay time of the system becomes thus fundamental in order to achieve the desired levels of pollutant emissions and to avoid excessive rates of heat release, which would result in engine knock. Achieving this goal requires a profound understanding of the physics of spray autoignition. Improving our knowledge of ignition phenomena in two-phase flows is also critical for the development and validation of advanced turbulent combustion models for spray combustion. Two popular models for simulating combustion in single-phase flows are the Flamelet model of Peters [3] and the Conditional Moment ⇑ Corresponding author. Current address: Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, MS 125-130, Pasadena, CA 911098099, USA. E-mail address: [email protected] (G. Borghesi).

Closure (CMC) of Bilger and Klimenko [4]. In the last decade, formulations of both these models have been derived to study combustion in two-phase flows [5–7]. Preliminary results have been encouraging, but several aspects need to be clarified; these include the choice of the spray flamelet library, or the modelling of terms describing droplet evaporation and scalar mixing in the CMC equations. Addressing these points requires a profound understanding of the interactions existing between evaporation, mixture formation and chemical reactions in sprays. Previous numerical work [8] on the autoignition of single-phase flows revealed that ignition occurs first in those regions of the flow where the mixture fraction has a well-defined value, which was named most reactive mixture fraction nMR, and the scalar dissipation rate is low. These results, obtained with two-dimensional DNS with simplified methane chemistry, were later confirmed by additional DNS results, including both three-dimensional domains [9,10] and more complex kinetic schemes [11]. Valuable information on the structure of the ignition kernels in single-phase flows has also been obtained from experiments; in particular, it was found that formaldehyde plays an important role in methane/air jet flames, acting as an ignition precursor [12]. Similar results were obtained for hydrogen/nitrogen jets [13]; in this case, however, the role of the ignition precursor was played by the HO2 radical.

0010-2180/$ - see front matter Ó 2013 The Combustion Institute. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.combustflame.2013.02.009

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

Despite the large body of studies available on the autoignition of single-phase flows, which have also been the subject of a recent review [14], considerably less information is available for spray autoignition. O’Loughlin and Masri [15] observed that, similarly to the purely gaseous case, formaldehyde appears early in the autoignition of turbulent dilute methanol sprays. A marked difference with single-phase flows autoignition was the frequent appearance of double reaction zones. Wang and Rutland [16,17] investigated numerically with 2D-DNS the ignition of n-heptane sprays at high temperature and moderate pressure levels. They found that an increase in the bulk gas temperature resulted in a faster ignition event, while an opposite behavior was observed for the global equivalence ratio. Ignition occurred at a few isolated locations at first, and then spread out rapidly to the whole mixture. Schroll et al. [18] used three-dimensional DNS with one-step global chemistry to analyze the behavior of spray autoignition in mixture fraction space. Their study revealed that the first spots to ignite occur in those regions where n ’ nMR and the scalar dissipation rate is small. These results are fully consistent with the observations made for autoignition in single-phase flows [8]. A close relationship between the initial droplet diameter and mixing was also noted. However, due to the use of a single-step chemistry, nothing could be inferred on the ignition behavior of sprays in the range of operating conditions relevant to HCCI and diesel engines. This work deals with the numerical simulation of autoignition in n-heptane sprays using three-dimensional DNS coupled with complex chemistry. The use of detailed chemistry in DNS has been hindered so far by its enormous computational requirements. Three-dimensional DNS simulations have often relied on global chemistry to describe the complex phenomena leading to ignition [10,18], while the use of reduced and skeletal kinetics has been confined to two-dimensional configurations [16,17,19]. The influence of the kinetic scheme on the numerical simulation of selfigniting n-heptane jets was recently investigated by Viggiano [20]. Her study revealed that, depending on the operating conditions investigated, there may be situations in which ignition is fluid-dynamics controlled, and others in which chemical kinetics plays a dominant role. In the latter case, the use of a detailed scheme is required to predict ignition properly, while in the former even global chemistry may be sufficient to obtain reliable results. The suitability of a two-dimensional configuration for studying autoignition in non-homogeneous mixtures was discussed by Sreedhara and Lakshmisha [10]. Their study revealed that the structure of the autoignition sites in non-premixed flows remains the same irrespective of whether the flow is represented in two or three dimensions. These results were obtained using a four-step global n-heptane mechanism, and their validity for more complex chemistry has not been assessed yet. Differences may be expected, especially in those situations where the interaction between chemistry and turbulence is strong and a correct representation of the turbulent eddies is mandatory. The above discussion justifies the choice of a three-dimensional configuration with detailed chemistry treatment for the simulations presented here. It should be emphasized that, to the authors’ knowledge, no 3D-DNS with detailed chemistry of spray autoignition has been performed in the past; previous work on two-phase flows with complex chemistry DNS was limited to the spark-assisted ignition of sprays [21,22]. The scope of this work is to provide some physical insight into the complex interactions between turbulence, chemistry and fuel/ air mixing that characterize hydrocarbon sprays autoigniting at low temperatures and moderate pressure levels, for which no fundamental computational study exists. More specifically, our aim is to verify whether a most reactive mixture fraction can still be found for sprays autoigniting at intermediate temperature conditions, and if the negative correlation between heat release and

1255

scalar dissipation rate that characterizes gaseous flow autoignition continues to hold for these conditions. We also want to investigate how changes in the operating conditions, such as the initial turbulence strength in the bulk phase, the global equivalence ratio in the droplet-laden region and the initial droplet size distribution of the spray affects mixture formation and the ignition process. The data obtained from the simulations will also be used to assess the predictive capabilities of the CMC method for spray autoignition. 2. Mathematical formulation 2.1. Liquid phase The liquid phase is treated using a Lagrangian point-source formulation. The temperature within each droplet is assumed uniform, and droplet–droplet interactions are neglected. The motion of the droplets is subjected to the aerodynamic drag force only. A thin film assumption (subscript f) is used to calculate the droplet evaporation rate and the amount of heat exchanged with the gaseous phase. Physical properties in the thin film are computed according to the 1/3 rule [23]. For each droplet d, its position xd, velocity vd, diameter ad and temperature Td are computed by solving the following system of equations [23]:

dxd ¼ vd dt dv d uðxd ; tÞ  v d ¼ dt svd

ð1Þ ð2Þ

2

a2 dad ¼  dp dt sd "  0:38 # dT d 1 Lv T crit  T d ¼ T Tðxd ; tÞ  T d  BT;d dt W F cp;f T crit  T ref sd

ð3Þ ð4Þ

u and T are the gaseous velocity and temperature, evaluated at the droplet location. Tref is the liquid boiling temperature at the reference pressure pref, Tcrit the liquid critical temperature, Lv the molar latent heat of evaporation at the reference temperature and WF the molecular weight of the fuel species. Note that the value of the heat of vaporization at the droplet temperature is evaluated using a power-law originally proposed by Watson [24]. This is not specific to n-heptane, and is also valid for other fuels. The Spalding number for mass transfer is given by:

BM;d ¼

Y sF;d  Y F ðxd ; tÞ 1  Y sF;d

ð5Þ

where Y sF;d is the fuel vapor mass fraction at the droplet surface, evaluated according to the Clausius–Clapeyron law, and YF(xd, t) is the fuel mass fraction at the computational cell where the droplet resides. The Spalding number for heat transfer BT,d is evaluated from BM,d following the procedure described by Abramzon and Sirignano [23]: cp;F Shc

BT;d ¼ ð1 þ BM;d Þcp;f Nuc LeF  1

ð6Þ

LeF is the Lewis number of the fuel species, while Shc and Nuc are the modified Sherwood and Nusselt numbers. They are evaluated here following [23]:

Sh0  2 ; Sh0 ¼ 1 þ ð1 þ Red PrÞ1=3 f ðRed Þ FðBM;d Þ Nu0  2 Nuc ¼ 2 þ ; Nu0 ¼ 1 þ ð1 þ Red ScF Þ1=3 f ðRed Þ FðBT;d Þ

Shc ¼ 2 þ

ð7Þ

where Red is the droplet Reynolds number, ScF the Schmidt number of the fuel species, and Pr the Prandtl number, set equal to 0.7.

1256

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

Additionally, f ðRed Þ ¼ minð1; Re0:077 Þ, and F(B) = (1 + B)0.7ln (1 + B)/B. d The relaxation times appearing in the droplet governing equations are given by [23]:

svd ¼

qL a2d



2 3

18lf 1 þ 16 Red



ð8Þ

qL a2d PrLef 1 4Shc lf lnð1 þ BM;d Þ q a2 PrLef BT;d cp;L sTd ¼ L d 6Shc lf lnð1 þ BM;d Þ cp;f

spd ¼

ð9Þ ð10Þ

where qL the liquid fuel density, and cp,L the heat capacity of the liquid. The fuel considered in this work is n-heptane (pref = 1 bar, Tres = 371.58 K, Tcrit = 540.15 K, Lv = 31.80 kJ/mol and cp,L = 2260.40 J/(kg K)).

! Ns X @Y a @Y b Ya V a;j Y a ¼ Da þ Db @xj @xj b¼1

The thermal conductivity is computed using the expression P s k = 2.58  105  cP(T/298)0.7 W/(m K), with cP ¼ Na¼1 Y a cP;a . The molecular viscosity is then given by l = Pr(k/cP). The diffusion coefficient of species a is obtained from Da = k/(qcPLea), where Lea is the Lewis number of the ath species, here assumed constant and equal to one. The liquid source terms appearing in Eqs. (11)–(14) are:

Cm ¼  Cui ¼ 

1 X dmd ad V d dt X 1 dmd v d;i V

ad

d

CE ¼ 

The Navier–Stokes equations, supplemented by Ns  1 transport equations for the species mass fractions (Ns being the number of species appearing in the chemical mechanism used), are solved to describe the temporal and spatial evolution of the gas phase:

CY a ¼ da;F Cm

ð11Þ ð12Þ

P ¼ qRT E¼

ð15Þ

Ns X p 1 Y a ha  þ um um q 2 a¼1

ð16Þ

In the above equations, R is the universal gas constant, Wa the molar mass of species a and ha the corresponding species enthalpy, defined as:

ha ðTÞ ¼

Z

T

T0

ad ðxÞ ¼

  A  exp  12 jx  xd j2 ; if jx  xd j < 2ad;0

ð25Þ

otherwise R Vd

  exp  12 jx  xd j2 dV, and Vd is the spherical vol-

ume centered at the droplet location and having radius equal to 2ad,0, where ad,0 is the initial Sauter mean diameter of the spray. The function ad redistributes the droplet source terms in the neighborhood of droplet d according to a distance function that decreases exponentially with increasing distance from the droplet center. The distance function has a cut-off length of 2ad,0 to avoid any influence of the droplet source terms in the droplet far field. This approach for assigning the liquid source terms to the gas phase is not new, and similar formulations have already been experimented successfully by other authors [17,25]. It is chosen here to improve numerical stability and to ensure that the evaporation time of the droplets is independent of the grid size, as the droplet source terms are now redistributed over a spatial region whose size only depends on ad,0, and not on the cell volume. The chosen cut-off length is the result of some numerical experiments. It ensures that the temporal evolution of a2d for a droplet at P = 0.1 bar, T = 741 K and ad,0 = 0.7 mm, for which experimental data are available [26], is correctly reproduced.

ð17Þ

where cP,a is the specific heat capacity of species a, and ha is the corresponding species enthalpy at the reference temperature T0. The fluid studied here is assumed to be Newtonian, and thus the viscous stress tensor is given by:

  @ui @uj 2 @um þ  di;j @xj @xi 3 @xm

ð18Þ

The heat flux vector is:

qj ¼ k

ð23Þ

2.3. Numerical procedure and computational parameters 0

cP;a dT þ ha 0

sij ¼ l

(

where A ¼ V=

Ns X Ya Wa a¼1

!

where md ¼ qL pa3d =6 is the mass of droplet d, hF(Td) is the fuel vapor enthalpy at the droplet surface, daF is 1 for the fuel species, 0 otherwise, and V is the volume of the Cartesian grid cell corresponding to the droplet location. The non-dimensional function ad is evaluated for each grid node according to the following expression:

ð13Þ

where sij is the viscous stress tensor and qj the heat flux vector. Va,j _ a are the diffusion velocity and the chemical reaction rate of and x species a. The last term in each equation is the liquid source term, detailed later. Note that the Einstein summation convention for repeated Roman indices has been used. This convention will be adopted in the remainder of this work, except where differently specified. Eqs. (11)–(14) are closed with the following state equations:

ð22Þ

dt

ð24Þ

0;

ð14Þ

ð21Þ

1X Tðx ; tÞ  T d dmd 1 dmd v 2d;i ad cLP md d T þ hF ðT d Þ þ V d 2 dt dt sd

2.2. Gas phase

@ q @ quj ¼ Cm þ @t @xj @ qui @ qui uj @P @ sij ¼ þ þ q F i þ Cu i þ @xi @xj @t @xj @ qE @ quj E @Puj @ sij ui @qj þ ¼ þ  þ CE @t @xj @xj @xj @xj @ qY a @ quj Y a @ qV a;j Y a _ a þ CY a ¼ þ qx þ @t @xj @xj

ð20Þ

Ns @T X þ qV a;j Y a ha @xj a¼1

ð19Þ

The diffusion velocity of species a, Va,j, is modelled according to a modified Fick’s law to ensure mass conservation:

The equations governing the evolution of the gaseous and liquid phases were solved numerically using the DNS code SENGA2 [21,27]. First and second order spatial derivatives were evaluated using a 10th-order explicit central difference scheme [28]. Periodic boundary conditions were applied in all spatial directions. Operator splitting between transport in physical space and chemistry was implemented to ease the constraints on the time-step size. DNS with operator splitting was studied in detail and validated in [29]. A fourth-order, low storage, explicit Runge–Kutta scheme [30] was used to advance in time the gas-phase transport equations without chemical source terms. Then, chemical reactions were computed using the implicit solver VODPK [31]. The timestep size Dt was chosen to be 5  109 s: this is more than 10 times smaller than the smallest non-chemical time scale in the simulation, that is the acoustic time-scale in a grid cell. The droplet

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

transport equations were advanced at the beginning of every time step using the same explicit Runge–Kutta scheme employed for the gas-phase. The present code has been previously used for DNS of spark ignition in sprays [21,22]. The cubic domain had a length of either L = 2.1 mm or L = 2.9 mm, depending on the initial turbulence level of the bulk phase. The grid resolution was equal to Dx = 22.8 lm, resulting in 92 and 128 computational nodes respectively along each spatial direction. The initial ambient pressure was P0 = 24 bar, and the initial ambient gas temperature was Tair = 1000 K. It should be noted that the chosen pressure value is well below the levels found in real diesel engines. However, it is also below the critical pressure of n-heptane. Simulations at higher pressure levels would require a detailed thermodynamic treatment of the species and energy fluxes, including Soret and Dufour effects, and the inclusion of off-diagonal diffusion coefficients in the diffusion matrix [32]. The turbulent velocity field was initialized according to a Batchelor–Townsend energy spectrum with integral length scale L11 = L/ 6. Turbulence was isotropically decaying, and there was no initial mean flow. Two values of the initial turbulent velocity fluctuations were considered, u00 ¼ 0:25 m=s and u00 ¼ 0:7 m=s. The corresponding values of the Kolmogorov length scales were gK = 39.7 lm and gK = 27.9 lm. Fuel droplets were initially distributed in a layer extending from x = 0.35L to x = 0.65L. A sketch of the computational domain grid is shown in Fig. 1. The initial droplet velocity was equal to that of the background carrier-phase. The overall equivalence ratio in the droplet-laden region was either U0 = 2 or U0 = 4 for all simulations. The group combustion number G0 was evaluated as [33]:

1257

was confirmed a posteriori, since, for all cases investigated, no collision event was detected. The initial fuel temperature and the initial Sauter mean diameter of the spray were equal to Tf = 450 K and aSMD,0 = 25 lm for all cases investigated. All but one of the sprays studied were monodisperse. For the polydisperse case, the droplet diameters were initialized according to a Rosin–Rammler distribution [34]:

  q  ad;0 Fðad;0 Þ ¼ 1  exp  r

ð27Þ

where nT is the total number of droplets in the droplet-laden region, and Red,0, ld,0 are the initial droplet Reynolds number, equal to zero here, and the initial mean inter-droplet distance. The sprays considered in this work are dilute, as the initial volume occupied by the droplets is approximately 1.5% (3% for the case with U0 = 4) of the droplet-laden layer volume. The initial mean inter-droplet distance is approximately 10 droplet diameters for the case with U0 = 4, and larger for the remaining cases. Given the nature of the sprays investigated, no algorithm for handling droplet–droplet collisions was implemented in the computational code. The validity of this choice

where F(ad) is the cumulative distribution function of the droplet diameter and represents the relative amount of liquid contained within droplets having diameter smaller than ad. q is the shape parameter and controls the width of the distribution. It was equal to 3.5 in our study. The value of the amplitude parameter r was adjusted to obtain aSMD,0 = 25 lm. Since the Rosin–Rammler distribution assumes an infinite range of droplet sizes, the tails of the distribution were not considered in the spray initialization. More specifically, only values of ad that satisfied the condition 0.05 < F(ad) < 0.95 were used to assign the droplet diameters. The resulting droplet size distribution of the polydisperse spray is shown in Fig. 2. It should be noted that, although the most probable droplet diameter fall around 20 lm, a large part of the liquid mass is associated with larger drops. There is a considerable presence of small droplets in the spray (e.g., ad,0 < 20 lm) but, due to their corresponding small size, their contribution to the total liquid mass is negligible. Details of the various simulations performed can be found in Table 1, together with the corresponding reference times for eddy turnover st, droplet evaporation sev and homogeneous mixture ignition sid,0. The reference evaporation time was defined as the time needed for a motionless droplet with initial diameter ad,0 to evaporate in a hot, quiescent medium. It was evaluated by solving numerically Eqs. (3) and (4) with constant conditions in the gas phase, e.g., T = T0 and Yf = 0. In the case of a polydisperse spray, sid,0 was taken to be the reference evaporation time of the largest drop initialized in the computational domain. The reference ignition delay time will be defined in Section 3.1. The operating conditions for the simulations presented in this work were chosen so as to obtain comparable values for at least two of the three reference times, and to investigate the response of spray autoignition to changes in turbulence strength, initial

Fig. 1. Sketch of the computational domain. Grey corresponds to initial droplet laden region.

Fig. 2. Initial droplet size distribution and cumulative liquid volume fraction for Case D in Table 1.

  2=3 ad;0 G0 ¼ 1:5LeF 1 þ 0:276Sc1=3 Re1=2 d;0 nT ld;0

ð26Þ

1258

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

Table 1 Turbulence and spray properties for the numerical simulations presented in this work. For all cases investigated, Tair = 1000 K, Tf = 450 K and lt = L/6, where lt is the turbulence integral length scale. Ret is the turbulent Reynolds number, Ret = u0 lt/m. Case

L (mm)

Ret ()

G0 ()

a0 (lm)

q

st (ms)

sev (ms)

U0

A B C D

2.1 2.9 2.9 2.9

19 73 73 73

15.2 15.2 30.2 19.6

25 25 25 25

– – – 3.5

1.52 0.75 0.75 0.75

0.82 0.82 0.82 1.60

2 2 4 2

mass loading, oxygen dilution and initial droplet size distribution. In particular, a comparison between Cases A and B reveals the effect of u0 ; between B and C the effect of the overall equivalence ratio (e.g., initial number of droplets in the spray); and between B and D the effect of polydispersity. In the data processing, a mixture fraction n based on elemental mass fraction was used. For a generic fuel species CmHn, this is given by:

n¼ b¼

b  bO bF  bO ( Ns X

ð28Þ )

Ya Wa

ð29Þ

1 ðm þ n=4ÞW O

ð30Þ

cC W C aC;a þ cH W H aH;a þ cO W O aO;a

a¼1

cC ¼

1 ; mW C

cH ¼

1 ; nW H

cC ¼ 

where al,a denotes the number of atoms of element l in species a, and Wl is the molar mass of element l. bO and bF correspond to the value of b in the oxidizer and fuel streams respectively. 2.4. Chemical mechanism Reaction rates were computed using the reduced n-heptane kinetic mechanism of Liu et al. [35]. The mechanism was derived from a skeletal one consisting of 43 species and 185 reactions. It contains 22 nonsteady-state species, reacting according to 18 global steps. The mechanism was validated in [35] in terms of ignition delays for homogeneous mixtures at different equivalence ratios, temperatures and pressures by mean of experimental data from [36]. It has also been used for non-premixed autoignition problems by Wright et al. [37] and Borghesi et al. [7]. 3. Results 3.1. Homogeneous and flamelet calculations Calculations of the ignition delay time in homogeneous systems are useful to understand the ignition behavior of turbulent sprays. For fixed values of the oxidizer and fuel stream temperatures (oxidizer stream: Y O2 ¼ 0:233, Y N2 ¼ 0:767. Fuel stream: Y C7 H16 ¼ 1), and assuming inert mixing between the two, the ignition delay of a perfectly stirred mixture has a minimum sid,0 corresponding to a specific mixture fraction value. This is referred to as the most reactive mixture fraction nMR [14], while sid,0 can be taken as a reference autoignition time. Comparison of the ignition delay of a spray with the corresponding sid,0 then gives indications on the influence of turbulence and of the spray parameters on the ignition process. The dependence of the ignition delay of an homogeneous mixture on n for the air/fuel temperatures of the sprays considered in this work is shown in Fig. 3. The quickest ignition is observed at rich mixture compositions, nMR = 0.166 (for n-heptane/air mixtures, nst = 0.062), and occurs at sid,0 = 0.84 ms. In these calculations, ignition was assumed to occur when the mixture temperature reached the threshold value of 1250 K. A similar criterion was used in the DNS simulations presented in the next

sections, with the ignition event being defined as the instant at which the gas-phase temperature reaches the value of 1250 K anywhere in the computational domain. The influence of mixing on non-premixed autoignition can also be quantified through the critical scalar dissipation rate N0,crit. This is the value of scalar dissipation that inhibits ignition in the laminar flamelet corresponding to the same ambient pressure and air/ fuel stream temperatures of the spray [14]. N0,crit was evaluated here by simulating the ignition of several transient flamelets corresponding to different values of N0, with N0 being the maximum scalar dissipation through the flamelet. In the calculations, the functional dependence between N and n in the flamelet equation was modelled according to the AMC model of O’Brien and Jiang [38], and unity Lewis number was assumed for all species to be consistent with the DNS. The computed dependence of sid on N0 for the air/fuel temperatures of the sprays considered in this work is also shown in Fig. 3. sid increases monotonically with increasing scalar dissipation. As N0,crit is approached, sid becomes very high, and eventually ignition no longer occurs. N0,crit is equal to 129 s1 here. A preliminary analysis of the ignition and subsequent flame propagation processes in the sprays studied in this work was performed by simulating the evolution of the corresponding laminar transient flamelet. In this calculation, N0 was kept constant and equal to 20 s1. This value corresponds to a ratio N0/N0,crit of 0.155, and is thus sufficiently far from the critical scalar dissipation rate. Results obtained from this simulation are shown in Fig. 4. Ignition is occurring at rich mixture compositions, as predicted by the homogeneous calculations. Once a burning kernel is established, the flame spreads quickly in mixture fraction space, and a fully burning solution is recovered shortly after. The flame propagation phase appears to be uneven; in particular, the flame moves faster towards leaner mixture compositions. This is a consequence of the increasing values of the flame temperature as nst is approached. Note the very long time needed to approach the steady-state burning flamelet (about twice the ignition time); this time, however, is a strong function of the scalar dissipation. 3.2. Main observations 3.2.1. General remarks The temporal evolution of the volume-averaged heat release rate for the different cases is shown in Fig. 5. To allow for a comparison between the different sprays investigated, the heat release rate of each case was divided by U0. The ignition sequence can be split into three different phases, each corresponding to a specific behavior of the mean heat release rate. We distinguish here between an inert mixing phase (the label IN will be used in the following), a pre-ignition phase (PI) and an ignition and flame propagation phase (IG). The inert mixing period (up to about 0.6 ms in Fig. 5a) is dominated by droplet evaporation and mixing of the resulting fuel vapor with the surrounding hot air. Heat release is negligible during this phase. Once a reactive mixture has been formed and the radical pool has grown, exothermic reactions start to be significant. This marks the beginning of the pre-ignition phase, which is initially characterized by a constant growth in the rate of heat release. As ignition is approached, the heat release rate stops growing and exhibits either a plateau or a sudden decrease to a local minimum, depending on the value of U0. This is followed shortly after by an abrupt increase occurring over a short period of time, which corresponds to the beginning of the ignition and flame propagation phase. This phase is characterized by the formation of several isolated ignition kernels and occurs after about 1.05 ms in Fig. 5b. These grow in size with time and, once they reach a critical size, they eventually start to merge, leading shortly to the combustion of the whole mixture. The present inert phase is

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

1259

Fig. 3. (a) Calculations of ignition delay time sid,0 for homogeneous n-C7H16/air mixtures. (b) Calculations of ignition delay time in unsteady n-C7H16/air laminar flamelets under constant scalar dissipation rate, parametrized by the maximum N0 across the layer. Data shown correspond to P = 24 bar, Tf = 450 K, Tair = 1000 K.

Fig. 4. Temperature versus mixture fraction for a transient autoigniting nonpremixed flamelet corresponding to P = 24 bar, Tf = 450 K and Tair = 1000 K. N0 was 20 s1, which corresponds to a ratio N0/N0,crit of 0.155.

conceptually similar to the conventionally-defined mixture formation phase often discussed in the diesel autoignition literature, such as in [39]. For the range of parameters examined here, this phase seems to have the same duration for all cases studied. The ignition delay time for all sprays investigated is about one millisecond, which is between 20% and 30% longer than sid,0(sid,0 = 0.84 ms). Therefore, the spontaneous ignition of a two-phase flow occurs later than the associated fastest homogeneous mixture, similarly to what was observed in DNS of autoigniting gaseous flows [8]. Turbulence promotes a faster ignition event, consistent with [8,14]. Less distinct conclusions can be made for the effects of inhomogeneities in the initial droplet size distribution. The pre-ignition phase appears to be faster for polydisperse sprays, with the corresponding peak in heat release rate being reached at earlier times. For fixed turbulence strength, however, ignition occurs earlier for those sprays that are initially monodisperse. A more detailed description of the phenomena occurring during the ignition sequence is presented in the following subsections. Due to the similarities observed between the cases investigated, the analysis is mostly restricted to Case A. Referenced instants of time, at which the numerical solution is investigated in more detail, are given in Fig. 5b.

Fig. 5. (a) Temporal evolution of volume-averaged heat release rate for Cases A, B, C and D in Table 1. (b) Detailed view of Case A, showing instants of time at which the solution is investigated further. Heat release rate of each case was divided by the corresponding value of U0 to allow for a comparison.

1260

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

3.2.2. Inert mixing phase The inert mixing phase is characterized by virtually zero fuel consumption, since the radical pool has just started to form. This can be observed in Fig. 6, which shows the scatterplots of heat release rate, temperature and n-C7H16, OH, H and CH2 O mass fractions plotted against mixture fraction at t = tIN = 0.6 ms. It is evident that OH and H are found in very small mass fractions (108 and 1011 respectively), while CH2O is seen at mass fractions of the order of 104. Contour plots of temperature, mixture fraction, scalar dissipation rate, heat release rate and HO2 and CH2O mass fractions at the same instant of time above are shown in Fig. 7. The plane x = 0.5L was chosen for visualizing the data. Since the sprays considered in this work are dilute, the evaporation process is spotty. This leads to the formation of small, cold fuel-rich islands surrounded by the hot ambient gas. The regions with the steepest gradients in mixture fraction are found near the droplets, where the fuel vapor is generated. Randomness in the initial conditions and turbulence may promote the formation of small clusters of drops, which are characterized by a more uniform mixture fraction distribution and lower values of the scalar dissipation rate. The complicated mixing pattern described above is responsible for the formation of a strongly inhomogeneous mixture, in which well-mixed spots alternate with poorly mixed ones. The most reactive regions start to develop close to either individual droplets or small groups of droplets, where sufficient fuel vapor is found. In particular, reactions during this phase are found to be most intense around n = 0.1. 3.2.3. Pre-ignition phase The pre-ignition phase starts after t = tIN = 0.6 ms, when the volume-averaged heat release rate becomes non-negligible. This phase is characterized by a slow temporal growth in heat release, culminating in a local maximum reached at t = tPI4 = 1.0 ms. Contour plots of temperature, heat release rate and selected species mass fractions at this instant of time are shown in Fig. 8. The mixture is still highly inhomogeneous, with few rich spots, corresponding to the locations of the evaporating droplets, surrounded

by leaner flow regions. Heat release reactions are again most intense in the proximity of either single or small groups of droplets. According to the droplet group combustion theory of Chiu and Liu [33], the combustion regime is intermediate between the single droplet and the internal group combustion modes. It was previously noted that the volume-averaged heat release rate decreases suddenly prior to ignition. To explain this behavior, the temporal evolutions of the mixture fraction PDF and of the conditional averages of the heat release rate and the mass fractions of selected species are shown in Fig. 9, the conditioning being on the mixture fraction. As time progresses, the consumption of oxygen and n-heptane becomes significant. The concentrations of the intermediate and final species, such as CH2O, CO and CO2 (CO2 is not shown in the figure), increase at all mixture fraction values, with their peak values shifting towards richer mixtures, and an increment in the conditional temperature is observed (temperature is not shown in the figure). The temporal evolution of the radical species, e.g., H, OH, HO2 and CH3, is slightly different. At early times, their conditional mass fractions increase for all mixture fractions. As ignition is approached, however, their concentration decreases for n between 0.05 and 0.15, while it continues to grow for richer mixture conditions. Since the lack of radicals inhibits the ignition reactions, the conditional heat release rate exhibits a similar temporal trend, decreasing where the radicals mass fractions are diminishing, and increasing otherwise. These results indicate that the sharp decrease in mean heat release rate prior to ignition is caused by a moderate decrease of the reaction intensity around n ’ 0.1, which is amplified by the relative abundance in the flow of regions having this composition. Figure 9 also shows that the value of n at which the heat release rate peaks undergoes a shift from about 0.1 to about 0.2 as ignition is approached, bracketing the reference value nMR. Ignition happens quite close to nMR, as will be discussed next. 3.2.4. Ignition phase The ignition phase begins at t = tIG1 = 1.075 ms, when the mean heat release rate starts to increase exponentially with time. This phase is characterized by local ignition and subsequent flame

Fig. 6. Scatterplots of heat release rate, temperature, and n-C7H16, OH, H and CH2O mass fractions against mixture fraction at t = tIN = 0.6 ms. Data shown correspond to Case A in Table 1.

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

1261

Fig. 7. Contour plots of temperature, mixture fraction, scalar dissipation rate, HO2 and CH2O mass fractions, and heat release rate corresponding to the axial plane x = 0.5L at t = tIN = 0.6 ms. Data shown correspond to Case A in Table 1.

Fig. 8. Contour plots of temperature, mixture fraction, scalar dissipation rate, HO2 and CH2O mass fractions, and heat release rate corresponding to the axial plane x = 0.5L at t = tPI4 = 1.0 ms. Data shown correspond to Case A in Table 1.

propagation. The ignition sequence is depicted in Fig. 10, which shows the temporal evolution of the T = 1250 K isosurface from t = tIG2 = 1.09 ms to t = tIG5 = 1.15 ms. The isosurface color corresponds to the local mixture fraction value. Ignition first occurs at few spatial locations where the mixture is locally rich

(n ’ 0.195). This value of mixture fraction is close to nMR(nMR = 0.166), suggesting that the most reactive mixture fraction concept continues to be valid also for sprays igniting at high pressure/low temperature conditions. As time progresses, the burning spots grow in size, expand toward leaner flow regions

1262

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

Fig. 9. Temporal evolution of conditional heat release rate, mixture fraction PDF and selected conditional species mass fractions. Referenced instants of time are given in Fig. 5b. Their values are: tPI1 = 0.70 ms, tPI2 = 0.80 ms, tPI3 = 0.90 ms, tPI4 = 1.00 ms, tPI5 = 1.05 ms. Data shown correspond to Case A in Table 1.

and eventually start to merge with each other. This process is accompanied by the appearance of independent, additional ignition kernels, which speed up the occurrence of combustion over the entire spray. These results are qualitatively similar to the findings of Wang and Rutland [16] for the autoignition process in n-heptane sprays at low pressures (P = 6 bar) and high oxidizer temperatures (Tair = 1100  1300 K) with 2D-DNS. Scatterplots of temperature and selected species mass fractions against n during ignition are shown in Figs. 11 and 12. At t = tIG1 = 1.075 ms, fuel and oxygen consumption around n = 0.195 has been considerable, and ignition appears to be imminent. Once

a flame kernel is formed, the reactants, the products, and some intermediate and radical species, among which OH, H2 and CH2O (OH is not shown in the figure), undergo a monotonic transition between their pre-ignition and burning states. Other species, however, exhibit a more complex behavior. This is the case of some radicals, such as HO2 and CH3 (CH3 is not shown in the figure), or intermediates like C4H8, whose concentrations first increase during the ignition transient, before being depleted quickly once fully burning conditions have been reached. The role of CH2O as a precursor to ignition in both single and two-phase flows has been discussed in detail in the experimental works of Gordon et al. [12] and

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

1263

Fig. 10. Temporal evolution of T = 1250 K iso-surface, colored according to the local mixture fraction. Referenced instants of time are given in Fig. 5b. Data shown correspond to Case A in Table 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

O’Loughlin and Masri [15]. The behavior of this species in the sprays considered in this section is fully consistent with these experimental findings, confirming the observed similarity between autoignition in single and two-phase flows. 3.3. Scalar dissipation at the autoigniting kernels Autoignition in single-phase flows has been found to occur first at those spatial locations where n ’ nMR and the scalar dissipation rate is low [8,10]. Similar behaviors have been observed for autoigniting n-heptane sprays [18], although the correlation existing between heat release and scalar dissipation rate was weaker than for purely gaseous mixtures. The picture emerging from the simulations presented in this work, while confirming some of these observations, reveals additional details of spray autoignition. Figure 13 shows the doubly-conditioned heat release rate at different instants, the conditioning being on the mixture fraction and the scalar dissipation rate, for Case A. It is clear that ignition is occurring first in regions with n = 0.195 and low values of scalar dissipation. The correlation existing between N and the heat release rate appears to be stronger than in previous works on spray autoignition [18]. This is indeed confirmed by Fig. 14a, which shows the temporal evolution of the conditional cross-correlation coefficient between these two quantities:

_ H  hx _ H jgiÞjgi hðN  hNjgiÞðx hrjgi ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ H  hx _ H ijgiÞ2 jgi hðN  hNjgiÞ2 jgi  hðx

ð31Þ

P _ a . Prior to ignition, hr—gi is negative around _ H ¼  a ha x where x n = 0.195, with values down to 0.7 being reached. These are lower than the minimum of 0.3 found by Schroll et al. [18] for autoigniting n-heptane sprays, although they remain higher than the value of 0.9 obtained for gaseous fuel autoignition [8]. The other cases

investigated in this work exhibit a similar behavior, as shown in Fig. 14b. These results indicate that the scalar dissipation rate is an important controlling parameter not only for the ignition of gaseous mixtures, but also for diesel sprays, and that the influence of N on sid is somewhat weaker for two-phase flows rather than for single-phase ones. The results presented so far indicates that ignition in two-phase flows is similar in many respects to autoignition in gaseous mixtures, at least for the range of conditions investigated in this work. Differences should be attributed to the mixture formation phase, which is the result of a complex interaction between turbulence and evaporation. In particular, the existence of nMR up to ignition may not be always guaranteed in two-phase flows, especially when the characteristic evaporation time is considerably lower than the reference ignition delay time. Under these conditions, the ratio between sid and sid,0 is expected to increase considerably. This situation has not been encountered in the present DNS database; however, it may appear when sid,0 is increased for a fixed value of sev, for example due to a decrease in the operating pressure or to a reduction of the oxygen concentration in the oxidizer stream [40]. 3.4. Parametric study 3.4.1. Influence of initial turbulence level Figure 15 shows the temporal evolution of the mean droplet evaporation rate, the volume-averaged scalar dissipation rate, the mean liquid temperature and the mean droplet diameter for Cases A and B. These simulations correspond to the same operating conditions, with the only difference being in the initial turbulence. Turbulent fluctuations in the background gas are responsible for faster droplet evaporation through augmented Shc, and lead to higher values of scalar dissipation during the transient droplet heating phase. As soon as the droplets reach their boiling

1264

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

Fig. 11. Scatterplots of temperature and n-C7H16, O2, H2 mass fractions against mixture fraction at different times during ignition. Referenced instants of time are given in Fig. 5b. Data shown correspond to Case A in Table 1.

temperature, these trends are reversed. The mean evaporation rate becomes larger for the low-turbulence case, due to the correspondingly larger size of the droplets. This shift is accompanied by a reversal in the observed trends for the volume-averaged scalar dissipation rate, with higher values being now associated with weaker turbulence conditions. The conditional scalar dissipation rate exhibits a similar behavior, as shown in Fig. 16. It is initially larger

for the high-turbulence case, before relaxing to lower values at later times. The role of turbulent fluctuations in promoting/delaying ignition in non-premixed systems has received considerable attention for both single- [8,41,42] and two-phase flows [17,43]. Turbulence affects the ignition process through the scalar dissipation rate. Depending on the particular set of conditions investigated, higher

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

1265

Fig. 12. Scatterplots of HO2, CH2O, C4H8 and CO2 mass fractions against mixture fraction at different times during ignition. Referenced instants of time are given in Fig. 5b. Data shown correspond to Case A in Table 1.

turbulence levels may result in lower values of scalar dissipation during most of the time interval preceding ignition. This, in turn, leads to shorter ignition delay times, as discussed in [8,14]. The behavior of the sprays that have been investigated in the present work is fully consistent with this physical picture. One should bear

in mind that the complex interaction existing between evaporation and mixing in sprays may lead to behaviors that are apparently different from those observed in gaseous flows. In particular, the role of u0 in promoting/delaying ignition may not be the same, as higher turbulence intensity promotes faster liquid evaporation, which

1266

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

Fig. 13. Doubly conditionally averaged heat release rate, plotted at different instants of time for the indicated range of scalar dissipation. Referenced instants of time are given in Fig. 5b. Data shown correspond to Case A in Table 1.

Fig. 14. (a) Temporal evolution of conditional correlation coefficient between heat release and scalar dissipation rate. Referenced instants of time are given in Fig. 5b. Their numerical values are: tPI2 = 0.80 ms, tPI3 = 0.90 ms, tPI4 = 1.00 ms, tIG3 = 1.10 ms. Data shown correspond to Case A in Table 1. (b) Conditional correlation coefficient between heat release and scalar dissipation rate at t = 1.00 ms for all cases listed in Table 1.

eventually results in better mixture formation due to the corresponding longer time available for the mixing processes to take place. In single-phase flows, turbulence seems to promote higher

values of hNjgi at first [42]; however, when u0 is large, well-mixed conditions are achieved quickly, and the conditional scalar dissipation rate then rapidly falls to zero.

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

1267

Fig. 15. Temporal evolution of (a) droplet mean evaporation rate and volume-averaged scalar dissipation rate, (b) droplet Sauter mean diameter and mean liquid temperature, and (c) droplet Stokes and modified Sherwood numbers. Data shown correspond to Case A and Case B in Table 1.

Fig. 16. Conditional scalar dissipation rate plotted at different instants of time. Data shown correspond to Case A (a) and Case B (b) in Table 1.

The literature on the evaporation of droplets in turbulent flows has been reviewed by Birouk and Gökalp [44]. The effect of turbulence on the evaporation rate of large hydrocarbon droplets (ad,0 = O(1) mm) was investigated either in the presence or the absence of a mean convective flow. With zero mean velocity of the carrier gas, increasing the strength of the turbulent velocity fluctuations always resulted in an enhancement of the droplet evaporation rate. This was attributed to the faster transport rate of the fuel vapor available in the proximity of the droplet surface. The picture that emerges from the simulations presented in this work is consistent with these results. The ability of a droplet to adapt to local variation of the carrier-phase velocity is given by its Stokes number St, which corresponds to the ratio between p the relaxation ffiffiffiffiffiffiffidroplet ffi time svd and the Kolmogorov time scale sg ¼ m= [45]. Turbulence can only affect the evaporation rate of those droplets whose characteristic relaxation time is large compared to the Kolmogorov time scale, so that a relative motion between the gaseous and the liquid phase can be established. As the initial turbulence strength is increased, sg diminishes, leading to an increase in the droplet Stokes number. This enhances the rate of heat and mass transfer between the gas and the droplet, resulting in faster evaporation, as seen in Fig. 15. As the droplet diameter decreases, its relaxation time becomes comparable to sg. St then falls below unity, and the droplet behaves as a particle tracer of the bulk

phase. No relative motion between the two phases exists anymore, and evaporation proceeds as in the absence of turbulence. 3.4.2. Influence of global equivalence ratio Previous investigations reported a longer ignition delay time for initially denser sprays [16]. Figure 17, which presents the temporal evolution of the mean droplet evaporation rate, the volumeaveraged scalar dissipation rate, the Sauter mean diameter and the mean liquid temperature for Case B and Case C in Table 1, suggests that the droplet evaporation is only slightly affected by U0. The mean properties of the liquid phase appear to be independent of the initial concentration of the spray. The volume-averaged scalar dissipation rate, however, is strongly affected by U0; in particular, higher values of the global equivalence ratio are responsible for higher values of N after t ’ 0.4 ms. This result appears to be in contrast with the fact that sid decreases with increasing U0. To solve this apparent contradiction, the profiles of P(g) and hNjgi for Case B and Case C have been plotted in Fig. 18 at two different instants of time prior to ignition. When plotted in mixture fraction space, no significative difference for the conditional scalar dissipation rate profiles is observed. However, the mixture fraction PDFs have different shapes; in particular, for the initially denser spray, a strong hump appears around n = 0.20, where h Njgi reaches its R1 maximum values. Since N ¼ 0 hNjgiPðgÞdg, this explains the

1268

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

Fig. 17. Temporal evolution of (a) droplet mean evaporation rate and volume-averaged scalar dissipation rate, and (b) droplet Sauter mean diameter and mean liquid temperature. Data shown corresponds to Case B and Case C in Table 1.

observed differences in the temporal evolution of the volumeaveraged scalar dissipation rate. These results also suggest that, as already noted by Schroll et al. [18], N is not a good indicator for predicting how delayed the ignition of the spray will be compared to the ignition of the corresponding homogeneous reactor, while hNjnMRi is. The observed differences in mixture composition also affect the ignition sequence. We first note that, differently from Cases A, B and D, the heat release rate of Case C reaches a plateau prior to ignition, instead of decreasing sharply (see Fig. 5a). Analysis of the evolution of the conditional species mass fractions, not shown here, reveals that the reaction intensity in the pre-ignition phase decreases for regions with n between 0.05 and 0.15. This behavior is identical to the one observed for Case A, and indicates that the different temporal evolution of the mean heat release rate is again a consequence of the different shapes of the mixture fraction PDFs. The abundance of flow regions where n ’ nMR in the dense spray case also affects the time needed to achieve ignition in the whole spray volume. This is revealed by Fig. 19, which shows the temporal evolution of the T = 1250 K isosurface for Case C at four selected instants of time during ignition. The spotty nature of the ignition event described for Case A is still evident; however, the number of independently burning kernels is considerably larger, due to the higher probability of finding flow regions with optimal ignition conditions.

Given the strong correlations existing between evaporation, mixing and ignition in sprays, it would be advantageous, from a practical viewpoint, to link the ignition delay time to a parameter describing the degree of stirring of the fuel/air mixture. Since ignition occurs first in those regions of the flow where N is low and n is close to its most reactive value, we consider here the PDF of the scalar dissipation rate at n = nMR, PSDR. Profiles of PSDR at t = 1.05 ms have been plotted in Fig. 20a for the different cases investigated in this work, while the correlation existing between PSDR and sid is shown in Fig. 20b. For all cases investigated, PSDR is maximum for low values of N, and decreases rapidly as N is increased. In general, sid increases if PSDR(0) is small, e.g., if there are few well-mixed nMR spots. The longest ignition delay (sid = 1.083 ms) occurs for Case A, for which PSDR(0) has the second lowest value (PSDR(0) = 0.072). Ignition is the fastest for Case C, with the corresponding value of PSDR(0) being the highest one (sid = 1.020 ms, PSDR(0) = 0.2284). Case B corresponds to intermediate conditions (sid = 1.035 ms,PSDR(0) = 0.1650). PSDR(0) is the lowest for Case D (PSDR(0) = 0.05), with the corresponding ignition delay time being the second longest one (sid = 1.052 ms). Cases B, C and D correspond to the same turbulence intensity of the bulk phase; additionally, they exhibit an almost identical correlation between heat release and scalar dissipation rate. For these cases, a negative linear dependence of sid on PSDR(0) is observed. These results suggest that, for sprays with similar turbulence levels, PSDR(0)

Fig. 18. Scalar dissipation rate (a) and mixture fraction PDF (b) at selected instants of time. Data shown correspond to Case B and Case C in Table 1.

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

1269

Fig. 19. Temporal evolution of T = 1250 K iso-surface, colored according to the local mixture fraction. Data shown correspond to Case C in Table 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 20. (a) PDF of scalar dissipation rate, PSDR, at n = nMR and t = 1.05 ms. (b) Ignition delay time plotted against PSDR(N = 0). Data shown for all cases listed in Table 1.

can provide an estimate on how delayed two ignition events will be with respect to each other. Case A is characterized by a lower turbulence than the other cases, which eventually result in a stronger correlation between heat release and scalar dissipation rate. This may partly explain why, for this spray, the point (PSDR(0), sid) does not satisfy the linear relation found for the high turbulence cases. Further investigation, covering a much wider range of operating conditions, is needed to extend these findings.

3.4.3. Influence of initial droplet size distribution In order to investigate the influence of spray polydispersity on autoignition, we performed a simulation (Case D in Table 1) where the initial diameter of the droplets was assigned according to a Rosin–Rammler distribution. The parameters of the Rosin–Rammler curve were chosen so that the initial Sauter mean diameter of the spray was aSMD = 25 lm and the droplet diameter varied between 5 and 40 lm. The resulting initial droplet size distribution is shown in Fig. 2.

1270

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

Previous numerical studies on single droplet combustion revealed that the ignition delay time is a strong function of the initial droplet diameter, and that changes in the operating conditions (e.g., pressure, initial air temperature) lead to different dependencies between sid and ad,0 [46]. A delaying effect of larger initial droplet size on spray ignition was observed experimentally by Gordon and Mastorakos [47] for atmospheric turbulent dilute diesel sprays, and numerically by Schroll et al. [18] for n-heptane droplets evaporating and reacting in a turbulent medium. An opposite trend was found by Wang and Rutland [17], where sid was seen to decrease with increasing ad,0. Differences in the liquid phase behavior due to spray polydispersity can be seen in Fig. 21, which shows the temporal evolution of the mean liquid temperature, the Sauter mean diameter, the residual liquid volume fraction / and the mean droplet evaporation rate for Case B and Case D. These simulations correspond to the same operating conditions, and only differ by the initial droplet size distribution. / is the ratio between the instantaneous and the initial value of the spray volume. The behavior of the polydisperse spray is initially dominated by the presence of the small droplets, which heat up and evaporate quickly. The mean liquid temperature and the mean droplet evaporation rate during this period have values similar to those observed for the monodisperse spray. At later times, however, the polydisperse spray starts to behave in a different way. In fact, once the smaller droplets have evaporated completely, the liquid phase properties depend on the behavior of the larger droplets only. Due to their large thermal inertia, these droplets heat up slowly; consequently, evaporation becomes less intense as compared to the monodisperse case, and the peak in _ is reached at later times. The mean evaporation rate becomes m larger again for the polydisperse spray as ignition is approached; this is a consequence of the slower temporal decrease of aSMD for the initially larger droplets, and of the quadratic dependence of _ on the droplet diameter. m The impact of the particular evaporation behavior of polydisperse sprays on the mixture fraction PDF, the conditional scalar dissipation rate and its standard deviation is shown in Fig. 22. During the inert mixing phase, evaporation is less intense for Case D, resulting in lower values of hNjgi. As ignition is approached, however, the situation changes, and the scalar dissipation is now lower _ in these for Case B, reflecting the different temporal evolution of m two simulations. The observed behavior of hNjgi is consistent with the findings of Schroll et al. [18], who showed that an increase in

the initial droplet size of a monodisperse spray leads to higher values of the conditional scalar dissipation rate in the whole mixture fraction space. Although the profiles of hNjgi for Case B and Case D exhibit large differences, the mixture fraction PDFs for these simulations are nearly the same at both instants of time considered, indicating that the initial droplet size distribution only plays a marginal role in determining the mixture composition. In the work of Schroll et al. [18], the shape of the mixture fraction PDF was found to be strongly affected by the initial value of the droplet diameter, with its peak shifting towards richer mixture compositions as a0 was increased. The lack of such a dependency here may be attributed to the simultaneous presence of both small and large droplets in the spray, whose effects on P(g) balance out. Despite the large values of hNjgi reached during the whole pre-ignition phase, ignition in Case D is still occurring at a relatively early time; this is due to the considerably larger fluctuations in scalar dissipation rate as compared to the monodisperse cases (see Fig. 22b), which indicate higher possibility of having some regions in the flow with low hNjnMRi where ignition occurs quickly. One should also note that the ignition delay time of the polydisperse spray is longer than those of the monodisperse ones with the same turbulence levels: this is consistent with the corresponding values of PSDR for these cases, as it can be observed in Fig. 20b. Scatterplots of temperature and selected species mass fractions are shown in Fig. 23. Ignition is now occurring over a wider interval of mixture fraction values as compared to the monodisperse sprays investigated in this work, with n ranging between 0.14 and 0.17. These values are close to nMR = 0.166, indicating that the concept of most reactive mixture fraction continues to hold also for polydisperse sprays. Another peculiar difference with the monodisperse spray cases lies in the presence of a considerable number of fuel pockets that have undergone little chemical reaction when ignition occurs. This is particularly evident from the scatterplots of n-C7H16 and CH2O mass fractions. These behaviors are a consequence of the differences in the evaporation rate of the small and large droplets in the spray. Small droplets evaporate quickly and are responsible for the formation of well-mixed spots. As time progresses, however, these regions continue to mix with the surrounding air and may become too lean for ignition to occur there. Those that remain rich enough during the whole pre-ignition phase eventually ignite. The maximum value of the mixture fraction in these regions depends on the initial size of the droplets that have evaporated in their neighborhood and this explains why

Fig. 21. Temporal evolution of (a) mean liquid temperature and Sauter mean diameter, and (b) residual liquid volume fraction and mean droplet evaporation rate, for different initial droplet size distributions. Data shown correspond to Case B and Case D in Table 1.

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

1271

Fig. 22. Conditional scalar dissipation rate (a), standard deviation of scalar dissipation (b) and mixture fraction PDF (c) at selected instants of time. Data shown correspond to Case B and Case D in Table 1.

ignition may now occur over a wider range of mixture fraction values. In fact, even though the value of nmax may differ from the one corresponding to the shortest ignition delay, early ignition in these regions is still possible as a result of the corresponding low values of the scalar dissipation during the pre-ignition phase. The evaporation of large droplets initially proceeds at a slow rate. However, as soon as they have been heated up, evaporation becomes substantial and, as discussed previously, this leads to the formation of very rich spots. The longer time needed for the formation of these regions and the strong fluctuations in the scalar quantities that characterize them explain why a considerable amount of fuel vapor has undergone only little chemical reactions as ignition is occurring elsewhere in the flow. This simultaneous presence of regions that have ignited and others that have remained almost unreacted is compatible with the notion that droplets of different sizes ignite at different times [18,46,47]. 3.5. Validation of the CMC method for spray autoignition 3.5.1. Mathematical formulation CMC has been successfully employed in the past to simulate a wide range of single-phase reacting flows [48–50]. Its application to spray combustion, however, has raised several questions. These are related essentially to the extension of CMC sub-models that were originally developed for single-phase flows to sprays, and to the treatment of terms describing droplet evaporation effects in the CMC equations [7]. The scope of this section is to exploit the DNS simulations listed in Table 1 to assess the CMC method for spray autoignition. This is done by comparing the conditionally-averaged DNS data against the numerical predictions obtained from the solution of the volume-averaged CMC equations:

  @Q a @2Q a @Q a _ a jgi þ daF  Q a  ð1  gÞ hPjgi þ h x ¼ hNjgi @t @ g2 @g

ð32Þ

Eq. (32) describes the temporal evolution of the conditionally averaged reactive species. A similar equation holds for the conditional temperature. hNjgi is the conditional scalar dissipation rate. This quantity is modelled here according to two different approaches, labelled as method A and method B in the following. In method A, hNjgi is closed according to the Amplitude Mapping Closure (AMC) model of O’Brien and Jiang [38]. This is the approach most commonly used in CMC simulations of both single [49,50] and two-phase flows [7,37,51]. The mixture fraction PDF, required by the AMC model, is constructed by computing the mean mixture fraction and its variance from the DNS data, and then assuming e which is also its shape to be described by a b-distribution. N, needed by AMC, is obtained by spatially averaging the DNS data.

In method B, the conditional scalar dissipation rate is extracted directly from the DNS data. hPjgi is the conditional evaporation rate term. This quantity is neglected in method A, and modelled in method B using the approach proposed by Borghesi et al. [7]:

qg hPjgi ¼

1 e gÞ V Pð

X _ d dðg  ns;d Þ m

ð33Þ

d

This model is based on the observation that evaporation occurs at the droplet surfaces only, where n = ns, and therefore the conditional evaporation rate term must be different from zero only at the saturation mixture fraction value. Additionally, it satisfies the condition R1  g is the conditional mean density, which is hPjgiPðgÞdg ¼ P. q 0 estimated from the ideal gas equation of state by assuming P to be constant in mixture fraction space. When hPjgi is closed using Eq. (33), estimates for the conditional scalar dissipation rate and the mixture fraction PDF at the saturation mixture fraction ns must be provided. This are obtained following [7]. The conditionally averaged chemical source term is closed at _ a ðY b ; TÞjgi ¼ x _ a ðQ b ; Q T Þ, first order in both approaches, e.g., hx where QT is the conditional temperature. The conditional heat reP _ H jgi ¼  a ha ðQ T Þx _ a ðQ b ; Q T Þ. Initializalease rate is modelled as hx tion of the conditional moments is performed by assuming an adiabatic frozen mixing distribution of the conditional scalars between the values corresponding to the oxidizer and fuel streams. 3.5.2. Main findings The validity of the AMC model to describe accurately scalar mixing processes in two-phase flows has often been questioned [52]. Similar concerns have been raised for the first-order modelling of chemical source terms in autoigniting mixtures [53]. To assess the validity of these approaches, we compare in Fig. 24 the _ H jgi against those obtained by modelled profiles of hNjgi and hx conditionally averaging the DNS data. The comparison is performed for Case A in Table 1 at selected instants of time prior to ignition. Similar results were obtained for the remaining DNS cases investigated in this work and are not shown here. We observe, however, that CMC was able to capture an earlier ignition delay time with increasing turbulence intensity and global equivalence ratio in the droplet laden region, consistent with the DNS results. This result was obtained for both CMC methods studied in this section. The modelled profiles of hNjgi closely follows the DNS for mixture fraction values up to 0.1. However, when richer mixture compositions are considered, large differences between these curves are observed, as shown in Fig. 24a. In particular, the modelled profiles overestimate the DNS ones around nMR, and is therefore expected that, for the cases considered here, usage of

1272

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

Fig. 23. Scatterplots of temperature and C7H16, O2, CH2O mass fractions against mixture fraction at different times during ignition. Data shown correspond to Case D in Table 1.

the AMC model instead of the DNS data will result in a longer ignition delay time. Additionally, we observe that, as a consequence of the evaporation process, hNjgi in sprays is not defined anymore over the entire mixture fraction space, and that AMC in its standard implementation is not able to capture this feature. Several methods

have been put forward to overcome this limitation of AMC [52]; their implementation in RANS or LES methods, however, is not straightforward, and therefore they have not been tested here. _ a jgi are shown in The limitations of a first-order closure for hx Fig. 24b. The comparison is limited here to the conditional heat

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

1273

Fig. 24. (a): comparison between conditional scalar dissipation rate as obtained from AMC (thick lines), and by conditionally averaging the DNS data (thin lines). (b): comparison between the conditional heat release rate as obtained using a first-order closure (thick lines), and by conditionally averaging the DNS data (thin lines). Data shown correspond to Case A in Table 1.

release rate. At times sufficiently far from the onset of ignition, fluctuations of the reactive scalars around their conditional means _ H ðQ b ; Q T Þ provides a good estimate of the condiare small, and x tional heat release rate. However, during the transition of the mixture between the unreacted and the burning state, fluctuations are _ H jgi is no longer no longer negligible, and a first-order closure of hx accurate. These results are in line with the findings of previous authors [53,54], who demonstrated that second-order closures of the chemical source terms in CMC are more appropriate in problems involving ignition and extinction phenomena. However, high-order closures are computationally expensive, and for this reason they have not been exploited in the CMC simulations presented in this section. To assess the capability of CMC to predict correctly ignition events in sprays, and to evaluate the impact that inaccurate modelling of the unclosed terms in the CMC equations has on the numerical predictions, a comparison between the results obtained from the CMC simulations and the DNS data is presented in Fig. 25. Data shown correspond to the temporal evolution of selected conditional moments at nig = 0.195 for Case A in Table 1. This corresponds to the mixture fraction value at which ignition is seen to occur in the DNS. CMC overpredicts the ignition delay time for both modelling strategies considered. The error is larger when modelling strategy A is considered, with sid being approximately 20% longer than the DNS value. Use of the hNjgi profiles extracted from the DNS data leads to a substantial improvement in the CMC results, highlighting once more the importance of an accurate modelling of this quantity. The discrepancy with the DNS is still in the order of 10%. This may be a consequence of the poor approximation of the chemical source terms during the onset of ignition that has been described earlier. It is interesting to observe that, at early times (e.g., t < 0.6 ms), CMC overpredicts the conditional temperature at nig, and that inclusion of the droplet source terms in the CMC equations does not lead to any change in the value of hTjnigi. This discrepancy is a consequence of the large separation in mixture fraction space existing between nig, where most of the pre-ignition reactions are occurring, and the saturation mixture fraction ns, at which droplets evaporate. In the CMC method, physical phenomena occurring at ns will affect the flow field at nig only if the model used for the scalar dissipation rate is able to accurately describe scalar mixing between ns and nig. The models used here do not provide such a description, and therefore effects due to droplet

evaporation do not propagate in mixture fraction space and remain confined at ns. It is clear that more accurate modelling of scalar mixing in the near droplet field is required in order to account correctly for droplet related effects in CMC. We finally observe that CMC predicts a two-stage ignition process, while single-stage ignition is found from the DNS data. Since this discrepancy occurs also when the conditional scalar dissipation rate is modelled according to the profiles extracted from the DNS simulations, we conclude that it is a consequence of the initially higher values of the conditional temperature in the CMC, and of the limitations of a first-order closure of the chemical source terms after ignition has started. 4. Conclusions Ignition in turbulent n-heptane sprays at P = 24 bar and Tair = 1000 K was studied by mean of DNS with complex chemistry. Ignition was characterized by a spotty pattern for all conditions investigated, with the ignition kernels developing first in those regions of the flow where the mixture fraction was close to the most reactive value, determined from homogenous autoignition calculations, and the scalar dissipation rate was low. The ignition delay time for the different cases investigated was between 1.020 ms and 1.083 ms, while the shortest possible time for homogeneous mixtures was 0.84 ms. The cross-correlation coefficient between heat release and scalar dissipation rate was approximately 0.7 around nMR for the low turbulence case, and approximately 0.5 for the remaining cases. These values, although not as low as for autoignition in single-phase flows, are lower than those found in previous works on autoigniting sprays. The spray operating conditions affected ignition through changes in the evaporation and mixing patterns. Initially higher turbulence intensity in the carrier gas was responsible for enhanced evaporation at first, leading to initially higher values of scalar dissipation. At later times, however, it promoted better air/fuel mixing, resulting in a reduction of the scalar fluctuation levels and, therefore, of the ignition delay time. The global equivalence ratio in the droplet laden layer was found to affect the number of locations at which ignition initially occurred. In particular, a higher value of U resulted in a larger number of ignition spots, due to the corresponding higher probability of finding regions where the mixture fraction was close to nMR. For the different sprays investigated in

1274

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

Fig. 25. Temporal evolution of selected conditional moments at n = nig = 0.195. Data shown correspond to Case A in Table 1. The conditional moments were obtained by conditionally averaging the DNS data (DNS) and from the numerical solution of the zero-dimensional CMC equations (CMC-A and CMC-B).

this work, changes in the initial turbulence level and the global equivalence ratio affected the PDF of the scalar dissipation rate at the most reactive mixture fraction, which provided a reliable indication of how delayed an ignition event will be compared to a reference case; in particular, higher values of PSDR(N = 0) corresponded to shorter ignition delay times. Ignition in the polydisperse spray occurred over a broader range of mixture fraction values as compared to the monodisperse spray cases. This result was interpreted as a consequence of the fact that, in the polydisperse spray, the mixing field is highly inhomogenous, with several poorly mixed rich spots alternating with leaner,

well-mixed ones. These variations in the fuel distribution are a direct consequence of the different evaporation behavior associated with the small and large droplets present in the spray. Despite this difference, the values of n at which ignition occurred were close to nMR, indicating that the concept of most reactive mixture fraction continues to be valid also for polydisperse sprays. Volume-averaged CMC simulations captured physical trends and provided reasonably accurate predictions of the ignition delay time. It was also shown, however, how accurate modelling of the unclosed terms in the CMC equations is fundamental for the success of the method. In particular, modelling of the scalar

G. Borghesi et al. / Combustion and Flame 160 (2013) 1254–1275

dissipation rate appears to be a critical point, especially concerning the description of scalar mixing in the near droplet field. Additionally, higher order closures of the chemical source terms may be necessary in those problems characterized by high levels of scalar dissipation rate fluctuations. Acknowledgments Professor J.B. Greenberg is acknowledged for many interesting discussions on the subject of spray autoignition. This work was sponsored by the EU Marie Curie Project MYPLANET. The authors kindly acknowledge the financial contribution received. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]

J. Dec, Proc. Combust. Inst. 32 (2009) 2727–2742. M. Yao, Z. Zheng, H. Liu, Prog. Energy Combust. Sci. 35 (2009) 398–437. N. Peters, Prog. Energy Combust. Sci. 10 (1984) 319–339. A. Klimenko, R. Bilger, Prog. Energy Combust. Sci. 25 (1999) 595–687. H.-W. Ge, E. Gutheil, Combust. Flame 153 (2008) 173–185. M. Mortensen, R. Bilger, Combust. Flame 156 (2009) 62–72. G. Borghesi, E. Mastorakos, C. Devaud, R. Bilger, Combust. Theory Modell. 15 (2011) 725–752. E. Mastorakos, T. Baritaud, T. Poinsot, Combust. Flame 109 (1997) 198–223. S. Sreedhara, K. Lakshmisha, Proc. Combust. Inst. 28 (2000) 25–34. S. Sreedhara, K. Lakshmisha, Proc. Combust. Inst. 29 (2002) 2051–2059. T. Echekki, J. Chen, Combust. Flame 134 (2003) 169–191. R. Gordon, A. Masri, E. Mastorakos, Combust. Flame 155 (2008) 181–195. R. Gordon, A. Masri, S. Pope, G. Goldin, Combust. Theory Modell. 11 (2007) 351–376. E. Mastorakos, Prog. Energy Combust. Sci. 35 (2009) 57–97. W. O’Loughlin, A. Masri, Flow Turbul. Combust. (2012) 13–35. Y. Wang, C. Rutland, Proc. Combust. Inst. 30 (2005) 893–900. Y. Wang, C. Rutland, Combust. Flame 149 (2007) 353–365. P. Schroll, A. Wandel, R. Cant, E. Mastorakos, Proc. Combust. Inst. 32 (2009) 2275–2282. C. Yoo, T. Lu, J. Chen, C. Law, Combust. Flame 158 (2011) 1727–1741. A. Viggiano, Combust. Flame 157 (2010) 328–340. A. Neophytou, E. Mastorakos, R. Cant, Proc. Combust. Inst. 33 (2011) 2135– 2142.

1275

[22] A. Neophytou, E. Mastorakos, R. Cant, Combust. Flame 159 (2012) 641–664. [23] B. Abramzon, W. Sirignano, Int. J. Heat Mass Transfer 32 (1989) 1605–1618. [24] B. Poling, J. Prausnitz, J. O’Connell, The Properties of Gases and Liquids, fifth ed., McGraw-Hill, 2000. [25] F. Mashayek, J. Fluid Mech. 405 (2000) 1–36. [26] H. Nomura, Y. Ujiie, H. Rath, J. Sato, M. Kono, Proc. Combust. Inst. 26 (1996) 1267–1273. [27] T. Dunstan, K. Jenkins, Int. J. Hydrogen Energy 34 (2009) 8389–8404. [28] C. Kennedy, M. Carpenter, Appl. Numer. Math. 14 (1994) 397–433. [29] O. Knio, H. Najm, P. Wyckoff, J. Comput. Phys. 154 (1999) 428–467. [30] C. Kennedy, M. Carpenter, R. Lewis, Appl. Numer. Math. 35 (2000) 177–219. [31] P. Brown, G. Byrne, A. Hindmarsh, J. Sci. Stat. Comput. 10 (1989) 1038–1051. [32] K. Harstad, J. Bellan, Int. J. Multiphase Flow 26 (2000) 1675–1706. [33] H. Chiu, T. Liu, Combust. Sci. Technol. 17 (1977) 127–142. [34] A. Lefebvre, D. Ballal, Gas Turbine Combustion, third ed., CRC Press, 2010. [35] S. Liu, J. Hewson, J. Chen, H. Pitsch, Combust. Flame 137 (2004) 320–339. [36] H. Ciezky, G. Adomeit, Combust. Flame 93 (1993) 421–433. [37] Y. Wright, O. Margari, K. Boulouchos, G.D. Paola, E. Mastorakos, Flow Turbul. Combust. 84 (2010) 49–78. [38] E. O’Brien, T. Jiang, Phys. Fluids A 3 (1991) 3121–3123. [39] J. Heywood, Internal Combustion Engine Fundamentals, first ed., McGraw-Hill, 1988. [40] G. Borghesi, Autoignition in Turbulent Two-Phase Flows, Ph.D. Thesis, University of Cambridge, 2012. [41] J. Blouch, C. Law, Combust. Flame 132 (2003) 512–522. [42] C. Markides, E. Mastorakos, Proc. Combust. Inst. 30 (2005) 883–891. [43] H. Koss, D. Brüggemann, A. Wiartalla, H. Bäcker, A. Breuer, Investigations of the Influence of Turbulence and Type of Fuel on the Evaporation and Mixture Formation in Fuel Sprays, Final Report of JOULE Project on Integrated Diesel European Action (IDEA), 1992. [44] M. Birouk, I. Gökalp, Prog. Energy Combust. Sci. 32 (2006) 408–423. [45] J. Rèveillon, F. Demoulin, J. Fluid Mech. 583 (2007). [46] R. Stauch, S. Lipp, U. Maas, Combust. Flame 145 (2006) 533–542. [47] R. Gordon, E. Mastorakos, Exp. Therm. Fluid Sci. 43 (2012) 40–46. [48] I. Kim, E. Mastorakos, Proc. Combust. Inst. 30 (2005) 911–918. [49] A. Triantafyllidis, E. Mastorakos, R. Eggels, Combust. Flame 156 (2009) 2328– 2345. [50] A. Garmory, E. Mastorakos, Proc. Combust. Inst. 33 (2011) 1673–1680. [51] Y. Wright, G.D. Paola, K. Boulouchos, E. Mastorakos, Combust. Flame 143 (2005) 402–419. [52] M. Zoby, S. Navarro-Martinez, A. Kronenburg, A. Marquis, Int. J. Heat Fluid Flow 32 (2011) 499–509. [53] G.D. Paola, I. Kim, E. Mastorakos, Flow Turbul. Combust. 82 (2009) 455–475. [54] S. Sreedhara, K. Huh, Combust. Flame 143 (2005) 119–134.