Direct numerical simulation of ignition front propagation in a constant volume with temperature inhomogeneities

Direct numerical simulation of ignition front propagation in a constant volume with temperature inhomogeneities

Combustion and Flame 145 (2006) 145–159 www.elsevier.com/locate/combustflame Direct numerical simulation of ignition front propagation in a constant ...

637KB Sizes 0 Downloads 62 Views

Combustion and Flame 145 (2006) 145–159 www.elsevier.com/locate/combustflame

Direct numerical simulation of ignition front propagation in a constant volume with temperature inhomogeneities II. Parametric study Evatt R. Hawkes ∗ , Ramanan Sankaran, Philippe P. Pébay, Jacqueline H. Chen Reacting Flow Research Department, Combustion Research Facility, Sandia National Laboratories, P.O. Box 969 MS 9051, Livermore, CA 94551-0969, USA Received 18 December 2004; received in revised form 20 August 2005; accepted 25 September 2005 Available online 5 January 2006

Abstract The influence of thermal stratification on autoignition at constant volume and high pressure is studied by direct numerical simulation (DNS) with detailed hydrogen/air chemistry. Parametric studies on the effect of the initial amplitude of the temperature fluctuations, the initial length scales of the temperature and velocity fluctuations, and the turbulence intensity are performed. The combustion mode is characterized using the diagnostic measures developed in Part I of this study. Specifically, the ignition front speed and the scalar mixing timescales are used to identify the roles of molecular diffusion and heat conduction in each case. Predictions from a multizone model initialized from the DNS fields are presented and differences are explained using the diagnostic tools developed.  2005 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Ignition; Direct numerical simulation; HCCI; Multizone model

1. Introduction Homogeneous charge compression ignition (HCCI) engines may provide efficiency gains over conventional spark-ignition engines due to a higher compression ratio and no throttling losses. Also improvements in emissions characteristics relative to a diesel engine are possible due to more complete mixing eliminating areas of excessively high temperatures or fuel concentration, and hence reducing emissions of NOx and particulates. However, control of the heat release rate becomes more difficult since the * Corresponding author. Fax: +1 925 294 2595.

E-mail address: [email protected] (E.R. Hawkes).

combustion process can become kinetically limited rather than mixing limited. For a purely homogeneous charge, options for control by varying the composition of the charge are limited, and there is considerable interest at present in the possibility of using a stratified charge in order to spread the heat release rate out over several crank angle degrees [1]. In addition, due to turbulence within the cylinder and heat transfer to and from the cylinder wall, a certain level of stratification is inevitable. At present, little is known about the fundamental combustion modes that can occur in weakly stratified charges. This knowledge may enable the development of mixing strategies to tailor the ignition timing and duration. Furthermore identifying the governing parameters would assist in the development of models for HCCI combustion.

0010-2180/$ – see front matter  2005 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.combustflame.2005.09.018

146

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

In this fundamental study direct numerical simulation is performed of the autoignition of a lean homogeneous H2 /air mixture with thermal stratification in order to shed light on the nature of the combustion mode, identify the controlling parameters, and comment on modeling, in particular multizone modeling [2–4], which has been demonstrated to provide very good predictions of many relevant quantities at HCCI conditions. In Part I of this study [5], numerical diagnostics were developed for the delineation of the combustion mode. Following theoretical arguments of Zel’dovich [6], the ignition front speed was identified as a key controlling parameter of the combustion mode. The utility of tracking the ignition front speed to understand the effects of diffusive transport within propagating ignition fronts was demonstrated and verified by comparison with other measures based on the temperature gradient [7] and an analysis of the importance of reactive and diffusive terms in the chemical species evolution equations. A passive scalar mixing timescale based on enthalpy fluctuations was identified as a controlling parameter for the effects of molecular mixing altering the distribution of temperature. Using a single 2D DNS case, the techniques were demonstrated. Predictions of the multizone model initialized from the DNS were compared with the DNS and explained using the diagnostics developed. This paper forms Part II of the study. Here, parametric studies are conducted by varying the initial root-mean-square (RMS) temperature fluctuations, the most energetic length scales of the initial temperature and velocity spectra, and the turbulent velocity fluctuations. The numerical diagnostics developed in Part I are employed to understand the influence of these parameters on the observed combustion mode. The application of these diagnostics in a broad range of conditions serves to confirm their value and robustness and develop confidence in future applications in less idealized settings. Multizone model predictions are presented for each case and the observed results are explained with reference to the diagnostic tools developed. It is worth noting at this point the goals with respect to the multizone model. The results shown here are not intended to provide validation of the model under HCCI engine conditions. Such validation has already been provided in comparisons with experiment. Rather, the work is intended to develop tools to define the valid regime, and explain why good predictions have been obtained. Thus, many of the conditions simulated here are deliberately beyond the intended realm of application. By applying the model in such regimes, where mixing effects are important, much will be learned about how these effects will degrade the results and how to predict a priori whether they are significant. Also, since the model does not include these effects, com-

