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

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

Direct numerical simulation of ignition front propagation in a constant volume with temperature inhomogeneities I. Fundamental analysis and diagnostics Jacqueline H. Chen a , Evatt R. Hawkes a,∗ , Ramanan Sankaran a , Scott D. Mason b , Hong G. Im c a Reacting Flow Research Department, Combustion Research Facility, Sandia National Laboratories, P.O. Box 969 MS 9051,

Livermore, CA 94551-0969, USA b Lockheed Martin Corporation, Sunnyvale, CA 94089, USA c Department of Mechanical Engineering, University of Michigan, Ann Arbor, MI 48109-2125, 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 with a view to providing better understanding and modeling of combustion processes in homogeneous charge compression-ignition engines. Numerical diagnostics are developed to analyze the mode of combustion and the dependence of overall ignition progress on initial mixture conditions. The roles of dissipation of heat and mass are divided conceptually into transport within ignition fronts and passive scalar dissipation, which modifies the statistics of the preignition temperature field. Transport within ignition fronts is analyzed by monitoring the propagation speed of ignition fronts using the displacement speed of a scalar that tracks the location of maximum heat release rate. The prevalence of deflagrative versus spontaneous ignition front propagation is found to depend on the local temperature gradient, and may be identified by the ratio of the instantaneous front speed to the laminar deflagration speed. The significance of passive scalar mixing is examined using a mixing timescale based on enthalpy fluctuations. Finally, the predictions of the multizone modeling strategy are compared with the DNS, and the results are explained using the diagnostics developed.  2005 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Ignition; Direct numerical simulation; HCCI; Multizone model

1. Introduction

* Corresponding author. Fax: +1 925 294 2595.

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

Homogeneous charge compression-ignition (HCCI) engines are considered a viable concept as an alternative to diesel engines [1]. In contrast to conventional compression-ignition (CI) engines, HCCI engines exploit a lean intake charge that is well mixed

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

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

prior to combustion. HCCI engines may thus provide efficiency gains over spark-ignition engines, and lower emissions as compared with CI engines. At present, the combustion mode in HCCI engines is not well understood, as both volumetric and frontlike combustion modes can occur [2]. A major challenge posed by this method of combustion is to control the heat release rate, and in particular, to provide a means to spread it out over several crank angle degrees, suppressing the occurrence of rapid rate of pressure rise. One possible control strategy is to introduce inhomogeneities in the temperature or mixture composition in order to produce the desired heat release rate [1,3]. Moreover, incomplete turbulent mixing and temperature stratification between the bulk gases and the cylinder wall lead to spatial nonuniformities that also contribute to a range of combustion modes that are distinct from homogeneous autoignition. Therefore, a better understanding of autoignition of an inhomogeneous charge at constant volume is crucial in the development of control strategies for HCCI engines. The overall objectives of the present study are to understand the influence of temperature inhomogeneities on characteristics of ignition and subsequent combustion in an HCCI-like environment and to provide information relevant to the modeling of HCCI combustion. The evolution of autoignition from an initial spectrum of “hot spots” at high pressure is simulated using two-dimensional direct numerical simulations (DNS) with a detailed hydrogen/air reaction mechanism. Although turbulence is inherently threedimensional, the prescribed two-dimensional random temperature distribution can still provide useful information regarding the chemical induction and front propagation behavior in response to different statistics of the imposed scalar and velocity fields. Due to the computational expense of resolving extremely thin ignition fronts at high pressure throughout their temporal evolution, the chemical kinetics considered are necessarily oversimplified to a hydrogen/air system. Hydrocarbon fuels, which build on hydrogen as an important subset of the mechanism, will be considered in future studies. One of the primary objectives in this study is to develop systematic and rational methods to identify and analyze various regimes of ignition and subsequent combustion processes. The theoretical basis stems from an earlier study by Zel’dovich [4], who described ignition of a nonuniform mixture in several distinct regimes of propagation, including deflagration, spontaneous ignition front propagation, and detonation. He further showed that the spontaneous ignition front speed is inversely proportional to the initial temperature gradient, suggesting a characterization of ignition regimes based on the speed of the igni-

129

tion front. For the parametric conditions applicable to HCCI engines, formation of a detonation wave is unlikely and thus the first two regimes are considered relevant in the present study. The concept has been applied to identifying ignition regimes in various onedimensional configurations [5–8]. The present work extends the idea into more realistic two-dimensional ignition problems, where details of multidimensional aspects including wrinkling and mutual interaction of fronts are fully accounted for. The second main objective is to provide information relevant to the modeling of HCCI. In the case of purely homogeneous combustion, the modeling problem reduces to simple homogeneous ignition calculations. However, stratification of the charge, whether in terms of temperature or mixture composition, introduces a range of ignition delay times into the cylinder, which must be accounted for in a model. At low levels of stratification or long stratification length scales, typical of HCCI engine conditions, molecular mixing effects can be neglected [4,9]. Presuming that length scales are not so long so that compressible flow effects are important, a reasonable approximation is then that different locations in the cylinder interact only by pressure-work. These assumptions are employed in a modeling strategy referred to as the multizone model, which has been advanced by Aceves et al. [10–12]. The model has been demonstrated to give good predictions of engine pressure traces and heat release rates at low computational cost in typical HCCI engine conditions. There is a need to understand why these predictions accurately reproduce the pressure and heat release rate traces and to develop techniques for estimating a priori the region of applicability. Using the DNS as a numerical experiment, we evaluate the performance of the model, and provide understanding of its limitations. In order to develop methods of determining the model’s validity, we will deliberately apply it in some regions outside the valid regime for which it was intended. Also, the multizone model serves as a useful foil to help validate the tools developed to characterize the diffusive transport processes precisely because the model does not include these processes and hence any discrepancies between the model and the DNS give an independent measure of their importance. Our previous DNS study Sankaran et al. [13] investigated the characteristics of turbulent ignition by considering three initial temperature fields characterized by different skewness of the fluctuation. In this study, a heuristic temperature gradient cutoff condition was used to quantify the relative dominance between spontaneous ignition and deflagration regimes. A criterion to predict the relative dominance of the two ignition regimes was also proposed.

130

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

The present study is composed of two parts. This paper, Part I, extends the previous work [13] to explore the fundamental issues with an emphasis placed on front propagation aspects of the phenomena. The main objective is to develop numerical diagnostics to identify and analyze various regimes of ignition and subsequent combustion processes. We also present comparisons of the DNS results with the multizone model, and identify the reasons for discrepancies. Using the tools developed in this paper, in Part II [14] a parametric study is performed on the influence of three key controlling parameters: the initial stratification amplitudes, the stratification length scales, and the turbulence intensity. Comparisons with the multizone model are presented for each case. In particular in the present paper, the utility of a scalar isosurface to track the propagation speed and behavior of the ignition fronts is assessed. The role of molecular diffusive effects including heat conduction and mass diffusion on ignition propagation is identified. Molecular diffusive effects are divided conceptually into passive scalar mixing and diffusion within propagating ignition fronts, with the latter controlling the transition between the regimes identified by Zel’dovich. Ignition front propagation evolution during constant volume combustion is monitored by evaluating the surface mean of the ratio of the local front speed to the local deflagration speed. This local measure is also used to distinguish front propagation in the different regimes. The results are compared with a volumetric measure based on a critical temperature gradient [13]. Predictions of the multizone model are presented and the reasons for discrepancies are explained.

