Mechanisms of strong pressure wave generation in end-gas autoignition during knocking combustion

Combustion and Flame 162 (2015) 1944–1956

Contents lists available at ScienceDirect

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

Mechanisms of strong pressure wave generation in end-gas autoignition during knocking combustion Hiroshi Terashima a,⇑, Mitsuo Koshi b a b

School of Engineering, The University of Tokyo, 2-11-16 Yayoi, Bunkyo, Tokyo 113-8656, Japan Graduate School of Environment and Information Science, Yokohama National University, 79-7 Tokiwadai, Hodogaya, Yokohama, Kanagawa 240-8501, Japan

a r t i c l e

i n f o

Article history: Received 10 October 2014 Received in revised form 20 December 2014 Accepted 20 December 2014 Available online 23 January 2015 Keywords: Knocking combustion Compressible flow Detailed chemistry Computational fluid dynamics Negative temperature coefficient

a b s t r a c t A knocking combustion modeled using a one-dimensional constant volume reactor is simulated in a manner of direct numerical simulations, in which large detailed chemical kinetic mechanisms for two premixed gases, n-butane (113 species) and n-heptane (373 species), are directly and efficiently introduced. Detailed mechanisms of strong pressure wave generation during end-gas autoignition are clarified. Comparison of n-butane and n-heptane shows that the presence of the negative temperature coefficient (NTC) region significantly influences not only the timing of knocking occurrence but also the amplitude of pressure oscillations. In the case of n-heptane with the condition of an adiabatic wall, there is one large peak produced in the strength of knocking intensity for initial temperature between 450 and 1000 K, whereas there is no peak produced in the case of n-butane. The peak generated around 650 K is attributed to a pressure wave intensified through propagation in the end-gas, which is locally generated near the wall with the influence of the NTC region. It is also found that there is a transition of the autoignition position in the end-gas region from the wall for higher initial temperatures to the region ahead of the flame front for lower initial temperatures, leading to different mechanisms of knocking intensity generation. In the case of the isothermal wall condition, the peak around 650 K is reduced due to the lack of a local temperature increase at the wall, demonstrating the influence of wall temperature conditions on the strength of knocking intensity. Ó 2014 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction Knocking is abnormal combustion in spark-assisted engines and is associated with severe pressure oscillations [1]. The risk of knocking has been a large obstacle to increasing compression ratios and thus operation efficiency. One of the primary causes of the severe pressure oscillations lies in the interaction of end-gas autoignition and the propagation of a flame front from a spark igniter, as sketched in Fig. 1. However, despite the numerous studies devoted to researching knocking combustion, some key physical mechanisms, such as those that determine whether mild or strong knocking combustion occurs, still remain unclear. This lack of understanding is mainly because end-gas autoignition, which initiates pressure oscillations, very rapidly occurs in an extremely localized region (called a hot spot). Such a small-scale transient phenomenon is considerably difficult to measure, even with various experimental techniques [2–5]. Although numerical simulation is a promising tool for the resolution of detailed flow ⇑ Corresponding author. E-mail address: [email protected] (H. Terashima).

structures in knocking combustion, a rigorous yet efficient coupling between fluid dynamics and detailed chemical kinetics, especially for large hydrocarbon fuels, remains a strenuous task even with the latest high-performance computers [6], which unfortunately limits the number of detailed numerical studies involving large hydrocarbon fuels for knocking combustion. As an earlier work on the interaction between a flame front and end-gas autoignition, Pitz and Westbrook [7] examined the transient behavior of a laminar flame propagating into auto-igniting end-gases using the chemical kinetic mechanisms for propane and n-butane. They showed that the high rate of heat release during end-gas autoignition produces strong acoustic waves propagating in the burnt gases. Focusing on an event of end-gas autoignition and the subsequent wave development without flame front propagation from a spark, several detailed numerical studies [8–10] have been carried out to investigate the behaviors of reaction front propagation from hot or cool spots associated with the inhomogeneity of end-gas. Gu et al. [8], following a pioneering work by Zel’dovich [11], identified five modes of reaction front propagation from a hot spot using one-dimensional (1-D) flow equations with detailed chemical kinetic mechanisms of

http://dx.doi.org/10.1016/j.combustflame.2014.12.013 0010-2180/Ó 2014 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

H2–CO–air and H2–air mixtures, showing the dependency of the initial temperature gradient of the hot spot. They also produced a map describing the limits of the developing detonation mode in terms of two dimensionless parameters. Subsequently, Bradley and Kalghatgi [12] and Bradley [13] discussed a theoretical relationship between one of the dimensionless parameters, overpressure, and noise produced by end-gas autoignition using a variety of fuel mixtures. Recently, Dai et al. [9] performed a 1-D simulation including a detailed chemical kinetic mechanism of n-heptane to account for the negative temperature coefficient (NTC) behavior of reaction front propagation. They classified the modes of reaction front propagation by using a new spot condition, i.e., a cool spot with a positive temperature gradient, in which four modes, including different supersonic ignition modes from the previous study [8], were identified with the importance of the NTC regime on knocking occurrence. However, due to the lack of modeling the flame front propagation from a spark in those studies, the mutual interaction between the flame front and end-gas autoignition is not totally considered. Thus, these studies need to assume pressure and temperature conditions in the end-gas region, while deliberately introducing temperature gradients for hot- or cool-spot generation, which are generally difficult to ascertain. Recently, a sequential process of knocking combustion without assuming temperature gradients, which includes flame front propagation, end-gas autoignition, and reaction front propagation, was simulated by Ju et al. [10] using a 1-D simulation with a detailed chemical kinetic mechanism of n-heptane at an equivalence ratio of 0.4, where the flame front dynamically interacted with the end-gas through the propagation. They found six different combustion regimes at near NTC temperatures in the combustion field with significant effects of NTC on flame dynamics, but they unfortunately did not address the transient oscillatory behavior induced by the end-gas autoignition. Martz et al. [14] performed a transient 1-D reaction front simulation during end-gas autoignition. However, the result represents no transient behaviors of pressure waves during end-gas autoignition seemingly due to the lack of consideration of the acoustic wave coupling in their model. Thus, there have been few studies on detailed mechanisms for generation of strong pressure waves caused by end-gas autoignition, which was also pointed out by a recent review paper [15]. The aim of this study, therefore, is to clarify detailed processes and mechanisms guiding strong pressure wave generation induced by end-gas autoignition. A sequential process of knocking combustion (flame front propagation, end-gas autoignition, and reaction front propagation, i.e., pressure wave generation) is fully considered with the compressible Navier–Stokes equations, where detailed chemical kinetic mechanisms for two premixed gases, nbutane and n-heptane, are applied in order to address the effects of NTC regions on the mechanisms of pressure wave generation in end-gas autoignition.

Flame front propagation to wall End-gas (unburnt)

Burnt gas

2. Numerical methods 2.1. Governing equations and models The flow fields are modeled with the compressible Navier– Stokes equations and the thermally perfect gas equation of state, which are written as follows:

@q þ r  ðquÞ ¼ 0; @t @ðquÞ þ r  ðqu  u þ pd  sÞ ¼ 0; @t @E þ r  ½ðE þ pÞu  s  u þ q ¼ 0; @t @ðqY s Þ _ s; þ r  ðqY s u  qDs rY s Þ ¼ x @t N X Ys T; p ¼ qR Ms s¼1

2 3

s ¼ lð2SÞ  lðr  uÞd;

Fig. 1. Schematic of knocking combustion.

ð2Þ ð3Þ ð4Þ ð5Þ

ð6Þ

where l is the viscosity of the mixtures and S is the symmetric strain rate tensor. The heat flux vector q is modeled as N X q ¼ jrT  q hs Ds rY s ;

ð7Þ

s¼1

where j is the thermal conductivity of the mixtures and hs is the enthalpy of each species. The single component viscosities and binary diffusion coefficients are calculated by the standard kinetic theory expression of Hirschfelder [16], and Warnatz’s model [17] is used for the single component thermal conductivities. For mixtures, the thermal conductivity j is modeled with the formula of Mathur et al. [18], and analogous to the thermal conductivity, the mixture-averaged viscosity l is purposely modeled with an empirical approximation [19] for the low computational cost. The mixture diffusion coefficient Ds is given by a mixture-averaged model based on Fick’s Law [20,21]. The Dufour effect for heat flux, the Soret effect and pressure diffusion for diffusion flux are neglected in this study. The Navier–Stokes Eqs. (1)–(4) are solved in the operator-splitting form [22] in order to efficiently handle a wide range of timescales. The fluid and chemical reaction parts are thus solved separately in terms of the time integration. The fluid part in Eqs. (1)–(4) is integrated under the assumption that the chemical reac_ s ¼ 0), while the chemical reaction part is treations are frozen (x ted under the assumption that the volume and internal energy of fluids are constant. The governing equations for the chemical reactions are thus derived as follows:

q

End-gas autoignition

ð1Þ

where q is the density, u is the velocity vector, p is the pressure, d is the unit tensor, s is the viscous stress tensor, E is the total energy (E ¼ qe þ 12 qu  u), q is the heat flux vector, e is the internal energy, _ s is the proY s is the mass fraction, Ds is the diffusion coefficient, x duction rate of each species s. R is the universal gas constant (R = 8.314 J mol1 K1), T is the temperature, and Ms is the molar mass of each species. Here the subscript s ¼ 1  N where N is the total number of species. The viscous stress tensor s is represented by

dY s _ s; ¼x dt N X dT qcv ¼  es x_ s ; dt s¼1

Spark ignition

1945

ð8Þ ð9Þ

where es is the internal energy of each species and cv is the specific heats at constant volume. The spatial derivatives in Eqs. (1)–(4) are

1946

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

neglected in the chemical reactions. The chemical reaction mecha_ s are generated by the Knowlnisms to compute the reaction rate x edge-basing Utilities for Complex Reaction Systems (KUCRS) [23,24]. CHEMKIN-II libraries [25,26] are used to evaluate variables related to thermodynamics, transports, and chemical reactions. Primitive variables such as the temperature, mass fractions, and pressure are exchanged between the fluid Eqs. (1)–(4) and chemical reaction Eqs. (8) and (9) at each time step. A detailed procedure describing the exchange of variables can be seen in [27].

L high-temperature spot: L/40

flame propagation

wall conditions symmetric condition Fig. 2. Schematic of a one-dimensional reactor.

2.2. Numerical method for fluid part For the Navier–Stokes equations without chemical reaction _ s ¼ 0), the numerical flux is evaluated using the Harten– terms (x Lax–van Leer–Contact (HLLC) scheme [28]. Higher-order spatial accuracy is achieved using the Monotone Upstream Centered Schemes for Conservation Law (MUSCL) with primitive variable interpolation and the minmod limiter [29]. The viscous, heat conductivity, and diffusion terms are evaluated by the second-order central differencing. The time integration is done with the thirdorder total variation diminishing (TVD) Runge–Kutta scheme [30]. These numerical techniques have been successfully used for various compressible reactive flow simulations [31,32]. In this study, no sub-grid scale model is used. 2.3. Species bundling for diffusion coefficients In the case of large reaction mechanisms, the cost of calculating transport properties may become a major expense, as compared to that of chemistry integration [33–35]. In this study, in addition to using the simple mixture models of viscosity and thermal conductivity, a species bundling technique [33] is introduced to efficiently compute diffusion coefficients. The relative error i;j of the diffusion coefficients of the ith and jth species is defined by

i;j ¼

     ln Di;k  < ; max   D j;k k ¼ 1; N T min < T < T max

ð10Þ

where Di;j is the binary diffusion coefficient for the ith and jth species. In Eq. (10), temperatures T min ¼ 300 K and T max ¼ 3500 K and a user-specified threshold value  ¼ 0:1 are set throughout this study. With those values, for example, a n-heptane mechanism with 373 species generated by KUCRS [23] is bundled into 21 groups and a n-butane mechanism with 113 species is bundled into 19 groups. 2.4. Numerical method for chemical reaction part In order to efficiently conduct the time integration of the reaction Eqs. (8) and (9) while avoiding possible stiffness, a robust and fast explicit time integration method [36] is introduced. The method is based on a quasi-steady-state assumption and a general formula that preserves conservation laws for any integration operators. The performance and accuracy are comprehensively validated with zero-dimensional ignition problems under a wide range of conditions [36]. A brief procedure of this integration method is given in Appendix A for reference and completeness.

unless otherwise noted. The initial pressure is fixed to p0 ¼ 5 atm in all conditions, while the initial temperature parametrically changes from 450 to 1000 K (maximum of 16 cases) to investigate its effects on knocking occurrence. Two premixed gases of nbutane/air and n-heptane/air mixtures with an equivalence ratio of 1.0 are considered as the working fluids; the n-heptane/air mixture shows a clear NTC region of ignition delay times compared to the n-butane/air mixture. A flat hot kernel of 1400 K and 5 atm with a length of L/40 cm is initially inserted at the left boundary to initiate flame front propagation. The flame front propagation elevates the pressure and temperature in the constant volume reactor, and autoignition may occur in the unburnt end-gas region between the flame front and the wall, followed by the generation and propagation of strong pressure waves. The detailed reaction mechanisms are generated by KUCRS [23]; n-butane consists of 113 species and 426 reactions and nheptane consists of 373 species and 1071 reactions, which are directly used in the present simulations with help of the highly efficient numerical methodology of fast explicit time integration and the species bundling technique. Figure 3 shows ignition delay times against initial temperatures for n-butane/air and n-heptane/ air stoichiometric mixtures under constant-volume adiabatic conditions, where the presence of the NTC region of n-heptane is more clearly highlighted in comparison with the result of n-butane. For validation of the reaction mechanisms by KUCRS, a comparison with available reaction mechanisms and experimental data is presented in Appendix B. KUCRS provides a satisfactory prediction of the ignition delay time for n-heptane, exhibiting a good agreement with experimental data and another detailed reaction mechanism. In the case of n-butane, KUCRS provides relatively small NTC behavior in the low-temperature region compared to other detailed reaction mechanisms, while exhibiting good agreement with experimental data in the high-temperature region. Note that the reaction mechanism for n-butane used in this study predicts the ignition delay times with a relatively small NTC in the lowtemperature region. A uniform grid with a minimum grid spacing of about 22.1 lm (1811 grid points for L = 4.0 cm) is used to resolve the flow fields. A grid convergence study was preliminarily conducted, demonstrating that a grid converged solution in terms of the knocking intensity, which is defined by Eq. (11) below, is obtained with a grid spacing of 22.1 or 15.6 lm. The result is shown in Appendix C.

4. Results and discussions 4.1. Effects of initial temperature on knocking occurrence

3. Computational model and conditions In this study, knocking combustion is modeled using a onedimensional reactor as shown in Fig. 2. Such a model has been used in several previous studies [8–10,14]. The length of the reactor is L ¼ 4:0 cm. The symmetric condition is used for the left boundary and the adiabatic wall condition is assumed at the right boundary

Figure 4 presents the time histories of pressure and temperature at the wall (x ¼ 4:0 cm) for n-butane/air and n-heptane/air mixtures, showing the effects of initial temperature on knocking occurrence. For both cases, the pressure and temperature in the reactor gradually increase due to the adiabatic compression caused by the flame propagation. If the end-gas autoignition occurs before

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

10 1

Ignition delay / s

10 0

5 atm 10 atm 15 atm

10 -1

10 -2

10 -3

10 -4 1.0

1.2

1.4

1.6

Temperature, 1000/T

(a) n-butane/air 10

Ignition delay / s

flame position when knocking occurs is determined as the position when a flame front reaches the maximum distance from the left boundary, which is identified by the trajectory of the maximum heat release rate. In the case of n-butane, the flame position when knocking occurs changes almost linearly with the initial temperature according to the behavior of the ignition delay times of nbutane. For example, with an initial temperature of 1000 K, the autoignition of the end-gas takes place when the flame front reaches approximately x ¼ 1:8 cm. Conversely in the case of n-heptane, the flame position does not change linearly for initial temperatures between 600 and 800 K, which is the range of the NTC region. The trajectory is associated with the curve of ignition delay times as shown in Fig. 3(b). For lower and higher initial temperatures outside the NTC region, linear trajectories are produced in the same way as in the n-butane case. 4.2. Knocking intensity

1

10 0

The results in the previous section showed that the oscillatory profiles are produced by end-gas autoignition and the timings of knocking occurrence are significantly influenced by the NTC region. Here, we focus on the amplitudes of pressure oscillations, i.e., how pressure waves with large amplitudes are produced. To quantify the amplitudes of pressure oscillations, the knocking intensity is defined as

5 atm 10 atm 15 atm

10 -1

10 -2

KI ¼ p 10 -3

10

1947

-4

1.0

1.2

1.4

1.6

Temperature, 1000/T

(b) n-heptane/air Fig. 3. Ignition delay times for stoichiometric premixed gases under constantvolume adiabatic conditions for which the ignition delay time is determined at the moment when the temperature exceeds 1500 K; the mechanisms are generated by KUCRS [23].

the flame front reaches the wall, the pressure and temperature in the end-gas region rapidly increase and pressure waves with non-negligible amplitudes are generated, which eventually lead to transient oscillatory profiles. In the case of n-butane shown in Fig. 4(a), the timing of knocking occurrence is almost proportional to the initial temperature values: the higher the initial temperature, the earlier the timing of the knocking occurrence. The two lower initial temperature cases of 500 and 600 K show knock-free combustion; the pressure and temperature reach constant values without significant oscillatory behaviors, indicating that the flame front reaches the wall before the end-gas autoignition takes place. Conversely, transient pressure and temperature oscillations are produced by the endgas autoignition for initial temperatures higher than 700 K. In contrast, in the case of n-heptane shown in Fig. 4(b), the presence of the NTC region significantly affects the timing of knocking occurrence and the characteristics of profiles; the timing to generate transient pressure and temperature oscillations is almost the same for temperatures between 700 and 800 K, and a two-stage temperature rise due to low temperature oxidation is distinctly observed in the temperature profiles of 600 and 700 K. Due to the shorter ignition delay times of n-heptane as compared to nbutane (see Fig. 3), n-heptane produces oscillatory profiles for most computational conditions. Figure 5 compares flame positions when end-gas autoignition takes place (i.e., when transient oscillations start in Fig. 4). The

pmax ; pe

ð11Þ

where pe is the equilibrium pressure determined from the time history of the maximum pressure and pmax is the first peak of the maximum pressure history as shown in Fig. 6. Note that the second or third peak of the maximum pressure history may be larger than the first peak pmax because of wave interactions at the symmetric or wall boundary, which are not taken into account in this study. Figure 7 shows the results of the knocking intensity estimated with Eq. (11) against initial temperatures for cases of both nbutane and n-heptane. For the case of n-butane, the knocking intensity constantly changes for temperatures between 700 and KI ¼ 1:45, showing a weak depen1000 K with an average value of p dency of the knocking intensity on the initial temperature. In contrast, there is one significant peak in the knocking intensity produced in the n-heptane case. The peak around 650 K with a KI  2:43 corresponds with the maximum knocking intensity of p NTC region, indicating that the presence of the NTC region significantly affects not only the timing of knocking occurrence (Fig. 5) but also the knocking intensity, i.e., the amplitude of the pressure oscillations. 4.2.1. Strong knocking intensity around 650 K To identify sources of large knocking intensity in the n-heptane case, detailed behaviors of pressure and temperature waves during end-gas autoignition are investigated. Figure 8 shows the pressure and temperature distributions in the x–t plane for an initial temperature of 650 K, where the time increases from the bottom to the top. Notice that only the results during a very short period of approximately 40 ls are presented. The results show that an autoignition event starts at a small region near the wall, generating distinct pressure and temperature waves that propagate toward the left boundary with an almost constant speed. The temperature distribution in Fig. 8(b) shows that the flame front gradually moves back due to end-gas autoignition, i.e., the expansion of the end-gas, and interacts with the temperature wave in the middle of the reactor. After the interaction, the pressure and temperature waves continue to propagate throughout the reactor, resulting in the generation of oscillatory pressure and temperature profiles.

1948

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

40

2500

Temperature / K

Pressure /atm

3000

500 K 600 K 700 K 800 K 900 K 1000 K

50

30 20

2000 1500 500 K 600 K 700 K 800 K 900 K 1000 K

1000 10

500 0

5

10

15

20

25

30

0

35x10-3

5

10

15

Time / s

20

25

30

35x10-3

Time / s

50

3000

40

2500

Temperature / K

Pressure /atm

(a) n-butane/air

30 450 K 500 K 600 K 700 K 800 K 900 K 1000 K

20

10 0

5

10

15

20

25

30

2000 1500

450 K 500 K 600 K 700 K 800 K 900 K 1000 K

1000 500

35x10

0

-3

5

10

15

Time / s

20

25

30

35x10-3

Time / s

(b) n-heptane/air Fig. 4. Effects of initial temperatures on transient histories at the wall (x ¼ 4:0 cm) for two premixed gases (left: pressure; right: temperature). (For interpretation of the references to color in this figure legend, refer to the web version of this article.)

2nd peak* 1000

pmax

1st peak

3rd peak*

pe

800 700

Maximum pressure

Temperature / K

900

600

n-C4 H10 500

n-C7 H16

400

*2nd and 3rd peaks from wave interactions at wall and symmetric conditions

0

1

2

3

4

x coordinate / cm Fig. 5. Flame positions when autoignition of the end-gas takes place.

Figure 9 is related to Fig. 8, in that it depicts a sequence of pressure and temperature profiles with numbers labeling each profile that correspond to those in Fig. 8. It can be seen clearly in Fig. 9(b) that, in the end-gas, autoignition takes place in a small region near the wall earlier than in the other regions, by which an inhomogeneity of temperature distributions is promoted and a considerable temperature gradient is formed near the wall (profile 1). The temperature near the wall then rapidly and locally increases, and pressure and temperature waves are generated (profiles 2 and 3). The generated pressure wave propagates to the left

Time Fig. 6. Schematic of the time history of the maximum pressure for the knocking intensity with pmax and pe indicated.

in the unburnt end-gas region, while being amplified due to the heat release of the chemical reaction behind it. The presence of chemical reactions associated with the pressure wave propagation is indicated by a sequential profile of Y OH and Y CH2 O in Fig. 10. The peak pressure value thus increases, reaching a maximum at t ¼ 8:481 ms (profile 5). After that, the peak pressure decreases because autoignition also occurs in the end-gas region ahead of the pressure wave (profile 6) and interacts with the burnt gases (profiles 7 and 8).

1949

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

8.50

8.49

7

3

t x10 / s

8

6 5

8.48

4 3

50

2

40 30

8.47

atm

KI against initial temperatures for n-butane and Fig. 7. Knocking intensity p n-heptane.

Pressure wave

1

20 10

ITDðtÞ ¼

1 T0

Z

t

ðT x¼4:0  T x¼3:8 Þdt;

ð12Þ

0

where T 0 is the initial temperature, T x¼4:0 is the temperature on the wall (x ¼ 4:0 cm), T x¼3:8 is that at x ¼ 3:8 cm, and t is the elapsed time. ITD is considered to be an indicator of the accumulation of temperature difference between the two points: the larger the ITD, the stronger the inhomogeneity of temperature distributions in the end-gas; the gas at the wall is more exposed with higher temperature than that at the other point. Figure 14 shows that the time histories of ITD for the cases of 650 and 850 K, and the result of 500 K is also plotted for reference

0

1

2

3

4

x /cm

(a) Pressure

8.50

Interaction

Temperature wave

3

t x10 / s

8.49

8.48 3200 2800 2400 2000

8.47

K

The wave speed, which is defined by the position of the maximum pressure peak, becomes supersonic at approximately t ¼ 8:478 ms (profiles 4 and 5) until the pressure wave reaches the left boundary as shown in Fig. 11, where the sound speed is measured from a point ahead of the pressure wave. However, the pressure wave is not developed as a detonation wave, since no more combustible gases exist in the reactor. This propagation mode may be termed as a developing detonation, as in a previous study [8]. In contrast, such a clear and strong pressure wave is not observed in the end-gas region for an initial temperature of 850 K as shown in Figs. 12 and 13. The pressure and temperature in the end-gas region almost uniformly increase with no appearance of a small hot spot, which is highlighted by the gradients dt=dx between ‘‘Autoignition’’ and ‘‘Interaction’’ drawn in Figs. 8(b) and 12(b). Thus, a steep pressure wave is not generated in the end-gas or further amplified with chemical reaction, resulting in a relatively weak knocking intensity. On the other hand, the rapid and almost uniform increase of pressure and temperature in the end-gas region produces pressure and temperature waves with high amplitudes (profiles 5 and 6), which propagate in the reactor, resulting in the oscillatory profiles. According to the previous study [8], this propagation mode is categorized as a rapid autoignitive deflagration, which is also seen in other cases with higher initial temperatures or using n-butane, in which the knocking intensity is relatively low as shown in Fig. 7. The above comparison between the results of 650 and 850 K indicates that large knocking intensity can be attributed to the generation of a strong pressure wave in the end-gas, which is caused by the rapid temperature increase of a small hot spot. Thus, a question arises concerning the appearance of such a hot spot near the wall. To address it, an integral of dimensionless temperature difference (ITD) for identifying mechanisms of hot spot generation in the end-gas is defined as:

Autoignition

1600 1200 800

Flame front

0

1

2

x /cm

3

4 End-gas

(b) Temperature Fig. 8. Pressure and temperature contours in the x–t plane for an initial temperature of 650 K in the n-heptane case. A short period between 8:462  103 and 8:503  103 s is plotted.

(the case of 500 K will be discussed later). The results show that ITD increases very slowly but steadily until end-gas autoignition, indicating that the temperature on the wall (x ¼ 4:0) becomes slightly larger than that on a representative position (x ¼ 3:8) in the end-gas. Therefore, since the gas on the wall is slightly more reactive than the other gases, autoignition in the end-gas generally starts on the wall, followed by a delayed ignition at other points. A comparison between the cases of 650 and 850 K reveals that, while the case of 850 K only shows a single rise of ITD with autoignition, a two-stage rise of ITD appears in the case of 650 K because of the presence of low temperature oxidation. The first rise by the low-temperature oxidation accelerates a further temperature

1950

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

60

Pressure / atm

50

40

10 -1

5 6

4 8 7

7

10

3

YOH

1: t = 8.473e-3 s 2: t+Δt 3: t+2Δt 4: t+3Δt 5: t+4Δt 6: t+5Δt 8 7: t+6Δt 8: t+7Δt Δt = 2.0e-6 s

-2

1: t = 8.473e-3 s 2: t+Δt 3: t+2Δt 4: t+3Δt 5: t+4Δt 6: t+5Δt 7: t+6Δt 8: t+7Δt

10 -3

2

30

10 20

-4

1

6

5

4

3

2

1

Δt = 2.0e-6 s

1.0

1.5

2.0

2.5

3.0

3.5

10 -5

4.0

1.0

1.5

2.0

2.5

x /cm

x /cm

(a) Pressure

(a) OH

3.0

3.5

4.0

10 -1 3200

1: t = 8.473e-3 s 2: t+Δt 3: t+2Δt 4: t+3Δt 5: t+4Δt 6: t+5Δt 7: t+6Δt 8: t+7Δt

10 -2

7

2800 10 -3 1: t = 8.473e-3 s 2: t+Δt 3: t+2Δt 4: t+3Δt 5: t+4Δt 6: t+5Δt 7: t+6Δt 8: t+7Δt

2400 2000 1600

4 6

3

YCH2O

Temperature /K

8

2

1.5

2.5

3.0

3.5

4

7

10 -6

2.0

5

Δt = 2.0e-6 s

1

10 -7

4.0

8

1.0

1.5

2.0

2.5

x / cm

x /cm

(b) Temperature

(b) CH2 O

Fig. 9. A temporal sequence of pressure and temperature profiles for an initial temperature of 650 K in the n-heptane case. The numbers labeling each profile correspond to those in Fig. 8.

increase on the wall before a second (primary) ignition, which promotes a larger temperature difference between the wall and other points in the case of 650 K than in the case of 850 K. Therefore, in the case of 650 K, the ignition at the wall takes place sufficiently earlier than that at the other points, which allows not only the generation of a steep pressure wave but also its development through the propagation in the end-gas, where sufficient combustible gases still exist. 4.2.2. Remarks on knocking intensity around 500 K In Fig. 7, it is recognized that there is a small bump in knocking intensity around 500 K, which is likely caused by a transition of the autoignition position in the end-gas region. The pressure and temperature contours in the x–t plane for the case of 500 K are shown in Fig. 15. In contrast to the previous cases, the pressure wave is not generated starting at the wall but in a region ahead of the flame front, which propagates in the end-gas region, resulting in an interaction with the wall. Sequential profiles of the pressure and temperature in Fig. 16 also demonstrate that the autoignition in the end gas takes place in the region ahead of the flame front (profiles 1, 2, and 3), not at the wall. The pressure wave, which propagates toward the wall, is then amplified due to the reflection at the wall (profile 4), resulting in the appearance of a small bump in the knocking intensity around 500 K. This propagation mode is attributed to an end-gas condition that the gas ahead of the flame front is more reactive than that on the wall. As shown in the previous cases, the temperature at the wall tends to be higher than that in the other regions and therefore the ignition tends to start at the wall. However, the pro-

3.0

3.5

4.0

Fig. 10. A sequential profile of the mass fractions of OH and CH2O for an initial temperature of 650 K in the n-heptane case. The numbers labeling each profile corresponds to those in Fig. 8.

2000 Wave speed Sound speed

1600

Speed / m/s

1.0

6

1

10 -5

5

Δt = 2.0e-6 s

1200

10 -4

2 3

1200

800

400 8.48

8.49

8.5x10

-3

Time / s Fig. 11. Comparison of pressure wave speed and the sound speed estimated from a point ahead of the pressure wave.

file of IDT for 500 K in Fig. 14 shows that, while the temperature at the wall is slightly higher than that at other points in the end-gas by approximately t ¼ 5:0 ms, IDT then begins to decrease gradually, indicating that the temperature at another point subsequently becomes higher than that at the wall. This is likely because the flame front propagation is so relatively slow that the heat conduction from the high-temperature flame can affect the temperature in the end-gas, especially the temperature near the flame front. Therefore, a local region ahead of the flame front can be more reactive than other regions including the wall with the influence of the

1951

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

60 5.96

1: t = 5.934e-3 s 2: t+Δt 3: t+2Δt 4: t+3Δt 5: t+4Δt 6: t+5Δt 7: t+6Δt

Pressure / atm

50

t x10 / s

5.95

40

Δt = 2.0e-6 s

30

3

7

20

6

7

Pressure wave

5.94

1.0

4

50 40 atm

30

5

1.5

2.0

2.5

3.0

x /cm

2

(a) Pressure

3.5

4.0

1

3200 6

7

10

5

1

2

3

Temperature /K

2800 0

2 1

3

4

3

20

5.93

6

5

4

x /cm

(a) Pressure

4

2400

1: t = 5.934e-3 s 2: t+Δt 3: t+2Δt 4: t+3Δt 5: t+4Δt 6: t+5Δt 7: t+6Δt

2000

5.96

3

1600

2

Δt = 2.0e-6 s

1

1200 1.0

1.5

2.0

2.5

3.0

3.5

4.0

x / cm

(b) Temperature

5.95

3

t x10 / s

Interaction

Temperature wave

5.94

Fig. 13. A temporal sequence of pressure and temperature profiles for an initial temperature of 850 K in the n-heptane case. The numbers labeling each profile correspond to those in Fig. 12.

1000

3200

650 K

2800

100

2400 K

2000

Autoignition

1600

ITD

5.93

850 K

10

1200 800

500 K

1 0

1

2

3

4

x /cm

0.1

(b) Temperature Fig. 12. Pressure and temperature contours in the x–t plane for an initial temperature of 850 K in the n-heptane case. A short period between 5:927  103 and 5:962  103 s is plotted.

0.01

0

5

10

15

20x10 -3

Time / s Fig. 14. Time histories of integrals of dimensionless temperature differences (ITD).

NTC behavior; the autoignition then takes place at some distance from the flame front. This type of propagation mode is not observed in the case of n-butane, indicating the influence of the NTC region on the knocking intensity. Another case, with an initial temperature of 475 K, shows a similar propagation mode with a relatively small knocking intensity as shown in Fig. 7, whereas the ignition starts at the wall in the case of 550 K and the flame front almost reaches the wall before autoignition in the end-gas in the case of 450 K. Figure 17 shows a sketch indicating the mechanisms of knocking intensity for varying initial temperatures: the transition of the autoignition position in the end-gas region. The strong knocking

intensity occurring at approximately 650 K is caused by a steep pressure wave generated at the wall, which is intensified by the chemical reaction through propagation in the end-gas. For the case of 500 K, although the knocking intensity is not so large, the mechanism of knocking intensity is different from those for higher initial temperatures including 650 K. The knocking intensity at 500 K is attributed to the wall reflection of a pressure wave, which is generated in the region ahead of the flame front, not at the wall. Both pressure waves are produced by a rapid increase in temperature in a very small region, which is promoted by the presence of lowtemperature oxidation.

1952

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

60

50

Pressure / atm

20.55

20.54

40

Δt = 2.0e-6 s

4 5 6

7 8

30

3 2 1

8

Pressure waves

7

3

t x10 / s

1: t = 20.523e-3 s 2: t+Δt 3: t+2Δt 4: t+3Δt 5: t+4Δt 6: t+5Δt 7: t+6Δt 8: t+7Δt

20

6 5

20.53

4

1.0

3

50

1.5

2.0

Reflection on wall

1

3.5

4.0

atm

30

3.0

(a) Pressure

40

20.52

2.5

x /cm

2

20

3200

10

5

20.51 0

1

2

3

Temperature / K

8

4

x /cm

(a) Pressure

7

4

6

2800 1: t = 20.523e-3 s 2: t+Δt 3: t+2Δt 4: t+3Δt 5: t+4Δt 6: t+5Δt 7: t+6Δt 8: t+7Δt

2400

2000

1600

3

Δt = 2.0e-6 s

2 1

20.55 1.0

1.5

2.0

2.5

3.0

3.5

4.0

x / cm

(b) Temperature

Temperature wave

3

t x10 / s

20.54 Fig. 16. A temporal sequence of pressure and temperature profiles for an initial temperature of 500 K in the n-heptane case. The numbers labeling each profile corresponds to those in Fig. 15.

20.53 3200 2800

Autoignition

2400 2000

K

20.52

1600 1200 800

20.51 0

1

Flame front

2

3

4

End-gas

x /cm

(b) Temperature Fig. 15. Pressure and temperature contours in the x–t plane for an initial temperature of 500 K in the n-heptane case. A period between 20:508  103 and 20:558  103 s is plotted.

4.3. Isothermal wall cases While the mechanisms for large knocking intensity have been clarified through the above discussion, the results also showed the importance of conditions at the adiabatic wall. As shown in Fig. 14, in case of higher initial temperatures, the temperature at the wall becomes slightly higher than that in the other region in the end-gas, thus promoting earlier autoignition at the wall. The temperature gap generated in the end-gas produces a strong pressure wave with the peak of knocking intensity at approximately 650 K. Here, in order to address the influence of wall conditions on knocking combustion, additional computations (10 cases for

n-heptane) are carried out using the condition of an isothermal wall where the wall temperature is fixed to the initial temperature. In Fig. 18, the results of the knocking intensity with the isothermal wall condition are plotted with the previous results of the adiabatic wall condition. In the case of isothermal wall, the peak around 650 K seen in the adiabatic wall condition is considerably reduced. Figure 19 shows a temporal sequence of pressure and temperature profiles for 650 K with the isothermal wall condition T w ¼ 650 K. Although the pressure wave is generated by the autoignition in the region near the wall (profiles 1 and 2), the amplitude is not as intensified, with a maximum peak of approximately 40 atm. This is because the autoignition also takes place in the end-gas region ahead of the pressure wave with very slight time delays. Thus, this propagation mode is likely to an intermediate mode between 650 and 850 K with the adiabatic wall condition, demonstrating that a temperature inhomogeneity in the end-gas is weakened with the isothermal wall condition. On the other hand, knocking intensities around an initial temperature of 500 K are almost the same as that with the adiabatic wall condition, showing no indication of the influence of wall conditions. As shown by a temporal sequence of pressure and temperature profiles in Fig. 20, the autoignition first starts in the region ahead of the flame front (profiles 2 and 3) followed by the pressure wave generation (profiles 3 and 4) and reflection at the wall, which is the same as in the results of adiabatic wall condition in Fig. 16. The results of the isothermal wall condition in comparison with those of the adiabatic wall condition demonstrate the influence of wall temperature on the knocking intensity with regards to the

1953

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

reflection on the wall

supersonic wave

pressure wave development

flame front propagation wall

wall

local end-gas autoignition

local end-gas autoignition

500 K

650 K

Fig. 17. Schematic of the mechanisms of knocking intensity for n-heptane with the adiabatic wall condition.

60 2.6

Isothermal wall Adiabatic wall

2.4

50

Pressure / atm

Pmax/P0

2.2 2.0 1.8 1.6

1: t = 8.49e-3 s 2: t + Δt 3: t + 2Δt 4: t + 3Δt 5: t + 4Δt 6: t + 5Δt 7: t + 6Δt 8: t + 7Δt

40

7

Δt = 2.0e-6 s

30

5

6

4

8

3

1.4 20

1.2

2 1

1.0 400

500

600

700

800

900

1.0

1000

1.5

2.0

2.5

3.0

3.5

4.0

x /cm

Temperature / K

(a) Pressure KI against initial temperatures with the isothermal wall Fig. 18. Knocking intensity p condition in the case of n-heptane. The results with the adiabatic wall condition, which are the same as those in Fig. 7, are also plotted for comparison.

3500 8

mechanisms of pressure wave generation. The comparison also suggests that the strong pressure wave around 650 K may be suppressed by controlling the wall temperature, while knocking intensities around 500 K with no dependency on the wall temperature may be reduced by preventing the wave reflection on the wall. On the other hand, it is noted that, if the wall temperature is fixed to a temperature, e.g., 400 K, not to the initial temperature discussed above, an NTC region may exist near the wall with higher initial temperatures, e.g., 1000 K, which potentially leads to the generation of hot spots and thus strong pressure waves with appropriate temperature gradients [8]. 5. Conclusions A detailed numerical simulation was performed for a knocking combustion model using a one-dimensional constant volume reactor and the detailed mechanisms of pressure wave generation in end-gas autoignition have been discussed. Two premixed gases, n-butane and n-heptane, were considered with the detailed chemical kinetic mechanisms. which are directly and efficiently introduced in the compressible flow solver. The results of comparing cases of n-butane and n-heptane under the condition of an adiabatic wall show that the presence of the negative temperature coefficient (NTC) region significantly influences not only the timing of the knocking occurrence but also the strength of the knocking intensity. In the case of n-heptane, one clear peak in the knocking intensity with varying initial temperature are captured. The peak, around an initial temperature of 650 K, is produced by the strong pressure wave generated at the wall, which is caused by the appearance of temperature gap pro-

Temperature / K

3000

7 6

2500 1: t = 8.49e-3 s 2: t + Δt 3: t + 2Δt 4: t + 3Δt 5: t + 4Δt 6: t + 5Δt 7: t + 6Δt 8: t + 7Δt

2000 1500 1000

4

5

3

2 1

Δt = 2.0e-6 s

500 0 1.0

1.5

2.0

2.5

3.0

3.5

4.0

x / cm

(b) Temperature Fig. 19. A temporal sequence of pressure and temperature profiles for an initial temperature of 650 K with the isothermal wall condition.

moted by the low-temperature oxidation (the presence of the NTC region). It is also found that there is a transition of the autoignition position in the end-gas region from the wall for higher initial temperatures to the region ahead of the flame front for lower initial temperatures. Thus, knocking intensities for initial temperatures below approximately 500 K are produced by the reflection of the pressure wave at the wall. The results with isothermal wall condition show that the peak around 650 K is significantly reduced due to the absence of a temperature increase at the wall, demonstrating the influence of wall temperature conditions on the strength of knocking intensity. In summary, this study distinctly demonstrates that a key for producing strong knocking lies in the appearance of temperature

1954

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

where C s is the production rate and ss is the characteristic time of each species. The superscript denotes a predictor step and m ¼ n as an initial step. Dt is the variable time step size, which initially corresponds to Dt used in the fluid part. The mass conservation error is then estimated with a predictor value as

60 1: t = 20.603e-3 s 2: t + ∆t 3: t + 2∆t 4: t + 3∆t 5: t + 4∆t 6: t + 5∆t 7: t + 6∆t 8: t + 7∆t

Pressure / atm

50

40

5

8

Δ t = 2.0e-6 s

30

7

6

4

1





s¼1

3 2

 >   1  107 , the time step size is reduced by 1.0

1.5

2.0

2.5

3.0

3.5

Dt ¼ Dt

4.0

x /cm

(a) Pressure

8

7

6

5

2500

1000

s¼1 Y s

4

1: t = 20.603e-3 s 2: t + ∆t 3: t + 2∆t 4: t + 3∆t 5: t + 4∆t 6: t + 5∆t 7: t + 6∆t 8: t + 7∆t

1500

1.0

1.5

2.0

1

2.5

3.0

3.5

Y s ;

4.0

x / cm

10 1 Ciezki and Adomeit, shock tube Mehl et al. KUCRS

(b) Temperature Fig. 20. A temporal sequence of pressure and temperature profiles for an initial temperature of 500 K with the isothermal wall condition.

inhomogeneities in the end-gas. The presence of the NTC region can promote the development of temperature inhomogeneities through heat release by low-temperature oxidation. Due to the inhomogeneous autoignition in the end-gas, a pressure wave can be generated, not only at the wall but also in the local region ahead of the flame front. The pressure wave generated at the wall can be intensified through the propagation in the end-gas, while one generated at the region ahead of the flame front can be amplified by reflection at the wall.

Ignition delay / s

10 0 10 -1 10 -2 10 -3 10 -4 10 -5 0.6

10 0

     Dt Dt m þ Cm ; Y s ¼ Y m s ss 1  exp  m s exp  m

ss

ðA:1Þ

Ignition delay / s

This work was partially supported by JSPS KAKENHI Grant No. 26390128. We are very grateful to Drs. Yu Daimon and Hiroumi Tani at the Japan Aerospace Exploration Agency for a fruitful discussion during the course of this work.

Although the concept and a detailed procedure of the method are given in [36], we briefly describe the procedure for the completeness of this paper. Consider a time integration from the nth to the ðn þ 1Þth step. Sub-iterations m are conducted between the nth and ðn þ 1Þth steps. A quasi-steady-state assumption (QSSA) [37] is applied to the time integration for the chemical reaction Eqs. (8). The discrete form is written as

1.0

1.2

1.4

1.6

1.8

1.6

1.8

(a) 13.3 atm 10 1

Appendix A. An explicit time integration method [36]

0.8

Temperature, 1000/T

Acknowledgments

ss

ðA:4Þ

which rescales the solutions at the predictor step so that perfect P mass conservation at the ðm þ 1Þth step ( Ns¼1 Y mþ1 ¼ 1) is guarans teed. Note that the procedures above are introduced because the integration using QSSA does not mathematically satisfy the mass P conservation ( Ns¼1 Y s – 1). Any possible instability caused by the mass conservation error can be eliminated. Since the time step size may become very small, it is limited to a user-defined value, e.g.,

3 2

Δ t = 2.0e-6 s

500

ðA:3Þ

1 Y mþ1 ¼ PN s

3000

2000

 

based on the mass conservation error, and applied to the time integration in Eq. (A.1). The solutions are then updated to the ðm þ 1Þth step with the operation

3500

Temperature / K

ðA:2Þ

s

If the error exceeds a user-specified threshold value, i.e.,

20

0

  1  PN Y   s¼1 s  ¼  PN :  Y 

Ciezki and Adomeit, shock tube Mehl et al. KUCRS

10 -1 10 -2 10 -3 10 -4 10 -5 0.6

0.8

1.0

1.2

1.4

Temperature, 1000/T

(b) 41.45 atm Fig. 21. Ignition delay times for stoichiometric n-heptane/air premixed gases under constant-volume adiabatic conditions; the mechanisms are generated by Mehl et al. [38] and KUCRS [23], and experimental data are obtained in shock tubes [39].

1955

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

1  1012 s, to avoid infinite loops. The temperature in Eq. (9) is simultaneously updated as follows:

62.4 μm 31.2 μm 22.1 μm 15.6 μm

2.4

N 1 X e m DY m s : m cv s¼1 s

ðA:5Þ

The above sub-iterative operation continues until the elapsed time reaches a target time, which may correspond to the time step size used in the fluid part. Then, Y mþ1 ! Y nþ1 ; T mþ1 ! T nþ1 and the s s procedure continues from the fluid part.

2.2

Pmax/P0

DT m ¼ 

2.6

2.0 1.8 1.6 1.4 1.2

Appendix B. Validity of the reaction mechanisms

1.0

The ignition delay time, which is a key component for knocking combustion, is used for validation of the reaction mechanisms used in this study. The computation is performed under assuming constant-volume adiabatic conditions and the ignition delay time is determined at the moment when the temperature exceeds 1500 K. Figure 21 shows the ignition delay times of n-heptane calculated by two reaction mechanisms [38,23] in comparison with shock tube data [39] for two initial pressure conditions. The results by the reaction mechanism used in this study, KUCRS, exhibit very good agreement with the experimental data and the computational data by the other reaction mechanism, demonstrating a satisfactory prediction of the ignition delay time with the NTC behavior.

10 1

Ignition delay / s

10 0 10 -1

400

500

600

700

800

900

1000

Temperature / K KI against initial Fig. 23. Grid convergence study in terms of the knocking intensity p temperatures in the case of n-heptane with the adiabatic wall condition.

Figure 22 also shows a comparison of the ignition delay times of n-butane among several reaction mechanisms [40–43,23] and shock tube data [42] for two initial pressure conditions. The overall trend in the ignition delay time by KUCRS is similar to those by the mechanisms of Cord et al. [43] and Healy et al. [42] for two pressure conditions, and the ignition delay times in the high-temperature region are in good agreement with the shock tube data, demonstrating sufficient performance of KUCRS for n-butane. On the other hand, KUCRS provides relatively small NTC behavior compared to the mechanism of Healy et al. [42] in the low-temperature region. In this study, although there is an uncertainty associated with the reaction mechanisms for n-butane in terms of the NTC behavior, the mechanism by KUCRS is used for consistency with n-heptane.

10 -2 Healy et al., shock tube Marinov et al. Wang et al. Healy et al. Cord et al. KUCRS

10 -3 10 -4 10 -5 0.6

0.8

1.0

1.2

1.4

1.6

1.8

Temperature, 1000/T

(a) 10 atm

Appendix C. Grid convergence study A grid convergence study was preliminarily conducted for the case of n-heptane with the adiabatic wall condition. Figure 23 shows the knocking intensities, which is defined by Eq. (11), with four grid spacings, 62.4, 31.2, 22.1 and 15.6 lm. The results demonstrate that a grid converged solution is almost achieved with a grid spacing of 22.1 or 15.6 lm. Coarser grids tend to generate spurious larger knocking intensities at initial temperatures of 475 and 500 K.

10 1

References

Ignition delay / s

10 0 10 -1 10 -2 Healy et al., shock tube Marinov et al. Wang et al. Healy et al. Cord et al. KUCRS

10 -3 10 -4 10 -5 0.6

0.8

1.0

1.2

1.4

1.6

1.8

Temperature, 1000/T

(b) 19 atm Fig. 22. Ignition delay times for stoichiometric n-butane/air premixed gases under constant-volume adiabatic conditions; the mechanisms are generated by Marinov et al. [40], Wang et al. [41], Healy et al. [42], Cord et al. [43], and KUCRS [23]. Experimental data are obtained in shock tubes [42]. Note that the mechanism for 10 atm is also used for the condition of 19 atm in the result of Cord et al. [43], while the mechanism is provided for 1 and 10 atm.

[1] J.B. Heywood, Internal Combustion Engine Fundamentals, McGraw-Hill, New York, 1988. [2] J. Gabano, T. Kageyama, F. Fisson, J. Leyer, Symposium (International) on Combustion, vol. 22, Elsevier, 1989, pp. 447–454. [3] M. Pöschl, T. Sattelmayer, Combust. Flame 153 (4) (2008) 562–573. [4] N. Kawahara, E. Tomita, Int. J. Hydro. Energy 34 (7) (2009) 3156–3163. [5] C. Strozzi, A. Mura, J. Sotton, M. Bellenoue, Combust. Flame 159 (11) (2012) 3323–3341. [6] T. Lu, C.K. Law, Prog. Energy Combust. Sci. 35 (2) (2009) 192–215. [7] W.J. Pitz, C.K. Westbrook, in: Dynamics of Reactive Systems Part II: Modeling and Heterogeneous Combustion, Progress in Astronautics and Aeronautics series, vol. 105, American Institute of Aeronautics and Astronautics, 1986, pp. 115–130. [8] X. Gu, D. Emerson, D. Bradley, Combust. Flame 133 (1) (2003) 63–74. [9] P. Dai, Z. Chen, S. Chen, Y. Ju, Z. Chen, in: Proceedings of the Combustion Institute, in press, http://dx.doi.org/10.1016/j.proci.2014.06.102. [10] Y. Ju, W. Sun, M.P. Burke, X. Gou, Z. Chen, Proc. Combust. Inst. 33 (1) (2011) 1245–1251. [11] Y.B. Zeldovich, Combust. Flame 39 (2) (1980) 211–214. [12] D. Bradley, G. Kalghatgi, Combust. Flame 156 (12) (2009) 2307–2318. [13] D. Bradley, Philos. Trans. Roy. Soc. A: Math., Phys. Eng. Sci. 370 (1960) (2012) 689–714. [14] J.B. Martz, G.A. Lavoie, H.G. Im, R.J. Middleton, A. Babajimopoulos, D.N. Assanis, Combust. Flame 159 (6) (2012) 2077–2086.

1956

H. Terashima, M. Koshi / Combustion and Flame 162 (2015) 1944–1956

[15] J. Pan, G. Shu, H. Wei, Combust. Sci. Technol. 186 (2) (2014) 192–209. [16] J. Hirschfelder, C. Curtiss, R. Bird, Molecular Theory of Gases and Liquids, John Wiley and Sons, Inc., New York, 1954. [17] N. Peters, J. Warnatz, Numerical Methods in Laminar Flame Propagation: A GAMM-Workshop, vol. 6, Informatica International, Inc., 1982. [18] S. Mathur, P. Tondon, S. Saxena, Molecul. Phys. 12 (6) (1967) 569–579. [19] J. Warnatz, U. Maas, R.W. Dibble, Combustion: Physical and Chemical Fundamentals, Modeling and Simulation, Experiments, Pollutant Formation, Springer, 2006. [20] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, vol. 2, Wiley, New York, 1960. [21] K. Sutton, P.A. Gnoffo, in: Proceedings of 7th AIAA/ASME Joint Thermophysics and Heat Transfer Conference, AIAA98-2575, 1998. [22] G. Strang, SIAM J. Numer. Anal. 5 (3) (1968) 506–517. [23] A. Miyoshi, KUCRS Software Library for update information. The program uses THERM program for thermodata generation. [24] E.R. Ritter, J.W. Bozzelli, Int. J. Chem. Kinet. 23 (9) (1991) 767–778. [25] R. Kee, F. Ruply, J. Miller, Chemkin-II: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical Kinetics, Sandia Report SAND89-8009, 1989. [26] R. Kee, G. Dixon-Lewis, J. Warnatz, M. Coltrin, J. Miller, A Fortran Computer Code Package for the Evaluation of Gas-Phase, Multicomponent Transport Properties, SAND86-8246, 1986. [27] R.P. Fedkiw, B. Merriman, S. Osher, J. Comput. Phys. 132 (2) (1997) 175–190.

[28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41]

[42] [43]

E.F. Toro, M. Spruce, W. Speares, Shock Waves 4 (1) (1994) 25–34. B. Van Leer, Flux-Vector Splitting for the Euler Equation, Springer, 1997. S. Gottlieb, C.-W. Shu, Math. Comput. Am. Math. Soc. 67 (221) (1998) 73–85. H. Terashima, M. Koshi, C. Miwada, T. Mogi, R. Dobashi, Int. J. Hydro. Energy 39 (11) (2014) 6013–6023. H. Tani, H. Terashima, M. Koshi, Y. Daimon, in: Proceedings of the Combustion Institute, in press, http://dx.doi.org/doi:10.1016/j.proci.2014.07.053. T. Lu, C.K. Law, Combust. Flame 148 (3) (2007) 117–126. T. Lu, C.K. Law, C.S. Yoo, J.H. Chen, Combust. Flame 156 (8) (2009) 1542–1551. H. Terashima, Y. Morii, M. Koshi, Int. J. Energ. Mater. Chem. Propul. (2015). Y. Morii, H. Terashima, M. Koshi, T. Shimizu, E. Shima, in: Proceedings of AIAA Propulsion and Energy Forum and Exposition 2014, AIAA2014-3920, 2014. D.R. Mott, E.S. Oran, B. van Leer, J. Comput. Phys. 164 (2) (2000) 407–428. M. Mehl, W.J. Pitz, C.K. Westbrook, H.J. Curran, Proc. Combust. Inst. 33 (1) (2011) 193–200. H. Ciezki, G. Adomeit, Combust. Flame 93 (4) (1993) 421–433. N.M. Marinov, W.J. Pitz, C.K. Westbrook, A.M. Vincitore, M.J. Castaldi, S.M. Senkan, C.F. Melius, Combust. Flame 114 (1) (1998) 192–213. H. Wang, X. You, A.V. Joshi, S.G. Davis, A. Laskin, F. Egolfopoulos, C.K. Law, USC mech version II. High-temperature combustion reaction model of H2/CO/C1–C4 compounds . D. Healy, N. Donato, C. Aul, E. Petersen, C. Zinner, G. Bourque, H. Curran, Combust. Flame 157 (8) (2010) 1526–1539. M. Cord, B. Sirjean, R. Fournet, A. Tomlin, M. Ruiz-Lopez, F. Battin-Leclerc, J. Phys. Chem. A 116 (24) (2012) 6142–6158.