parison with the DNS serves as an independent check on whether the developed diagnostics can predict robustly different combustion modes.

2. Numerical implementation and initial conditions Two-dimensional direct numerical simulations of the autoignition of a thermally stratified fuel-lean homogeneous H2 –air mixture are performed using a detailed reaction mechanism [8]. The turbulent flow field is prescribed by an initial turbulent kinetic energy spectrum function [9] and a random temperature field superimposed on a mean temperature field. The temperature spectrum, which is similar to the kinetic energy spectrum, is used to specify the characteristic scales of the initial exothermic spots. The numerical technique is described in Part I of this study [5]. Parametric studies on the effect of the initial amplitude of the temperature fluctuations, the initial length scales, and the turbulence intensity are performed. In all cases the temperature fluctuations are superimposed on a uniform mean temperature of 1070 K at 41 atm, and the equivalence ratio is initially uniform and equal to 0.1. The 2D DNS used in the first part of this study [5] and as Case A in [7] forms a baseline case about which the parametric studies are conducted. This case, which will be referred to here as Case A, has an RMS temperature fluctuation of 15.0 K and an initial most energetic length scale of the temperature fluctuations of 1.32 mm. The initial turbulence intensity is 0.5 m/s and the initial most energetic length scale of the turbulence fluctuations is 1.0 mm. The first parametric study investigates the effect of the initial amplitude of the temperature fluctuation. Case A is used as a baseline and three other cases are simulated, having RMS amplitudes of 3.75, 7.5, and 30 K. Otherwise parameters are identical to Case A. These cases are referred to in the paper by these amplitudes as T  = 3.75, T  = 7.5, and T  = 30.0. In the context of this parametric study, Case A is referred to as T  = 15.0 for clarity. Note that in these cases, it is possible to use an identical initial temperature fluctuation field, which is then scaled by the appropriate amount according to the specified amplitude—this eliminates statistical errors associated with random initialization in the comparison between cases. In the second parametric study, the effect of the length scale of the temperature fluctuations is investigated in the absence of initial turbulent velocity fluctuations. Starting from the most energetic length scale of Case A, 1.32 mm, the length scale is then doubled in a second case, and doubled again in a third. These cases are referred to as u0L = 0.41, u0L = 0.82, and

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

147

Table 1 Initial turbulence and grid parameters used in DNS cases

Case A

T (K) 15.0

T  = 3.75 T  = 7.5 T  = 30.0

3.75 1.32 7.50 1.32 30.0 1.32

Case

LTe (mm) 1.32

u (m/s) 0.5

Lue (mm) 1.0

L (cm) 0.41

0.5 0.5 0.5

1.0 1.0 1.0

0.41 0.41 0.41

720 720 720

N 960

u0L = 0.41 15.0 u0L = 0.82 15.0 u0L = 1.64 15.0

1.32 2.64 5.28

0.0 0.0 0.0

– – –

0.41 0.82 1.62

960 720 720

L = 0.82 L = 1.64

2.64 5.28

1.0 2.0

2.0 4.0

0.82 1.64

1440 1920

15.0 15.0

u0L = 1.64, where u0 refers to the fact that there are no turbulent velocity fluctuations, and L is the size of computational domain in cm, which is also doubled in conjunction with the most energetic length scale. In the third parametric study, the effect of the length scale of the temperature fluctuations is investigated in the presence of turbulent velocity fluctuations. When changing the length scale of scalar fluctuations, it is also logical to change the length scale of the turbulent velocity fluctuations. However, some quantity must be kept fixed. Here, the turbulent integral time scale is chosen. Therefore, in this parametric study, since the length scale is increased, but time scale is to be kept constant, the turbulence intensity is increased. This also leads to an increase in the Reynolds number. Starting from Case A, the length scale is then doubled in a second case, and doubled again in a third. These cases are referred to as L = 0.82 and L = 1.64. For clarity Case A is referred to as L = 0.41 in the context of this parametric study. The final parametric study is on the effect of turbulence fluctuations, comparing cases with and without turbulence. Table 1 summarizes the relevant parameters for each case. The RMS temperature fluctuation is given by T  , the most energetic length scale of temperature fluctuations by LTe , the turbulence intensity by u , the most energetic length scale of turbulence fluctuations by Lue , the size of the computational domain in both directions by L, and the number of grid nodes in each direction by N . The resolution requirement differed somewhat between cases and was determined by comparison with solutions on coarser grids. 3. DNS results—variation of T  The first parametric study examines the influence of the amplitude of the temperature fluctuations—

Fig. 1. Normalized mean heat release rate versus time for Cases T  = 3.75, 7.5, 15.0, and 30.0 K. Table 2 Reference values corresponding to homogeneous ignition delay at 1070 K and 41 atm for a hydrogen–air mixture at equivalence ratio 0.1 τ0 (ms)

HR0 (MJ/m3 /s)

2.9

2.9408 × 104