2. Numerical implementation and initial conditions The present study is focused primarily on the effect of temperature inhomogeneities on autoignition in the presence of compression heating resulting from combustion heat release rate due to ignition front propagation. In an HCCI engine, these effects are important near top dead center. It is assumed that the isentropic compression stage, where turbulent mixing of inert mixtures prevails, has already occurred. Therefore, the present simulations focus on autoignition of a lean hydrogen–air mixture in a closed volume in the presence of an inhomogeneous temperature field. A uniform hydrogen/air mixture with a fuel-air equivalence ratio of 0.1 is chosen. This mixture has been chosen such that the pressure rise is sufficient for studying ignition under compression heating conditions, but is insufficient for shock wave formation. It is known that as the volumetric heating

increases there exists a potential for the coalescence of ignition fronts to form shock waves [6,7], which the present numerical grid cannot resolve. In all cases presented, the mean initial temperature is 1070 K, and the uniform initial pressure is 41 atm. These conditions correspond to a compression ratio of 15:1 starting from 1 atm and 400 K. At this pressure and temperature an insignificant amount of heat release rate has occurred for the hydrogen/air mixture. The primary source of heat release rate is through the thermal dissociation of hydrogen peroxide which is negligible at this pressure until the temperature exceeds 1100 K.1 Hence, it is valid for the present H2 /air mixture to decouple the initial isentropic compression stage from the combustion-induced compression heating stage. The ignition behavior is first analyzed by a simple one-dimensional model in order to provide a reference case to help interpret the two-dimensional simulation results. The 1D constant-volume cases were initialized from a sinusoidal temperature profile. Boundary conditions are periodic, such that heat release rate in the computational domain leads to a pressure rise and compression heating of the reactants as in the 2D calculations. Parametric studies are conducted by varying the length scale and amplitude of the temperature fluctuation. In the latter part of the paper, two-dimensional simulation results, previously studied as Case A in Ref. [13], are analyzed further. In the two-dimensional turbulent ignition case, the turbulent flow field is prescribed by an initial turbulent kinetic energy spectrum function by Passot–Pouquet [15], 32 E(k) = 3




      2 u 2 k 4 k 2 , exp −2 π ke ke ke

(1)

where k is the wave number magnitude, ke is the most energetic wavenumber, and u is the RMS velocity fluctuation. A random temperature field is 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 shown in Fig. 1. Turbulent mixing modulates the scalar gradients imposed by the initial inhomogeneities, affecting the local dissipation rates of key radicals and thermal energy. Therefore, the specification of the initial tem1 However, this would not be the case for a multistage ignition of a hydrocarbon fuel such as n-heptane where low-temperature reactions involving the thermal dissociation of ketoheptylhydroperoxide generate significant heat release rate at temperatures greater than approximately 750 K at this pressure.

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

131

reactive species and 21 reversible reactions. The domain is 4.1 mm in each of the two spatial directions. A uniform grid spacing of 4.3 µm was required to resolve the ignition fronts at this pressure.

3. Ignition front tracking and species budget

Fig. 1. Initial temperature spectrum for 2D DNS. Temperatures vary between 1106 K (blue) and 1203 K (red).

perature and turbulence fields affects both the ignition timing and the subsequent heat release rate. The initial autocorrelation integral scale of the velocity fluctuations is 0.34 mm, and the most energetic length scale is 1.0 mm. This length scale corresponds to the wavelength in the spectrum where the kinetic energy is maximum. The velocity RMS is 0.5 m/s. These parameters were selected so that the turbulence integral time scale based on the most energetic length scale is comparable to the homogeneous ignition delay of the mixture at the same mean temperature and pressure, thus promoting strong interaction between the turbulent mixing and ignition chemistry. The homogeneous ignition delay is τ0 = 2.9 ms. The temperature field has an RMS fluctuation of 15.0 K, the autocorrelation length scale is 0.15 mm, and the most energetic length scale is 1.32 mm. These parameters are comparable to those in an engine operating at low load, with the exception that the length scale of the exothermic regions and the velocity fluctuations are approximately five times smaller than observed in an engine [2]. The resulting turbulence integral time scale and ignition delay are comparable to those timescales in an engine. The conservation equations for a compressible reacting flow are solved in the DNS using an explicit eighth-order accurate finite differencing scheme [16] and a fourth-order accurate Runge–Kutta scheme for time advancement [17]. The boundary conditions are periodic and lead to compression heating of the reactants as the mixture ignites. The mixture-averaged thermal conductivity is temperature-dependent [18], and the individual species specific heats are obtained using the Chemkin thermodynamic database [19]. The diffusion coefficients are obtained by prescription of Lewis numbers for individual species, fitted from mixture-averaged transport coefficients [19]. The pressure-dependent hydrogen–air chemistry is based on a detailed chemical mechanism [20] with 8

According to Zel’dovich [4], the speed of the advancing combustion wave, given an inhomogeneous initial condition for temperature, is related to the gradients of temperature and to the mode of combustion. Zel’dovich [4] identifies the following main regimes in order of decreasing speed and increasing initial temperature gradient, which are also discussed in [6]: (i) a nearly instantaneous thermal explosion, (ii) an autoignitive spontaneous ignition front propagating at speeds greater than the normal detonation rate, (iii) a developing detonation, (iv) a subsonic spontaneous ignition front propagation, in which molecular diffusion effects ahead of the front are unimportant, and (v) a normal deflagration, in which mass and heat transport at the molecular level control the propagation. The focus of the present study is on the transition between the normal deflagration regime and the subsonic spontaneous ignition front propagation regime. In the subsonic spontaneous ignition front regime, the front speed sig may be related to the gradient of the ignition delay time τ [4] as follows: sig =

1 . |∇τ |

(2)

Note that, in the constant volume case, the reactivity of a mixture depends not only on its initial state, but also on the amount of compression heating it receives during the progress of ignition. Therefore, τ depends not only on the local mixture condition but also on the distribution of mixture conditions in the volume which affect the amount of pressure work. When the composition is uniform, the gradients of ignition delay are solely due to temperature variations, and the equation may be expanded as 1 , sig =  dτ   |∇T0 | dT

(3)

0

where T0 is the initial temperature. For low initial temperature gradients, it is evident that the spontaneous ignition front speed becomes large. In this situation there is insufficient time for molecular transport toward the reactants to be significant, and the front speed is well approximated by Eq. (3). However, as the temperature gradient is increased, the front speed slows down to the point where the molecular transport to the fresh gases becomes significant. The observed speed of the propagating front is then

132

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

increased relative to the prediction of Eq. (3). For sufficiently high temperature gradients, the propagation is controlled by a reactive–diffusive balance and is similar to a normal laminar premixed flame, although the flame is propagating into an unsteady autoigniting mixture. These theoretical predictions underline the importance of both the local temperature gradient, and the front speed as indicators of the combustion regime. In this study, we explore the utility of these quantities as delineators of the regime. To obtain the speed of the advancing combustion wave, a surface tracking technique is employed. While it is clear that in these conditions reaction is not always confined to a thin surface, it does provide a very convenient means for analyzing the DNS results. If an ignition delay time is defined for each fluid parcel as the time taken for a scalar φ to reach a critical level φc , then the speed sd of the ignition front may be identified with the displacement speed [21] given by  ±Dφ/Dt  , sd = (4) |∇φ| φ=φc where D/Dt is the material derivative, and the sign convention depends on whether the scalar increases (+) or decreases (−) as the critical value is approached. We have developed analytic front speed expressions for a temperature isosurface and a scalar isosurface based on the fuel mass fraction. While both definitions give very similar conclusions, the temperature evolution equation contains multiple terms representing different physical effects. Therefore, in this study only results from the mass fraction isosurface are presented. A density-weighted speed can be defined as sd∗ = ρsd /ρ0 , where the ρ0 is a representative density of the reactants. This is calculated from the local enthalpy and unburnt mixture condition assuming pressure and enthalpy do not vary as the front passes. Statistics were extracted from the YH2 = 8.5 × 10−4 isocontour, which corresponds approximately to the location of maximum heat release rate through most of the simulation. A nominal measure of the front speed of a diffusion controlled deflagration is required in order to compare with the speed of the propagating front. Here, we employ the steady premixed flame speed at the local conditions. It is recognized that because of the highly transient conditions with temperature gradients and autoigniting reactants, this is simply an approximation, however it will be later shown that it provides a good estimate of the transition between regimes that is consistent with several other measures. Laminar flame speeds are obtained with the PREMIX code [22], and are tabulated as a function of the reactants enthalpy and pressure.

In the following sections, the budget terms in the species continuity equation for hydrogen are examined to determine the relative importance of diffusion to reaction in determining the front propagation regime. The evolution of the hydrogen mass fraction, YH2 , is given by DYH2 ∂YH2 1 ∂ ω˙ = ρDH2 + , Dt ρ ∂xi ∂xi ρ

(5)

where the first term on the right-hand side represents diffusion, and the second term represents the net reaction rate of hydrogen. DH2 is the diffusivity of hydrogen, ω˙ is the reaction rate, and ρ is the density.

4. Multizone model The multizone model [10–12] has been demonstrated to accurately predict pressure traces and heat release rates in engine-derived HCCI conditions, with a very low computational overhead. It is an approximation that is applicable in the subsonic spontaneous ignition front propagation regime. It ignores effects associated with molecular diffusion and with the finite speed of sound. An initial spatial distribution of mixture states obtained from a prior nonreacting computational fluid dynamics calculation which fully includes turbulent mixing is approximated by a finite number of zones such that a joint probability density function (PDF) of temperature and species mass fractions approximately matches that of the original distribution. These zones are then treated as homogeneous reactors that are coupled together only by pressure, which is assumed to equilibrate instantaneously. In HCCI engine conditions, the pressure coupling is expected to be the most significant interaction between different parts of the cylinder. In the constant volume case, the evolution equations for the species mass fraction and temperature in each zone may be summarized as follows, ωi,j dYi,j = , dt ρj

  Nz dTj 1 k=1 Ak Bk Aj + , =  z dt ρj cp,j 1− N k=1 Bk

(6)

(7)

where the subscript j refers to the zone to which the equations apply, the subscript i refers to the chemical species, the subscript k refers to any zone, Nz is the total number of zones, Yi,j is the mass fraction, ωi,j is the reaction rate, ρj is the density, Tj is the temperature, and cp,j is the specific heat at constant pressure. For each zone the quantities Aj and Bj are defined as follows,

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

Aj =

Ns 

Nz 

i=1 k=1

dYi,k Vk ρk Ri Tk V dt

  Vj 1 1− , Bj = V γj



133

dYi,j − ρj hi,j , dt (8) (9)

where Ns is the number of chemical species, Vj is the  z volume of a zone, V = N k=1 Vk is the constant total volume of all zones, γj is the ratio of specific heats, hi,j is the species enthalpy, and Ri is the species gas constant. The uniform pressure is determined from the equation of state applied to all zones, and the individual zone volumes are determined from the equation of state applied to each zone. The multizone model is initialized from the DNS solution at a specified time by sampling a number of points. In the 1D cases, each DNS grid point is treated as a zone, whereas in the 2D case a reduced number is employed (57,600). The equations are integrated with the identical numerical time integration method, reaction mechanism and thermodynamic approximations as in the DNS. It is noted here that fewer zones, typically around 10–20, are found to be sufficient to satisfactorily reproduce experimental pressure traces. Thus, multizone calculations have a low computational cost which allows the treatment of very detailed chemical kinetic models. Here, a very large number of zones was used in all cases simply so that results are ensured to be independent of the number of zones and to allow an exclusive focus on the physical issues rather than those of numerical convergence. Yet the cost of the multizone model calculations was negligible when compared to the complete DNS.

5. One-dimensional test cases To understand the combustion process in the presence of temperature inhomogeneities, 1D test calculations were performed. These test cases provide a reference case and a basis for understanding the more complicated 2D turbulent DNS results. In a first parametric study, the RMS temperature fluctuation is fixed at 15.0 K, matching the twodimensional Case A. To help determine the role of temperature gradient and the wavelength of the initial fluctuations, cases were run with different wavelengths (4.1, 1.5, 0.75, 0.56, and 0.38 mm) representing the spectrum of length scales observed in the 2D simulations. These cases will be referred to henceforth by these wavelengths. In a second parametric study, the wavelength is fixed at 3.0 mm, and the RMS temperature fluctuation is varied. RMS temperature fluctuations of 3.75, 7.5, 15.0, and 30.0 K were employed. These cases will be henceforth referred to by these amplitudes. Comparison of the results of the two

Fig. 2. (a) Temperature versus distance for a sequence of times; (b) heat release rate normalized by the maximum heat release rate of the homogeneous ignition at the mean temperature (HR0 ) versus distance. The equally spaced time sequence is numbered from 1 to 10, starting at 1 ms with an increment of 0.25 ms. Table 1 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