all other parameters remain fixed. The case study is particularly interesting because while the change in amplitude affects the temperature gradient and thus affects the importance of deflagration in the domain, it does not affect the volumetric-averaged instantaneous mixing time-scale, defined as τmix =

h 2 , 2α|∇h|2

(1)

where h is the mixture enthalpy, and h is the RMS enthalpy fluctuation. Hence the importance of passive scalar dissipation is unchanged between cases. Fig. 1 shows the mean heat release rate as a function of time for the different cases, referred to by the level of the initial temperature fluctuation. In this figure and throughout the paper, the heat release rate is normalized by the maximum heat release rate of the constant volume ignition at the mean temperature, HR0 , and time is normalized by the homogeneous ignition delay time, τ0 , defined as the time to the point of maximum heat release rate summarized in Table 2. Two points are immediately obvious. First, the larger temperature fluctuations spread out the heat release rate and reduce its peak. A less obvious result is that the peak heat release rate is advanced for larger T  . This is primarily because early-igniting regions compress the later-igniting regions, accelerating the onset of ignition. A secondary reason is that effects of molecular transport within propagating fronts also accelerate the combustion of later igniting kernels. These results demonstrate that there is a strong influence of the initial temperature distribution on the timing and duration of combustion. For practical applications, the result also implies that increasing in-cylinder temperature fluctuations near top dead center would be an

148

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

effective way of smoothing the heat-release profile [1,10]. Next, the nature of the combustion modes for the different cases is assessed. Fig. 2 shows isocontours of heat release rate at the time of maximum heat release rate for each case. It is evident that T  has a large effect on the character of the combustion. Lower values of T  result in a more homogeneous, volumetric combustion process, while larger values tend to promote frontlike structures. Applying the diagnostics developed in Part I [5], the density-weighted mean displacement speed, sd∗ [11], is given by sd∗ =

ρ DYH2 /Dt , ρ0 |∇YH2 |

(2)

where D/Dt is the material derivative and YH2 is the hydrogen mass fraction. The mean displacement speed on the isosurface, YH2 = 8.5 × 10−4 , is extracted. In Part I [5], it was shown that this isosurface coincides approximately with the location of maximum heat release rate. The temporal evolution of this quantity for each case is presented in Fig. 3. The curves all show the same qualitative U shape that was observed in Part I. As expected, the minimum speed decreases as T  increases, due to the higher initial temperature gradients and therefore lower ignition front speeds. It is interesting to note that Case T  = 15.0 K and Case T  = 30.0 K have similar speeds at one point in the simulation. This points to the possible importance of molecular transport effects in the front limiting the minimum attainable speed. These results can be compared with corresponding 1D test cases in Fig. 4. These 1D test cases are initialized from a sinusoidal temperature profile, with a wavelength of 1.5 mm, close to the most energetic length scale in the DNS. The RMS temperature fluctuation is also varied, as in the DNS study. Qualitatively, the results are similar, validating the use of the simplified 1D test cases to understand the characteristics of combustion for these conditions. The fraction of the front length flagged as deflagrations using the criterion sd∗ /sL < 1.1 as in Part I [5] is compared between the different cases in Fig. 5. It was shown in Part I [5] that the choice of the value of 1.1 is a good indicator of the transition from deflagration to spontaneous ignition, although, of course, transition occurs more gradually. Case T  = 3.75 has no flame elements propagating in a diffusion-controlled regime—it burns purely by spontaneous ignition front propagation. Comparing Cases T  = 7.5, T  = 15.0, and T  = 30, it is noted that the amount of deflagration observed is strongly dependent on the thermal stratification, and increases with increased thermal stratification. Fig. 6 shows the fraction of production rate of area of “burnt” gases (having YH2 <

Fig. 2. Normalized heat release rate isocontours for Cases T  = 3.75, 7.5, 15.0, and 30.0 K, shown in this order from top to bottom (rainbow color scale—white denotes maximum heat release rate). The size of the computational domain is 4.1 × 4.1 mm.

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

Fig. 3. Density-weighted displacement speed sd∗ versus time for DNS Cases T  = 3.75, 7.5, 15.0, and 30.0 K.

Fig. 4. Density-weighted displacement speed sd∗ versus time for 1D reference cases having T  = 3.75, 7.5, 15.0, and 30.0 K.

8.5 × 10−4 ) due to deflagrations, according to the criterion. If dl is an incremental segment length associated with a given contour segment, the rate of production of burnt gases is calculated by summing the quantity sd dl for isosurface segments flagged as deflagrations, divided by the total for all segments. The values for the fraction of burnt area production rate are lower than the fractions of front length because the speed of the front is lower, as described in Part I. Only Cases T  = 15.0 and T  = 30 have significant burnt area production due to deflagrations. The other controlling parameter for the importance of molecular diffusion effects is the mixing timescale, as identified in Part I. The mixing timescale is plotted versus time for these cases in Fig. 7. At the beginning of the simulation, since the initial temperature fluctuation fields differ only by a scale factor,