parametric studies will help delineate between molecular diffusion effects occurring in ignition fronts, and the effects of passive scalar dissipation. The overall progression and character of the ignition and combustion process in this 1D configuration is illustrated in Fig. 2. Fig. 2a shows the temperature profile and Fig. 2b shows the heat release rate at several instants in time. All results presented in the present paper are normalized by their corresponding reference values for the maximum heat release rate, HR0 , and ignition delay, τ0 , corresponding to homogeneous ignition at the mean temperature, 1070 K and initial pressure of 41 atm summarized in Table 1. The time increments in these figures are 0.250 ms, or approximately 0.1 times the homogeneous ignition delay τ0 . Combustion in this configuration proceeds as follows. Ignition occurs first at the location of highest temperature, in the middle of the domain, moderated by the effects of scalar dissipation in the ignition kernel. Subsequently, a propagating combustion wave emanates from this location and travels to the left and right of the domain. As fronts propagate, the remaining charge is heated by compression, accelerating the ignition of the end-gas. Finally the end-gas

134

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

Fig. 3. (a) Front speed (sd∗ ) and (b) |∇T | versus time for wavelengths ranging from 0.56 to 4.10 mm in the 1D test cases.

is consumed as the front arrives at the boundaries of the domain. The presence of a propagating front does not imply that this is a normal deflagration. The nature of the front and the role of molecular diffusion need to be identified. Fig. 3a shows the extracted front speed sd∗ versus time for the various wavelengths. Also shown is a representative deflagration speed sL that has been determined based on a freely propagating premixed flame [22] for the enthalpy and pressure conditions at the front surface. Fig. 3b shows the corresponding evolution in the magnitude of the temperature gradient at the surface location, |∇T |. Comparison of the two figures indicates that the “U-shape” behavior in the front speed is due to the influence of the gradient term in Eq. (3). Both the initial ignition and final consumption points have a zero initial temperature gradient and the speed becomes unbounded.2 Comparison of the different wavelength cases shows that the initial ignition time is delayed for shorter wavelengths due to diffusive losses from the ignition kernel. On the other hand, the final consumption occurs earlier for shorter wavelengths, due to diffusive gains in the end gas parcel. In the longer wavelength limit where the spontaneous ignition front speed is much larger than the deflagration speed (but is still subsonic), diffusion is insignificant and the ignition delay and combustion duration are governed simply by constant-volume homogeneous ignition, moderated by the effects of compression heating.

2 It is worth noting that while the speed does become un-

bounded, this does not mean that consumption of reactants is unbounded in the present case. In fact, the description of combustion as a front breaks down at the early stage of ignition, when combustion occurs in a kernel, and at the final stage, when the fronts annihilate. Unbounded front speeds are theoretically possible in the case of a stratification occurring over an unbounded length scale. Practically, of course, this would never occur.

In the other limit of small wavelengths, the initial temperature fluctuations dissipate before significant reaction occurs and the entire domain again ignites homogeneously according to the mean temperature. Therefore, spontaneous ignition front propagation is expected for both very long and short wavelengths. Between the two limits, a situation may exist where the local temperature gradients at the ignition time are sufficient to develop a deflagration wave whose propagation speed is higher than the spontaneous ignition front speed at this local condition. Therefore, for an intermediate hot-spot size, the ignition process may undergo transition from a spontaneous propagation to deflagration, and then back to spontaneous propagation again. It is important to emphasize that diffusive effects can be conceptually divided into two different phenomena: (1) passive scalar dissipation—essentially the modification of the initial distribution of temperature within the domain, independently of reaction, and (2) diffusive effects occurring within the front, as in a normal premixed deflagration. The two effects are governed by different nondimensional parameters. Passive scalar dissipation may be important if the scalar dissipation timescale is comparable to the ignition delay, introducing a Damköhler number Da = τmix /τ0 . On the other hand, the effects of diffusion in the front become important when sig /sL ≈ O(1). In the present 1D parametric study, the RMS temperature fluctuation is held fixed while the wavelength is varied. This results in a change in both the Damköhler number and the front speed ratio. Due to the range of wavelengths used in the 1D cases, both effects were observable. The wave-propagation stage of the combustion process occurs near the bottom of the U-curve in Fig. 3a. For the largest wavelength of 4.1 mm, the temperature gradient never becomes large enough to develop a deflagration wave, such that the entire process occurs in the spontaneous propagation regime at a substantially higher speed. As the wavelength

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

135

Fig. 5. Mean temperature versus time for 1D test cases and multizone model (MZ 0%). MZ 0% refers to the multizone model initialized by the initial DNS field at 0% pressure rise.

Fig. 4. (a) YH2 budget terms and (b) temperature gradient in fronts, 1D test cases with 4.1 and 0.75 mm wavelengths. Dashed line: reaction. Solid line: diffusion. Dotted line: temperature gradient.

is decreased, the temperature gradient increases and the deflagration regime starts to occur, as evidenced by the existence of the minimum propagation speed. Note that further decrease of the wavelength from 1.5 to 0.75 mm does not result in a decreased front speed. Even further decrease in the wavelength (0.56 mm) approaches the small-wavelength limit behavior, where passive scalar dissipation dominates, and the speed increases again. The lower limit of the front speed corresponds to the critical initial temperature gradient for which diffusion becomes important, and the minimum front speed can be identified with a deflagration speed. This speed is approximately 50 cm/s. The transition evidently occurs near the 1.5mm wavelength, which corresponds approximately to the most energetic wavelength of the initial temperature spectrum in the DNS, indicating that both deflagration and spontaneous ignition front propagation regimes are expected in the 2D simulations. Furthermore, it is observed that the heuristic temperature gradient cut-off of 1800 K/mm used in Sankaran et al. [13] is consistent with this transition. Note that the speed independently obtained from PREMIX agrees with the lower limit. The quantitative agreement suggests that it can serve as a nominal measure of the reference deflagration speed to be used in the analysis of the 2D results, even under highly transient conditions. The above observation is further supported by examining the structure of the fronts in the different ignition regimes. Fig. 4 shows the reaction and diffusion budgets for the H2 mass fraction given by Eq. (5) along with the temperature gradient for the combus-

tion wave at the time of minimum propagation speed. The 4.1-mm case at this point is propagating as a spontaneous ignition front, and diffusion is nearly insignificant. For the 0.75-mm case, for which the front speed is limited from below by the effects of diffusion, a greater influence of diffusion is observed. These results suggest that the relative magnitude of diffusion is an alternative indicator to delineate combustion in the deflagrative and spontaneous ignition regimes.

6. Comparison of 1D cases with multizone model To understand the roles of passive scalar dissipation and in-flame diffusion, we have performed a number of the multizone model calculations initialized from the 1D DNS cases presented in Section 5. A sufficiently large number of zones has been used so that results are independent of the number of zones. The comparison between the 1D DNS results in the previous section and the multizone model results will thus reveal the effect of molecular mixing and dissipation during the ignition event. Fig. 5 shows the mean temperature in the domain for various wavelength cases of the 1D DNS and the prediction from the multizone model initialized from the initial field of the 1D DNS (namely 0% pressure rise). Since diffusive transport between the zones is not considered, the results of the multizone model for all wavelengths are identical because the initial temperature PDF is identical. It is clearly seen that the DNS results for longer wavelengths converge to the multizone model predictions. As the wavelength is decreased, however, the DNS results deviate significantly from the multizone model results because the enhanced molecular mixing of the scalar field decreases the burn duration and increases the maximum rate of heat release rate. It is interesting to note that for the shortest wavelength cases (0.38 and 0.56 mm),

136

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

Fig. 6. Mean temperature versus time for 1.50- and 0.56-mm 1D test cases and multizone model. MZ 0% and MZ 10% represent the multizone models initialized by the DNS field at 0% and 10% pressure rise, respectively.

although the front speeds are larger than the deflagration speed, they still do not agree with the multizone prediction. These cases exemplify the significant role of passive scalar mixing in the induction phase of ignition, thereby strongly modifying the mean temperature evolution. Therefore, it is evident that the multizone model can adequately predict the ignition of stratified mixtures only if the length scales of the stratification are sufficiently large. To better understand the role of passive scalar dissipation, alternative multizone model calculations were made by using the initial condition of the 1D DNS solutions at the point of 10% of the final pressure rise (called MZ 10%). Fig. 6 shows such results for two different wavelength cases, 0.56 and 1.50 mm, compared with the 1D DNS results. The MZ 0% result shown in Fig. 5 is also overlaid. For both wavelengths, a much better agreement between the DNS and multizone model results is found, suggesting that a substantial portion of passive scalar dissipation occurs during the earlier phase prior to the 10% pressure-rise point. Comparing with the MZ 0% result, it is also clearly seen that the dissipation effect is much stronger for shorter wavelengths.

Fig. 7. RMS temperature fluctuation T  versus time for 1D cases and multizone model (MZ 0%).

Fig. 7 shows the RMS temperature fluctuation T  as a function of time for various wavelength cases, along with the multizone model. This further demonstrates that the short wavelength cases exhibit significant dissipation prior to ignition. The extent of this effect is represented by the Damköhler number defined earlier. To unravel the distinct effects of the in-flame transport in contrast to the passive scalar dissipation, a second parametric study is conducted in which the length scale is fixed at 3.0 mm but the amplitude of temperature fluctuations is varied. The choice of the 3.0 mm wavelength is made because the passive scalar dissipation timescale is long compared with the ignition delay. Therefore, any discrepancies with the multizone model may be attributed to transport within the ignition front. Specifically, Da = τmix /τ0 ≈ 8, where the mixing time scale is defined as τmix =

h 2 , 2α|∇h|2

(10)

in which h is the mixture enthalpy, h is the RMS enthalpy fluctuation, and α is the thermal conductivity. Enthalpy is employed rather than temperature because it is less sensitive to changes during reaction. Discounting pressure work and differential diffusion effects, enthalpy would be a conserved scalar. Note that changes in the amplitude of temperature fluctuation, T  , have a negligible effect on the mixing timescale since its contribution in the numerator and denominator cancels. Fig. 8 shows the comparison of the multizone model and DNS results for cases with T  = 7.5, 30, and 60 K. There is near-perfect agreement in the 7.5 K case, but the agreement degrades as the amplitude of temperature fluctuation increases. In this situation, where the passive scalar dissipation is not important, the spontaneous front propagation speed is slow enough due to the increase in the temperature gradient (see Eq. (4)), such that ignition occurs in the

Fig. 8. Mean temperature versus time for rms temperature fluctuation between 7.5 and 60 K in the 1D test cases and multizone model (MZ).

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

137

deflagration regime where heat and mass transport to the front becomes important. In summary, the comparison of the 1D and the multizone predictions reveals that both passive scalar dissipation and molecular diffusion in deflagration fronts may play an important role in constant volume ignition of an inhomogeneous mixture. Moreover, the new diagnostics, based on the front speed, mixing timescales and species evolution budgets, provide a consistent delineation between regimes.

7. 2D DNS results In this section, the utility of the ignition front speed and species budget to delineate between deflagrative versus spontaneous ignition front propagation is further assessed using 2D DNS data of constant volume ignition in the presence of initial temperature inhomogeneities. The simulation data for Case A in our previous study [13] is used. Although details of the ignition processes can be found in [13], they are briefly summarized here for completeness. After an initial ignition delay, which is governed by the competition between reaction and diffusion in several candidate kernels, the strongest ignition kernels ignite and form combustion fronts. The thickness of these fronts varies significantly and determines their propagation speed and the importance of diffusive transport. As fronts propagate, the remaining charge is heated by compression, accelerating the ignition of later-igniting kernels. Near the peak heat release rate, fronts begin to merge and annihilate, and eventually the end-gas is consumed in a more volumetric process. Instantaneous solution fields are analyzed to investigate the characteristics of the front propagation. The 1D test cases showed that the front speed in the deflagration regime is very close to the laminar flame speed at the corresponding local enthalpy and pressure conditions. Therefore, we define a cut-off criterion, sd∗ = 1.1sL , to be a reasonable boundary to distinguish between the two ignition regimes. Fig. 9 shows the heat release rate (color) and YH2 = 8.5 × 10−4 (white/black lines) isocontours at 0.83τ0 when ignition kernels are forming and ignition fronts have begun to propagate. This instantaneous field is representative of the events during the interval of major heat release rate. The figure shows that this isocontour selection provides good tracking of the maximum heat release rate location. Similar results are obtained from the inception of front propagation through the burnout of the end-gas. Along the YH2 isocontour lines, the black and white segments denote the deflagration (sd∗
1.1sL ) and spontaneous ignition regime

Fig. 9. Heat release rate isocontours (rainbow color scale— red denotes maximum heat release rate), YH2 surface isocontour (white lines: sd∗ > 1.1sL , black lines: sd∗
1.1sL ), and location of representative front structures. (A) Spontaneous ignition front; (B,C) deflagration fronts.

Fig. 10. Ignition front length (dashed line) and volumetric mean normalized heat release rate (solid line) versus time for 2D DNS.

(sd∗ > 1.1sL ), respectively, thereby indicating the relative importance of the two regimes at a given instant in time. The evolution of the ignition front length and the volumetric mean heat release rate (normalized by the maximum heat release rate of the homogeneous ignition at the mean temperature) is presented in Fig. 10. The ignition front length is defined simply as the length of the hydrogen mass fraction isoline representing approximately the location of maximum heat release rate. After the initial ignition delay, both the heat release rate and front length begin to rise, reach the maximum at approximately the same time, and diminish as remaining parcels of end-gas are consumed. While it should not be inferred that heat release rate occurs only on an isolated surface, the two quantities do show a qualitative match, validating the use of this isosurface to study the rate of combustion. Note that

138

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

Fig. 12. Density weighted ignition front speed (sd∗ , thick lines) and |∇T | (thin lines) for 2D DNS (solid lines) and the 1D 1.5-mm reference case (dashed lines). Fig. 11. (a) Structure for YH2 budget terms and (b) temperature gradient for locations A, B, and C in Fig. 9. Propagation is to the left. Dashed line: reaction. Solid line: H2 diffusion. Dotted line: |∇T |.