149

Fig. 5. Fraction of front length identified as a deflagration by the criterion sd∗ /sL < 1.1 for Cases T  = 3.75, 7.5, 15.0, and 30.0 K.

Fig. 6. Fraction of burnt area production rate identified as a deflagration by the criterion sd∗ /sL < 1.1 for Cases T  = 3.75, 7.5, 15.0, and 30.0 K.

Fig. 7. Mixing timescales for Cases T = 3.75, 7.5, 15.0, and 30.0 K.

150

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

Fig. 8. Normalized mean heat release rate versus time for the DNS and the multizone model started from the DNS initial condition (MZ 0%) and at 10% pressure rise (MZ 10%): (a) Case T  = 3.75 K, (b) Case T  = 7.5 K, (c) Case T  = 15 K, (d) Case T  = 30 K.

the timescales are identical. Later they diverge due to the different combustion process encountered in each case, but always remain comparable. The spread of mixing timescales at the point of maximum heat release rate in each case is approximately 15%. The mixing timescales are initially comparable to the ignition delay but decrease with time due to the effects of turbulent straining. Hence passive scalar dissipation effects are projected to be significant and nearly equally so in all cases. However, it was already observed that the importance of deflagrations differed between the cases. Having examined the front speeds and mixing timescales, we now present the predictions of the multizone model for each case. As in Part I, the multizone calculations are initialized from the DNS field at the initial time of the simulation (denoted MZ 0%) and from the time of 10% pressure rise (denoted MZ 10%). Figs. 8a–8d show the multizone predictions for the cases with different amplitudes of temperature fluctuations. In addition, predictions for the burn duration (dt), defined as the time between the 10% pressure rise and the 90% pressure rise, and the predictions of the maximum heat release rate averaged over the domain are summarized in Table 3. Overall, the predictions starting from the initial condition of the DNS are poor in all cases, worsening for cases

with larger temperature gradient. The large discrepancy is a result of not accounting for changes in the effective PDF of the initial temperature fluctuations due to scalar mixing. This is evident because it was found that for Cases T  = 3.75 K and T  = 7.5 K deflagrations are insignificant. When the calculations are started from a 10% pressure rise, and the most significant scalar mixing has already occurred, the predictions improve in all cases. Note that in practical multizone applications, much of the scalar mixing before ignition is accounted for in a prior CFD calculation, and so perhaps these cases are more representative of a practical application. For Case T  = 3.75 K, the prediction is very good. However, it is noted that the prediction progressively deteriorates as the temperature fluctuation is increased. This is attributed to the increasing importance of deflagrations in the calculations. It is noted that the magnitude of the errors is qualitatively consistent with the predicted fraction of deflagration. Conditions having significant deflagrations are outside of the intended valid region of the multizone model and of course then good performance is not expected. Since good performance has been obtained with the multizone model under practical HCCI conditions, the present results can be used to infer that deflagrations were not significant at those conditions.

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

151

Table 3 Comparison of burn duration (dt), defined as the time between the 10% pressure rise and the 90% pressure rise, and normalized heat release rate (HRmax /HR0 ) averaged over the domain for DNS with multizone model for the parametric study in T  Case

T  = 3.75

T  = 7.5

T  = 15.0

T  = 30

T

3.75

7.50

15.0

30.0

DNS dt (ms) MZ-0% dt (ms) MZ-0% dt error (%) MZ-10% dt (ms) MZ-10% dt error (%)

0.51 0.56 9.80 0.52 1.96

0.54 0.71 31.5 0.56 3.70

0.64 1.14 78.1 0.72 12.5

0.95 2.10 121 1.40 32.1

DNS HRmax /HR0 MZ-0% HRmax /HR0 MZ-0% HRmax /HR0 error (%) MZ-10% HRmax /HR0 MZ-10% HRmax /HR0 error (%)

0.76 0.48 37.4 0.71 6.84

0.55 0.29 48.1 0.50 9.63

0.34 0.16 54.8 0.29 16.4

0.18 0.08 58.0 0.12 33.5

Fig. 9. Normalized heat release rate isocontours: Cases u0L = 0.41, u0L = 0.82, and u0L = 1.64 (rainbow color scale—white denotes maximum heat release rate).

It is interesting to note that the multizone model always predicts a longer burn duration, later maximum heat release rate timing, and lower peak heat release rate relative to the DNS. This fact may be useful in determining whether mixing effects are important in practical applications of the multizone model.

4. DNS results—variation of LTe without turbulence In the previous parametric study, the length scale was kept fixed while T  was varied. This resulted in a fixed initial scalar dissipation timescale. In the present study, the initial most energetic length scale of the