early in time, heat release rate increases during thermal runaway in ignition kernels, but no surface exists. Note that the two curves do not collapse exactly and deviate from each other in the earlier phase. This is attributed to the fact that a larger fraction of the total heat release rate comes from the ignition kernels during the initial ignition process. To underscore the importance of the local temperature gradient in determining the propagation regime, Fig. 11 shows the reaction–diffusion budget for the evolution of H2 mass fraction. Cut A in Fig. 9 represents a spontaneous ignition front, having a speed of approximately 2.2sL , with a relatively low temperature gradient. These data correspond to a spontaneous ignition front, as in the 4.1 mm 1D case in Fig. 4. The effects of diffusion are found to be much smaller compared to reaction, as was observed in the 1D test case. Cut B, on the other hand, shows a front which has a speed approximately equal to sL . In this case the temperature gradient and diffusion–reaction balance are comparable to the 0.75-mm 1D case (Fig. 4), which was found to be diffusion-controlled. Finally, cut C shows an area with a much larger temperature gradient. At this location, where sd∗ ≈ sL , diffusion nearly counterbalances reaction, and the front behaves much like a normal deflagration. For the mean and RMS temperature fluctuation considered in the 1D test cases, such a high gradient could not be obtained, pointing to a possible influence of turbulence in modifying the local temperature gradient, as will be clarified later in this study. The evolution of the ignition front speed and the magnitude of the temperature gradient on the front is monitored by tracking sd∗ and |∇T |, averaged on the

YH2 = 8.5 × 10−4 isocontour. A comparison of the front propagation for the 2D DNS and the 1.5-mm 1D test case is presented in Fig. 12. It is first noted that the 2D case ignites earlier. This is due to higher initial temperatures that occur owing to the random nature of the initialization. The front speed curve for the 2D and 1D test cases both show the same qualitative Ushape, with the 2D case showing a longer duration of combustion because of a wider distribution in the temperature. The 2D case attains higher mean speeds compared to the 1D test case. This is expected since, throughout the duration of combustion, new kernels form and ignite, and fronts merge and annihilate as the ignition front density reaches a maximum at the time of volumetric peak heat release rate. Both ignition and annihilation events involve very high speeds [23], so the mean will also be higher. Finally, the temperature gradients extracted at the isosurface show that the 2D case shows a higher average gradient relative to the 1D test case, indicating the influence of turbulent straining, as will be addressed later in this paper. Fig. 12 suggests that the variation in the instantaneous displacement speed is correlated with the local temperature gradient. To confirm this observation, the instantaneous solution field for the 2D DNS at 0.9τ0 is analyzed and sd versus |∇T | at all points along the isosurface is plotted in Fig. 13 on a log scale. A power fit to these data reveals an exponent of −1.005, very close to the value of −1. This is analogous to the theory by Zel’dovich [4], where spontaneous front propagation speed is inversely proportional to the temperature gradient. The dependence of sd∗ on |∇T | was also confirmed by the 1D test case study in the spontaneous ignition regime with moderate values of |∇T |. The results suggest that the inverse linear correlation between the two variables is a good indication of the spontaneous propagation regime.

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

139

Fig. 13. Displacement speed (sd ) shown on a log scale against |∇T | (inset: (sd∗ /sL ) plotted on a regular scale), and power-law fit at 0.9τ0 of the 2D DNS data.

On the other hand, the inset in Fig. 13 further reveals that at approximately |∇T | = 1800 K/mm and above [13], sd∗ does not decrease indefinitely but rather levels off to a minimum value, consistent with the 1D test result (recall the 1.5- and 0.75-mm cases). This is attributed to the reactive–diffusive balance that is intrinsic to a laminar deflagration wave. Note that the spontaneous front speed is determined by the ignition delay differential along neighboring mixture parcels, which is directly related to the local temperature gradient. In a deflagration wave, however, an increase in scalar gradient leads to an increase in diffusive transport and subsequent reaction rate due to the reaction–diffusion balance mechanism. Therefore, the front can self-adjust and the speed remains relatively constant over a wide range of temperature gradients, unless perhaps an extreme amount of transport can cause flame extinction. The weaker sensitivity of the front speed to the temperature gradient appears to be another distinct characteristic of the deflagration front. To determine what fraction of the volume is undergoing combustion in the deflagration versus spontaneous ignition regimes, the local normalized front speed is used as a conditioning variable for the front length and the rate of production of burnt gases. Fig. 14 shows the fraction of the ignition front area exhibiting a speed comparable to the deflagration speed, specifically having a front speed <1.1sL . The choice of the value of 1.1 is meant to give a qualitative sense of the transition from deflagration to spontaneous ignition. In fact, transition occurs more gradually. This figure also shows the fraction of production rate of area of “burnt” gases (having YH2 < 8.5 × 10−4 ) due to deflagration for this case. 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.

Fig. 14. Fraction of ignition front length undergoing deflagration (solid lines) and fraction of burnt gas area production due to deflagration (dashed lines).

Fig. 15. Evolution of volumetric mean absolute value of diffusion (solid lines) and reaction (dashed lines) for the 2D DNS data.

The fraction of ignition front length initially begins at zero, when ignition kernels dominate, and increases as combustion in fronts becomes important, reaching a maximum of around 40%. Eventually, the fraction decreases again during the burnout of the end-gas parcels. The production of area of burnt gases is perhaps a more relevant quantity to consider than the fraction of ignition front length, if heat release rate characteristics are of interest. This quantity is less than the length fraction, simply because the speed of the deflagration wave is lower. A comparison can be made of the production of area of burnt gases in this figure with volumetric heat release rate of Fig. 6 in [13], attained from a different measure based on a critical temperature gradient of 1800 K/mm. The comparison yields good qualitative agreement supporting the use of both measures. To further assess the global effect of diffusion without confining attention to a particular isosurface choice, the evolution of the volume average of the absolute value of the YH2 budget terms for the 2D DNS in Fig. 15. The use of the absolute value is necessary since diffusion is conservative in the mean. The

140

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

Fig. 16. Curvature-induced stretch (sd ∇ · N, solid lines) and curvature (∇ · N, dashed lines) versus time for the 2D DNS data.

figure shows that, for the 2D DNS data under study, diffusion is not negligible. At the time of maximum reaction, it is approximately 27% of reaction, roughly consistent with the previous discussion in relation to the front speeds. Qualitatively similar results are obtained considering the balance terms in the energy equation for temperature, showing that diffusion at maximum heat release rate is approximately 17% of reaction. The nature of the ignition front is further examined based on the flame stretch due to the combined effect of curvature and propagation, sd ∇ · N, where ∇ · N is the isosurface curvature, and N is the normal vector to the surface, pointing toward the unburned mixture. (Curvature-induced stretch dominates over tangential strain at these conditions.) The mean curvature of the surface and the associated stretch are plotted against time in Fig. 16. Initially, very large values of stretch are observed as the first kernels ignite: curvature and flame speed are both large. As fronts subsequently propagate outward, the curvature reduces and the mean speed also reduces as higher temperature gradients are encountered. Near the maximum heat release rate time (≈ 0.9τ0 ), the front length reaches a maximum and stretch and curvature both switch sign. Subsequently fronts propagate inward, and finally annihilate leading to large negative stretch. It is of interest to examine what fraction of the stretch is due to spontaneous ignition fronts and deflagrations waves. The fraction of stretch is defined as ratio of the sum of the absolute values of stretch associated with deflagration fronts to that with the entire fronts. In Fig. 17, the fraction of front length marked as deflagration is presented along with the fraction of stretch. Interestingly, while a significant fraction of front length is in the deflagration mode, there is practically no contribution to the total stretch rate—the maximum is less than 5%. To understand this result

Fig. 17. Fraction of ignition front length (solid line) and curvature-induced stretch (sd ∇ · N, dashed line) due to deflagration for 2D DNS at 0.9τ0 .

Fig. 18. Magnitude of the temperature gradient (|∇T |) against curvature (∇ · N) for 2D DNS at 0.9τ0 .

further, the local temperature gradient, on which the speed has been shown to be highly dependent in the ignition front regime, is plotted versus curvature for the time 0.9τ0 in Fig. 18. The data have been averaged over discrete intervals of curvature for clear presentation. The figure shows that low curvatures are typically associated with high temperature gradients, and therefore promote deflagrations, whereas highly curved elements typically have lower temperature gradients, and propagate as ignition fronts. Note also that, because of the low temperature gradients, speeds are very high in these regions, leading to very large stretch rates and the result observed in Fig. 17. This result could have also been qualitatively inferred from Fig. 9, where it was evident that the flame sections marked as deflagrations typically did not occur in highly curved elements. The increased likelihood of low gradients in highly curved elements can be explained by considering the effects of turbulent straining. Compressive strain, which decreases the temperature gradient, leads to curved fronts, whereas extensive strain, which increases temperature gradients, tends to straighten the fronts. This is verified in Fig. 19, showing the positive correlation of tempera-

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

141

Fig. 19. Magnitude of the temperature gradient (|∇T |) against strain (at τ0 ) for 2D DNS at 0.9τ0 . Fig. 21. Mean heat release rate against time for the 2D DNS (solid line) and the multizone model started from the 2D DNS initial condition (MZ 0%, dash–dot line) and at 10% pressure rise (MZ 10%, dashed line).

Fig. 20. Mixing timescale for 2D DNS as defined by Eq. (9).

ture gradient with tangential strain rate on the surface. This also explains why larger temperature gradients are observed in the 2D turbulent case relative to the 1D case. Thus far, discussion has focused on the front propagation characteristics, which mostly involve the inflame diffusion processes. However, it is also important to examine the effect of passive scalar dissipation. Fig. 20 shows the temporal evolution of the mixing time scale (Eq. (10)) for the 2D data. Both quantities are normalized by the homogeneous ignition delay time, τ0 . It is observed that the mixing time scale is initially comparable to the ignition delay, and subsequently reduces as turbulence increases the scalar gradients. Therefore, it is expected that turbulent mixing is important in changing the distribution of initial temperatures in the period prior to ignition.

8. Comparison of 2D DNS with multizone model To assess the effect of passive scalar mixing on the prediction of overall heat release rate, two multizone model simulations were conducted using a large number of zones having the thermochemical state sampled

from points in the 2D DNS data at 0% and 10% pressure rise, denoted as MZ 0% and MZ 10%, respectively. Fig. 21 shows the mean heat release rate history obtained from the DNS and the multizone model calculations. Due to the complete neglect of the initial mixing process and diffusive transport in fronts, the MZ 0% results overestimated the duration of burning by 78% and underestimated the peak heat release rate by 55%. Significantly better predictions are obtained from the MZ 10% case, where the burn duration is overpredicted by 12%, and the peak heat release rate is underpredicted by 16%. This demonstrates the importance of the passive scalar mixing during the early phase of ignition when a considerable level of turbulence is present. The remaining errors in the MZ 10% cases are mainly attributable to transport within fronts. The error values are qualitatively consistent with the fraction of burnt gas production that was identified to occur in deflagrations (cf. Fig. 14). The primary contribution of the present work is to identify the relevant parameters to determine the regime of validity of the models, rather than to make an absolute statement regarding the validity of the model across the range of possible engine operating conditions. The results show that it is possible using the diagnostics developed to explain the model performance. Indeed the results can be used to infer why good performance has been obtained to date with the multizone model. At typical operating conditions, temperature gradients in the bulk gases are low, leading to rapid ignition front propagation relative to slower deflagration waves, and turbulent mixing occurs on longer timescales than the range of ignition delays within the cylinder. However, near walls and crevices, the thermal gradients are likely to be larger than in the bulk gases, and hence, the possibility of increased mixing influences and deflagrative propa-

142

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

gation exists. The diagnostics developed here will be useful in assessing near-wall model performance. From the modeling perspective, the present results are quite encouraging in that, even though diffusive effects are found to be important, the main influence appears to be due to passive scalar mixing—that is, modifying the PDF of initial temperatures. In actual multizone model implementations, these effects are already largely accounted for in the CFD calculation, which includes mixing, that occurs prior to the start of the multizone calculation. Subsequent turbulent mixing may or may not be important, and depends on the Damköhler number. In the case where mixing is important, it is worth noting that modeling of turbulent scalar mixing is relatively well understood and there already exist a number of viable strategies. For example, the unsteady laminar flamelet approach [25] accounts for passive mixing by using a time-varying PDF of mixture fraction fluctuations. The transported PDF approach [24] also adopts the mixing time-scale as a key component of the model. In addition, the passive scalar mixing effect can be easily incorporated in the multizone model. This could be achieved by reweighting the mass in the different zones, similarly to the typical implementation of an unsteady laminar flamelet model. Indeed, a hybrid approach using CFD complete with standard turbulent mixing models to continuously advance the mixing fields coupled with a multizone model to calculate reaction rates [26] has been shown to provide improved predictions of hydrocarbon and CO emissions compared with the original multizone model. In contrast, if deflagrations were found to be significant, modeling would become very challenging since mixing may be controlled through deflagration rather than through the turbulence time scale. Moreover, the deflagration cannot be adequately modeled through the usual flamelet approaches for premixed flames because it occurs in the highly unsteady state corresponding to large temperature gradients. Which type of diffusive effects if any will be important in HCCI applications depends on the turbulence time scales relative to the ignition time under practical operating conditions.

9. Conclusions Direct numerical simulations of lean hydrogen– air ignition at high pressure at constant volume in the presence of temperature inhomogeneities have been conducted to begin to understand the combustion process in HCCI engines. The behavior of the speed and structure of an ignition front was demonstrated in simple one-dimensional constant-volume calculations, illustrating the importance of the initial temperature gradient and