temperature fluctuation is varied while the amplitude of the fluctuation is fixed, so that both the importance of passive scalar dissipation and the prevalence of transport effects in fronts changes in the simulations. The initial turbulent velocity fluctuations are zero in these cases. Fig. 9 shows contours of heat release rate for Cases u0L = 0.41, u0L = 0.82, and u0L = 1.64. The images are shown approximately to scale relative to each other. It is evident that smaller length scale cases lead to thinner front structures on average in the domain. The evolution of the mean front speed is presented for the cases with different initial temperature stratification length scales in Fig. 10. As expected, the mean

152

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

Fig. 10. Density-weighted displacement speed sd∗ versus time for Cases u0L = 0.41, u0L = 0.82, and u0L = 1.64.

front speed tends to decrease with decreasing length scale, due to the increasing mean temperature gradient. Note also that larger length scale cases, having lower scalar dissipation, tend to ignite earlier. These trends are consistent with the 1D parametric study performed in Part I [5]; see in particular Fig. 2. Note also that both the initial ignition and final consumption points have a zero initial temperature gradient and the speed becomes unbounded.1 The fraction of front length flagged as a deflagration using the criterion sd∗ /sL < 1.1 is shown in Fig. 11, and the fraction of burnt area production in Fig. 12. These figures reveal that only in Case u0L = 0.41 are deflagrations very significant—in the other cases temperature gradients are too low to produce any deflagration. Fig. 13 shows the evolution of the mixing timescale for the cases in this parametric study. Since there are no imposed turbulence fluctuations, the mixing time scales do not vary strongly in time. For Case u0L = 0.41, the mixing time scale is on the order of the ignition delay time and passive scalar mixing effects will be important, whereas the other cases have much longer mixing time scales, and passive scalar mixing is not significant. Multizone model predictions are presented for the cases where the initial length scale of the temperature fluctuation is varied in Figs. 14a–14c. The multizone calculations are initialized from the DNS initial condition. Results for the burn duration defined at maximum heat release rate are summarized in Table 4. The predictions for Cases u0L = 0.82 and u0L = 1.64 are very good, as expected, since both deflagrations 1 As explained in Part I [5], this does not imply that the reactant consumption rates are unbounded.

Fig. 11. Fraction of front length having sd∗ /sL < 1.1 for Cases u0L = 0.41, u0L = 0.82, and u0L = 1.64.

Fig. 12. Fraction of burnt area production having sd∗ /sL < 1.1 for Cases u0L = 0.41, u0L = 0.82, and u0L = 1.64.

Fig. 13. Mixing timescales for Cases u0L = 0.41, u0L = 0.82, and u0L = 1.64.

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

153

Fig. 14. Normalized mean heat release rate versus time for the DNS and the multizone model started from the DNS initial condition (MZ): (a) Case u0L = 0.41, (b) Case u0L = 0.82, (c) Case u0L = 1.64. Table 4 Comparison of burn duration (dt) and normalized heat release rate (HRmax /HR0 ) for DNS with multizone model for the parametric study in LTe without turbulence Case

u0L = 0.41

u0L = 0.82

u0L = 1.62

LTe (mm)

1.32

2.64

5.28

DNS dt (ms) MZ-0% dt (ms) MZ-0% dt error (%)

0.90 1.14 26.7

1.05 1.12 6.67

1.15 1.17 1.74

DNS HRmax /HR0 MZ-0% HRmax /HR0 MZ-0% HRmax /HR0 error (%)

0.20 0.16 22.8

0.17 0.16 6.89

0.15 0.15 0.96

and passive scalar mixing are unimportant. This is encouraging for the multizone model since these cases contain integral length scales more representative of typical engine length scales. The agreement in Case u0L = 0.41 is slightly worse due to the significance of transport within ignition fronts and passive scalar mixing. Overall, this parametric study offers no surprises; however, the inclusion of turbulence fluctua-

tions, presented in the next section, will result in some different conclusions.

5. DNS results—variation of LTe with turbulence The final parametric study consists of variation of the integral scale of the turbulence and temperature

154

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

Fig. 15. Normalized heat release rate isocontours: Cases L = 0.41, L = 0.82, and L = 1.64 (rainbow color scale—white denotes maximum heat release rate).

fluctuations. The turbulence integral time scale is kept constant and therefore the turbulent velocity fluctuation is changed by the same factor as the length scale. This also results in a change in the Reynolds number by the square of the factor by which length scale is changed. The initial temperature fields are identical to the corresponding length scale cases without turbulence. This allows an additional comparison of the effects of turbulent straining. Fig. 15 shows the heat release rate contours for Cases L = 0.41, L = 0.82, and L = 1.62. The images are shown approximately to scale relative to each other. Consistent with increasing Reynolds number, there is a larger range of scales for larger length scales. In contrast to the cases without turbulence fluctuation, there is no obvious pattern of decreasing thickness of the fronts, on average, as the length scale is decreased. The turbulent velocity fluctuations, particularly in the higher Reynolds number cases, strain the initially nearly isotropic hot spots into highly elongated filamentary structures, increasing the temperature gradient along the long edges, resulting in significant small-scale temperature fluctuations in all cases. The density-weighted mean front speed is shown for the different length scale cases in Fig. 16. Although the trend of decrease of the initial ignition delay is reproduced, with cases having larger length scales igniting earlier, the result during the main portion of heat release rate is very different from the results in Fig. 10 for the cases without turbulence. Here, despite initially larger length scales, the min-