temperature fluctuation length scale on the propagation. The effects of diffusion were shown to be twofold. The first effect is in limiting the minimum attainable front speed, resulting in propagation in a deflagration wave. Based on the theory of Zel’dovich the nondimensional controlling parameter for this transition is identified as the ratio of front speed to the deflagration speed sd∗ /sL . The second diffusive effect is in the dissipation of initial temperature gradients, which on the other hand promotes spontaneous ignition front propagation rather than deflagration. In turbulent flows the latter effect may be curtailed owing to the expected continuous cascade from larger to smaller scales. The importance of this effect can be ascertained with a Damköhler number τmix /τ0 , representing the ratio between the mixing timescale and the ignition delay time. The performance of the multizone model in predictions of the constant volume combustion was evaluated in one-dimensional cases. By testing the model both inside and outside of its expected valid regime we are able to better understand what factors control its validity. The levels of transport by passive scalar mixing and transport in ignition fronts were found to be controlling factors in the performance of the model. Good agreement was obtained in the expected valid regime in which mixing effects are unimportant, and it was noted that the multizone model always predicts a longer burn duration and lower peak heat release rate than the one-dimensional simulations. Two-dimensional simulations were performed with a random temperature distribution and the utility of the front speed as an indicator of the transition between deflagration and spontaneous ignition front propagation was assessed. Analysis of local front structures and global averages showed that this criteria is consistent with the observed importance of diffusion. The case considered here was found to correspond to mixed mode propagation, with a significant fraction of the front length propagating as a deflagration wave. However, since deflagration waves propagate slower, they result in a lower fraction of burnt gas production. It is noted that if the local conditions are below the lean limit, then no importance of deflagrations is expected according to this criteria, since the laminar flame speed is zero. An analysis of the flame stretch showed that a dominant fraction of flame stretch is associated with the spontaneous ignition fronts due to the straining mechanism by turbulence. The conceptual picture that emerges is of rapidly propagating ignition fronts moving along lines of low temperature gradient and leaving behind longer surfaces that show more important diffusive effects, but consume less fuel owing to their lower propagation speed.

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

The comparison between the multizone model and 2D DNS results revealed that the most significant diffusive effects were due to passive scalar dissipation during the early phase of ignition. Transport within ignition fronts was relatively less important but still significant. The explanation of the model predictions demonstrated the utility of the tools developed to understand the combustion process for these conditions. Application of these tools to actual engine conditions will help explain why previous good predictions have been obtained with this model, and help ensure that future applications will be within the expected valid regime. Part II of this study [14] applies the understanding and methods developed here to parametric studies regarding the influence of the initial temperature distribution and turbulence parameters on the combustion mode and the predictions of the multizone model. Future work needs to extend the present study to more realistic chemistry and to engine-relevant turbulence parameters. Fuels exhibiting multistage ignition and negative temperature coefficient regimes may result in different conclusions to those presented here. More information is needed regarding the characteristics of scalar and turbulence fluctuations occurring in real engines operating in HCCI mode in order to design numerical experiments that will be representative. It is anticipated that larger initial integral length scales of the temperature fields more representative of observations in HCCI engines [2] will tend to favor the mode of spontaneous ignition front propagation more than observed here. Considering the wide spectrum of length scales in turbulence, however, there may be some small regions in which diffusive transport is important. Further studies are needed to determine whether they are important for the overall combustion process, including the duration of heat release rate and emissions characteristics.

Acknowledgments The work at UM was supported by the Consortium on HCCI Engine Research directed by the UM and funded by the Department of Energy, and also by Department of Energy, Office of Basic Energy Sciences, SciDAC Computational Chemistry Program. 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

143

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, John Dec, and Magnus Sjöberg.

References [1] K. Epping, S. Aceves, R. Bechtold, J. Dec, SAE Technical Paper 2002-01-1923, 2002. [2] A. Hultqvist, M. Christenson, B. Johansson, M. Richter, J. Nygren, J. Hult, M. Alden, SAE Technical Paper 2002-01-0424, 2002. [3] M. Sjöberg, J.E. Dec, N.P. Cernansky, SAE Technical Paper 2003-01-0113, 2003. [4] Ya.B. Zel’dovich, Combust. Flame 39 (1980) 211–214. [5] D. Bradley, C. Morley, X.J. Gu, D.R. Emerson, SAE Technical Paper 2002-01-2868, 2002. [6] X.J. Gu, D.R. Emerson, D. Bradley, Combust. Flame 133 (2003) 63–74. [7] J.H. Chen, S.D. Mason, J.C. Hewson, in: Proceedings of the Third Joint Sections Meeting of the U.S. Sections of the Combustion Institute, 2003. [8] R. Sankaran, H.G. Im, Combust. Theory Modelling (2005), in press. [9] A. Liñan, F.A. Williams, Combust. Sci. Technol. 105 (1995) 245–263. [10] 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. [11] 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. [12] 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. [13] R. Sankaran, H.G. Im, E.R. Hawkes, J.H. Chen, Proc. Combust. Inst. 30 (2005) 875–882. [14] E.R. Hawkes, R. Sankaran, P.P. Pébay, J.H. Chen, Combust. Flame 145 (2006) 145–159. [15] J.O. Hinze, Turbulence, McGraw–Hill, New York, 1975. [16] C.A. Kennedy, M.H. Carpenter, Appl. Numer. Math. 14 (1994) 397–433. [17] C.A. Kennedy, M.H. Carpenter, R.M. Lewis, Appl. Numer. Math. 35 (2000) 315–357. [18] M.D. Smooke, V. Giovangigli, in: M.D. Smooke (Ed.), Reduced Kinetic Mechanisms and Asymptotic Approximations for Methane–Air Flames, in: Lecture Notes in Physics, vol. 384, Springer-Verlag, New York, 1991. [19] R.J. Kee, F.M. Rupley, E. Meeks, J.A. Miller, Chemkin III: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical and Plasma Kinetics, Tech. Rep. SAND96-8216, Sandia National Laboratories, 1996.

144

J.H. Chen et al. / Combustion and Flame 145 (2006) 128–144

[20] M.A. Mueller, T.J. Kim, R.A. Yetter, F.L. Dryer, Int. J. Chem. Kinet. 31 (1999) 113. [21] T. Echekki, J.H. Chen, Combust. Flame 118 (1999) 308–311. [22] R.J Kee, J.F. Grcar, M.D. Smooke, J.A. Miller, A Fortran Program for Modeling Steady Laminar OneDimensional Premixed Flames, Tech. Rep. SAND858240, Sandia National Laboratories, 1985.

[23] J.H. Chen, T. Echekki, W. Kollmann, Combust. Flame 116 (1999) 15–48. [24] Y.Z. Zhang, E.H. Kung, D.C. Haworth, Proc. Combust. Inst. 30 (2005) 2763–2771. [25] G. Blanquart, H. Pitsch, Proc. Combust. Inst. 30 (2005) 2745–2753. [26] D. Flowers, S. Aceves, J. Martinez-Frias, R. Hessel, R. Dibble, SAE Paper 2003-01-1821, 2003.