Fig. 16. Density-weighted displacement speed sd∗ versus time for Cases L = 0.41, L = 0.82, and L = 1.64.

imum flame speeds are quite comparable for all the cases. The fraction of front length flagged as a deflagration is plotted versus time in Fig. 17, and the fraction of burnt gas production in Fig. 18. All cases show roughly equal importance of deflagrations, with the large length scale cases even at times showing more importance of deflagration. These surprising results deserve further investigation. In Part I [5], we showed that the temperature gradient was a major factor in determining the speed of the front. Fig. 19 shows the mean temperature gradient on the surface for Cases L = 0.41, L = 0.82, and L = 1.64. Despite the fact that the initial mean tem-

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

155

Fig. 17. Fraction of front length having sd∗ /sL < 1.1 for Cases L = 0.41, L = 0.82, and L = 1.64.

Fig. 19. Mean spatially averaged magnitude of temperature gradient (|∇T |) versus time for Cases L = 0.41, L = 0.82, and L = 1.64.

Fig. 18. Fraction of burnt area production having sd∗ /sL < 1.1 for Cases L = 0.41, L = 0.82, and L = 1.64.

Fig. 20. Mixing time scales for Cases L = 0.41, L = 0.82, and L = 1.64.

perature gradients in the domain differed by a factor of 2, all three cases show qualitatively the same values, especially near the time of peak heat release rate, around 0.9τ0 in Case L = 0.41. This can be explained by considering the effect of turbulence modifying the temperature gradients, and is most easily understood by examining the mixing time scales, which are presented for these three cases in Fig. 20. Initially there is a factor of 4 difference in mixing time scale between Cases L = 0.41 and L = 0.82, and a further factor of 4 between Cases L = 0.82 and L = 1.64. However, the effects of turbulence straining rapidly increase the temperature gradients so that, by the time of ignition, the temperature gradients, and thus the mixing time scales, are all comparable. It is well known that scalar mixing time scales at the molecular level are essentially driven by

the large-scale turbulence integral time scale due to the classical turbulence cascade from large to small scales. The consequence here is that since the initial T  is equal, and the scalar mixing time scale rapidly becomes comparable, temperature gradients, which appear in the denominator of the expression for the mixing time scale are also comparable. Since the importance of molecular transport within the ignition front is controlled by the temperature gradient, this explains the observed insensitivity of the results to the initial length scale. Also, in contrast to the equivalent cases without turbulence fluctuations (cf. Fig. 13), the effects of turbulent straining reducing the scalar mixing time scale imply that passive mixing is important for Cases L = 0.82 and L = 1.64. Finally the multizone model predictions for L = 0.82 and L = 1.64 are shown in Figs. 21a and 21b.

156

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

Fig. 21. Normalized mean heat release rate against time for the DNS and the multizone model started from the DNS initial condition and at 10% pressure rise: (a) Case L = 0.82, (b) Case L = 1.64.

Table 5 Comparison of burn duration (dt) and normalized heat release rate (HRmax /HR0 ) for DNS with multizone model for the parametric study in LTe with turbulence Case

L = 0.41

L = 0.82

L = 1.64

LTe (mm)

1.32

2.64

5.28

DNS dt (ms) MZ-0% dt (ms) MZ-0% dt error (%) MZ-10% dt (ms) MZ-10% dt error (%)

0.64 1.14 78.1 0.72 12.5

0.75 1.12 49.6 0.86 15.2

0.87 1.17 34.6 1.02 17.1

DNS HRmax /HR0 MZ-0% HRmax /HR0 MZ-0% HRmax /HR0 error (%) MZ-10% HRmax /HR0 MZ-10% HRmax /HR0 error (%)

0.34 0.16 54.8 0.29 16.4

0.28 0.16 43.3 0.22 21.4

0.24 0.15 35.2 0.19 22.0

Case L = 0.41 is the baseline case and was already shown in Fig. 8c. Results for the burn duration and maximum heat release rate are reported in Table 5. For the multizone predictions started from the DNS initial condition, the prediction worsens with decreasing length scale. This is due to differences in the passive scalar mixing prior to ignition. Although mixing time scales are comparable throughout the main interval of heat release rate, they are initially very different, leading to differences in the cumulative scalar mixing prior to ignition. However, when the multizone model is started from the 10% pressure rise point, there is no clear trend of larger length scales improving the predictions; on the contrary, there appears to be a weak trend of worsening predictions. This is consistent with the observed fraction of deflagrations and the Damköhler number based on the mixing time scale.

6. Suggestions for practical application of the results Many of the detailed quantities available to the DNS are not available in a practical setting, particularly those relating to the front speed criterion. However, using basic turbulence scaling arguments, it may be possible to estimate whether transport effects are important within ignition fronts. These are indicated by a front speed ratio sig /sL , where sig is the expected speed a spontaneous ignition front. According to Zel’dovich’s [6] theoretical suggestion and as demonstrated in Part I [5], the spontaneous ignition front speed is given by     dτ  −1  sig = |∇T0 | (3) , dT0  where τ is the ignition delay time. Temperature gradients will likely be difficult to measure; however, if

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

a turbulence time scale τt is available, perhaps from CFD calculations, they can be estimated from standard turbulence scaling laws as  −1/2 , |∇T0 | ≈ T  2αcφ−1 τt (4) where cφ ≈ 2 is the ratio of the scalar to turbulence mixing timescales, and α is the diffusivity. These fluctuating gradients should be added to any mean gradients. The quantity dτ/dT0 could be obtained from single-zone calculations, or more accurately from a multizone calculation. Thus if estimates of T  , τt , dτ/dT0 , and sL are available, the ratio sig /sL can be estimated. The front speed ratio indicator should also apply in the boundary layer, but in this case the standard expression for the gradients should be replaced with the appropriate expression in the boundary layer.

7. Conclusions A parametric study of the autoignition of a thermally stratified lean hydrogen–air mixture has been performed. The parametric study was focused on the influence of the amplitudes and length scales of the initial temperature fluctuations and the presence of turbulence. The primary goal was to understand the effect of these parameters on the mode of combustion and on the performance of a model that is applicable to stratified ignition at typical HCCI conditions, the multizone model. The combustion mode was understood by visual comparison of the observed heat release rate fields and by two diagnostic techniques that were developed in Part I of this study [5]. The first diagnostic technique was based on tracking the speed of the advancing combustion wave, and comparison of this with a nominal deflagration speed to determine the significance of molecular transport effects within the ignition front. The second technique was to use the comparison of the mixing time scale with the ignition delay time to determine the influence of passive scalar mixing changing the probability density of temperature in the domain. These diagnostic techniques were applied to the data and used to explain the performance of the multizone model. In all cases the observed performance could be adequately explained using these simple arguments. The parametric study of the RMS temperature fluctuation revealed that this parameter has a strong influence on the observed combustion mode, and on the timing and duration of heat release rate. Higher T  was found to increase the combustion duration and advance the ignition timing. More volumetric, homogeneous combustion was observed for lower T  , whereas higher T  cases showed evidence of combustion in fronts. Using the flame speed diagnostic, it was found that higher T  levels led to a greater prevalence

157

of molecular transport within the fronts (i.e., deflagrations). The performance for the multizone model was linked to the effects of passive scalar dissipation and the prevalence of deflagrations. Very good performance was obtained in cases without significant deflagrations and when passive scalar mixing has been mostly accounted for, which are typical of HCCI conditions. Predictions deteriorated as the prevalence of deflagrations increased, i.e., at higher T  . The parametric study in the initial temperature fluctuation length scale without turbulent velocity fluctuations showed that larger length scales lead to thicker front structures, and a decrease in the importance of molecular transport effects both in deflagrations and through passive scalar mixing. The multizone model showed excellent predictions in the cases with larger length scales. The parametric study of the initial length scales with turbulent velocity fluctuations, where the turbulence time scale was kept constant, showed an increased range of scales in the larger length-scale cases. Turbulence had a strong effect in elongating and folding the initial hot spots. Surprisingly it was found that the importance of deflagrations was roughly independent of the initial length scale. This was explained considering the fact that the turbulence time scale was the same for each case and using the well-known characteristic of turbulence that molecular dissipation time scales, and thus temperature gradients, tend to scale with the turbulence integral time scale. The multizone model predictions for the cases with turbulence were significantly worse than those without turbulence, and furthermore showed a counterintuitive trend that predictions did not improve with increasing length scale, which was explained considering the roughly equal importance of deflagrations occurring in the ignition fronts. The main contribution of this study is the development of the tools needed for understanding the combustion mode in a weakly stratified charge, and understanding the valid regime of the multizone model. However, more information is needed to evaluate whether the chosen conditions are representative of those in a real engine. More information is needed on typical scalar and turbulence fluctuation amplitudes and length scales, both in the core of the charge and in boundary layers which may have very different characteristics. This basic information, along with classic turbulence scaling laws and the diagnostic arguments presented in this study, may go quite far toward predicting the combustion mode in HCCI conditions, and in conditions where a more significant stratification is deliberately introduced into the charge. Advanced optical diagnostic measurements or large eddy simulation could be particularly useful in providing such fundamental characterization.

158

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159

Regarding the multizone model, the results demonstrate, using a high-fidelity simulation approach, that the model performs very well in its intended regime of validity—that of a weakly stratified charge. However, molecular transport effects can cause the model to fail for more stratified charges. Whether or not these are significant in an engine is not assessed in the present work, rather, tools are provided to address that question. These tools can be used to understand why this model has been successful in applications at HCCI conditions to date and can guide future applications. Also, a comment can be made about the way in which mixing degrades the results so that these effects could be identified if they were occurring. In all cases considered, both passive scalar mixing and transport within ignition fronts tend to (1) increase the peak heat release rate and (2) decrease the burn duration relative to the multizone prediction. In order to distinguish the effects of passive mixing versus transport within ignition fronts it is worth noting that transport within ignition fronts will act to advance the timing of the maximum heat release rate, while passive mixing will tend to retard the timing of maximum heat release rate. A very encouraging result of the present work is that for most conditions that were simulated, the main molecular transport effects were manifested in the change of the scalar PDF in time—i.e., turbulent mixing effects. In turbulent flows, this effect may be curtailed in certain cases due to the expected cascade from larger scales. To assess this question, careful thought needs to be given to the source of the scalar fluctuations and whether that is of a continuous nature, as in boundary layer entrainment of cooler fluid, or decaying, as in mixing of fresh charge with residual gases. If it is expected that the PDF changes with time, and the time scale is comparable to the duration of combustion after the multizone model is started, then turbulent mixing needs to be accounted for. In some circumstances, turbulent mixing may also need to be accounted for even if the scalar PDF is not changing in time, since a strong enough turbulent mixing occurring on a timescale comparable to the range of ignition delay times could cause diffusion from burned to unburned zones even without significant propagation by deflagrations. Simple modifications to the multizone model could potentially account for these effects [12]. Also, unsteady laminar flamelet [13,14], conditional moment closure [15], and the transported PDF [16,17] approach all have mechanisms that account reasonably well for these effects. The presence of normal premixed deflagrations on the other hand may be poorly accounted for by these approaches, and since the deflagrations exist in an igniting medium, also by conventional premixed flame modeling approaches such as the G-equation [13].

It would be very interesting to assess the performance of these alternative approaches against the present DNS data.

8. Acknowledgments Sandia National Laboratories (SNL) is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under Contract DE-AC04-94AL85000. The work at SNL was supported by the Division of Chemical Sciences, Geosciences and Biosciences, the Office of Basic Energy Sciences, the U.S. Department of Energy. Calculations were performed at SNL and at the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC03-76SF00098. The High Performance Computing and Networking department at SNL provided access to a 256-processor Infiniband testbed. The authors acknowledge fruitful discussions with Drs. John Hewson, Hong Im, John Dec, and Magnus Sjöberg.

References [1] K. Epping, S. Aceves, R. Bechtold, J. Dec, SAE Technical Paper 2002-01-1923, 2002. [2] S.M. Aceves, D.L. Flowers, C.K. Westbrook, J.R. Smith, W. Pitz, R. Dibble, M. Christensen, B. Johansson, SAE Paper 2000-01-0327, 2000. [3] S.M. Aceves, D.L. Flowers, J. Martinez-Frias, J.R. Smith, C.K. Westbrook, W. Pitz, R. Dibble, J.F. Wright, W.C. Akinyemi, R.P. Hessel, SAE Paper 2001-011027, 2001. [4] S.M. Aceves, D.L. Flowers, F. Espinosa-Loza, J. Martinez-Frias, R.W. Dibble, M. Christensen, B. Johansson, R.P. Hessel, SAE Paper 2002-01-2869, 2002. [5] J.H. Chen, E.R. Hawkes, R. Sankaran, S.D. Mason, H.G. Im, Combust. Flame 145 (2006) 128–144. [6] Ya.B. Zel’dovich, Combust. Flame 39 (1980) 211–214. [7] R. Sankaran, H.G. Im, E.R. Hawkes, J.H. Chen, Proc. Combust. Inst. 30 (2005) 875–882. [8] M.A. Mueller, T.J. Kim, R.A. Yetter, F.L. Dryer, Int. J. Chem. Kinet. 31 (1999) 113. [9] J.O. Hinze, Turbulence, McGraw–Hill, New York, 1975. [10] M. Sjöberg, J.E. Dec, N.P. Cernansky, SAE Technical Paper 2003-01-0113, 2003. [11] T. Echekki, J.H. Chen, Combust. Flame 118 (1999) 308–311. [12] D. Flowers, S. Aceves, J. Martinez-Frias, R. Hessel, R. Dibble, SAE Paper 20003-01-1821, 2000. [13] N. Peters, Turbulent Combustion, Cambridge Univ. Press, New York, 2000.

E.R. Hawkes et al. / Combustion and Flame 145 (2006) 145–159 [14] G. Blanquart, H. Pitsch, Proc. Combust. Inst. 30 (2005) 2745–2753. [15] A.Y. Klimenko, R.W. Bilger, Prog. Energy Combust. Sci. 25 (1999) 595–687.

159

[16] S.B. Pope, Prog. Energy Combust. Sci. 11 (1985) 119– 192. [17] Y.Z. Zhang, E.H. Kung, D.C. Haworth, Proc. Combust. Inst. 30 (2005) 2763–2771.