- Email: [email protected]
Computational fluid dynamics modeling of hydrogen ignition in a rapid compression machine
Combustion and Flame 155 (2008) 417–428 www.elsevier.com/locate/combustflame
Computational fluid dynamics modeling of hydrogen ignition in a rapid compression machine Gaurav Mittal ∗ , Mandhapati P. Raju, Chih-Jen Sung Department of Mechanical and Aerospace Engineering, Case Western Reserve University, Cleveland, OH 44106, USA Received 3 March 2008; received in revised form 4 June 2008; accepted 9 June 2008
Abstract In modeling a rapid compression machine (RCM) experiment, a zero-dimensional code is commonly used along with an associated heat loss model. However, the applicability of such a zero-dimensional modeling needs to be assessed over a range of accessible experimental conditions. It is expected that when there exists significant influence of the multidimensional effects, including boundary layer, vortex roll-up, and nonuniform heat release, the zero-dimensional modeling may not be adequate. In this work, we simulate ignition of hydrogen in an RCM by employing computational fluid dynamics (CFD) studies with detailed chemistry. Through the comparison of CFD simulations with zero-dimensional results, the validity of a zero-dimensional modeling for simulating RCM experiments is assessed. Results show that the zero-dimensional modeling based on the approach of “adiabatic volume expansion” generally performs very well in adequately predicting the ignition delay of hydrogen, especially when a well-defined homogeneous core is retained within an RCM. As expected, the performance of this zero-dimensional modeling deteriorates with increasing temperature nonuniformity within the reaction chamber. Implications for the species sampling experiments in an RCM are further discussed. Proper interpretation of the measured species concentrations is emphasized and the validity of simulating RCM species sampling results with a zero-dimensional model is assessed. © 2008 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Rapid compression machine; Hydrogen; Ignition delay; CFD modeling
1. Introduction A rapid compression machine (RCM) typically simulates a single compression stroke of an engine and provides a direct measurement of ignition delay time based on the experimental pressure history. Recognizing the importance of the temperature distribution inside an RCM, many investigators (e.g., [1–5]) * Corresponding author.
E-mail address: [email protected] (G. Mittal).
have highlighted the effect of roll-up vortex caused by the piston movement during the compression stroke and the temperature nonhomogeneity induced by such fluid motion. On the basis of the suggestion of Park and Keck [6], Lee and Hochgreb [2] computationally showed that a properly designed piston crevice can suppress the vortex formation, leading to better definition of reacting core conditions by confining the cold gases to the wall regime. Recently, Mittal and Sung [7] conclusively demonstrated through experiments and simulations that the resulting temperature
0010-2180/$ – see front matter © 2008 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.combustflame.2008.06.006
418
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
nonhomogeneity inside an RCM is very severe when a flat piston is used, whereas a much better homogeneous environment is obtained with the use of an optimized creviced piston. Mittal et al. [8] further conducted a computational fluid dynamics (CFD) study using nonreactive mixtures to investigate the effect of physical and operating conditions on the performance of an RCM with creviced pistons. The work of Mittal et al. [8] emphasized the importance of assessing the creviced piston performance over the range of operating conditions to obtain reliable kinetic data using an RCM. In addition, the individual effects of compressed pressure, clearance volume, and stroke length on effective suppression of vortex are highlighted in [8]. Specifically, the influence of vortex was noted to increase at lower compressed pressures and by using longer stroke length and shorter clearance [8]. To validate chemical kinetic mechanisms with experimental data obtained from RCMs, proper modeling of RCM experiments is required. Typically, modeling of RCM experiments has been conducted by using a zero-dimensional code, such as SENKIN [9], while accounting for the effect of heat loss through empirical parameters based on the approach of Newtonian heat loss or the adiabatic volume expansion model [10–12]. Taking our earlier investigations [13– 15] as examples, for each reactive experiment, a corresponding nonreactive experiment is first conducted by compressing an inert mixture with the same specific heat and under the same operating conditions as the reactive mixture. The empirical heat loss parameters are then derived from the measured pressure history of the nonreactive mixture by using a specific heat loss model. Moreover, the heat loss modeling is considered adequate if the simulated and experimental pressure traces for the nonreactive case match. Subsequently, the modeling of the reactive experiment is conducted by using the same heat loss parameters. Although various heat loss modeling approaches (e.g., [10–12]) can be used to simulate and match the nonreactive pressure history, different approaches can lead to different temperature histories during the postcompression period [8,15]. The approach of Newtonian heat loss implicitly assumes uniform temperature field inside the reaction chamber and does not account for the fact that the boundary layer region is at a lower temperature and the core region is at a higher temperature. In contrast, the approach based on volume expansion assumes that there is no mixing between the cold boundary layer and the hot core region, and the only way the effect of near-wall heat loss penetrates the core region is through expansion of core caused by boundary layer cooling. Therefore, even though the geometric volume of reaction chamber re-
mains unchanged after the piston reaches the top dead center (TDC), the core region experiences an expansion that can be modeled as adiabatic volume expansion. For nonreactive experiments in particular, our CFD simulations have demonstrated the fallacies associated with using the Newtonian heat loss approach, which overcompensates for the heat loss effect, and the adequacy of the volume expansion approach in predicting the maximum temperature history when a creviced piston is used [7,8,15]. The validation of the modeling approach based on “adiabatic volume expansion,” however, has not been performed under reactive conditions. Considering that the phenomena inside an RCM are generally multidimensional in nature and that certain amounts of heat can be released prior to ignition for reactive mixtures, the adequacy of using a zero-dimensional model along with the approach of adiabatic volume expansion needs to be further scrutinized. Especially, the zero-dimensional modeling is based on the assumption of uniform heat release throughout the entire reaction chamber, and could overlook the effect of nonuniform heat release due to the presence of boundary layer and the piston crevice zone. Consequently, it is unclear a priori whether applying the heat loss parameters derived from a corresponding nonreactive experiment to model a reactive experiment is still valid when the extent of pre-ignition reactions become significant. This study aims to address the above-mentioned issues by conducting a CFD-based combustion modeling of a H2 /O2 /N2 /Ar mixture in an RCM. H2 was chosen as a representative fuel because its chemical kinetics is fairly well-established and the small size of the associated reaction mechanism facilitates CFD modeling with detailed chemistry. It is also noted that for the conditions above the extended second limit, the chemistry associated with the formation of HO2 and H2 O2 during the induction period can lead to significant amount of heat release prior to ignition. For such conditions, the influence of pre-ignition heat release on the hot ignition event will be systematically examined. Recognizing that the ignition characteristics of hydrocarbon fuels are different from hydrogen ignition, especially those exhibiting two-stage ignition phenomena, similar CFD studies on hydrocarbon ignition are warranted. In the following, we first describe the RCM setup to be simulated, including the configuration and the associated piston velocity profile. Specifications of the present numerical simulations are detailed next. Subsequently, the CFD results for nonreactive and reactive simulations are discussed, with special emphasis on assessing the validity of zero-dimensional simulations based on the adiabatic volume expansion approach of [13–15].
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
419
Fig. 1. Schematic of (a) the entire RCM, (b) creviced piston head configuration, and (c) typical computational grid distribution at the end of compression stroke. All length dimensions are in mm. In (b), A = 0.50, B = 4.0, C = 0.15, D = 20.0, and E = 1.50.
2. Experimental apparatus The specifications of the RCM we simulate here have been detailed in [7,13,14] and a schematic of this RCM is shown in Fig. 1a. Briefly, this RCM system consists of a driver piston, a reactor piston, a hydraulic motion control chamber, and a driving air tank. The machine is pneumatically driven and hydraulically stopped. Moreover, the driver cylinder has a bore of 127 mm and the reactor cylinder bore is 50.8 mm. This RCM also allows variations of stroke length and clearance volume. The reaction chamber is equipped with sensing devices for measuring pressure and temperature, gas inlet/outlet ports for preparing the reactant mixture, and quartz windows for optical access. Additionally, the machine incorporates an optimized creviced piston head design, as shown in Fig. 1b, to promote a homogeneous reaction zone. A homogeneous mixture of reactants is prepared manometrically inside a mixing tank equipped with a magnetic stirrer. Here we simulate RCM experiments for the stoichiometric mixture of H2 /O2 /N2 /Ar = 12.5/6.25/18.125/63.125 by mole. RCM ignition experiments using this mixture composition were performed previously over a range of compressed pressures and temperatures [13]. Because the mixture composition is fixed and the initial mixture temperature is kept at room temperature, the compressed temperature at TDC is varied through variation of compression ratio.
Based on the measured pressure traces, the procedure to deduce the velocity of the piston during the compression process is described as follows. This piston velocity profile is required by simulations. Because the compression process can be approximated as polytropic in nature, the relation of PC /P0 = (CR)n is used to determine the index of polytropic compression, n, which is taken as a constant. Here, PC is the measured pressure at the TDC, P0 is the initial pressure, and CR is the geometric compression ratio. Once the constant index of polytropic compression is calculated, the volume of the combustion chamber during compression, V (t), can be expressed as V (t) = V0 [P0 /P (t)]1/n , where V0 is the initial volume before compression begins and P (t) is the experimental pressure trace. Subsequently, the velocity of the piston during compression, Vel(t), can be obtained by Vel(t) = [dV (t)/dt]/(π D 2 /4), where D is the diameter of the reaction chamber and dV (t)/dt is the time rate of change of the chamber volume. To obtain dV (t)/dt, a piecewise continuous polynomial (including continuous first derivatives) is fitted to the time history of the calculated volume of the reaction chamber. This polynomial fit for V (t) is then used to determine the derivative dV (t)/dt. Based on this approach, the derived piston velocity is shown in Fig. 2. It is seen that after the initial rapid acceleration, the piston accelerates at a slower rate and reaches a peak velocity. When the TDC is approached, there is rapid but uniform deceleration. For simulations, the piston velocity is approximated by the four step piecewise
420
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Fig. 2. Piston velocity profiles deduced from experimental pressure history (dotted line) and the one used in simulations (solid line).
continuous linear profiles, which are also shown in Fig. 2.
3. Numerical specifications 3.1. Computational fluid dynamics simulations Simulations are conducted for an axisymmetric configuration using the FLUENT CFD package [16]. Computations are performed from the beginning of the compression stroke with a compression time of 33 ms, as illustrated in Fig. 2. Initially (at time = 0), the gas mixture at rest is specified with uniform temperature of 298 K and pressure. A fixed temperature of 298 K and no slip conditions are specified at the cylinder wall boundary and the piston surface. The temperature-dependent thermal conductivity and viscosity of the given gas mixture are calculated using the CHEMKIN subroutines [17, 18] and are specified as a temperature-dependent polynomial in the FLUENT calculations. The kinetic mechanism of H2 and the associated thermodynamic parameters are taken in CHEMKIN format from the model of Li et al. [19], consisting of 11 species and 19 reversible reactions. FLUENT simulations use the segregated implicit solver with the PISO (pressure-implicit split-operator) algorithm for pressure–velocity coupling, the PRESTO (pressureimplicit split-operator) scheme for pressure, and the second-order upwind discretization for density, momentum, and species. In the simulation, the piston starts from rest, and its motion during the compression stroke is given in a manner specified in Fig. 2. At the top dead center, the piston comes to rest and remains there for the subsequent time steps. Throughout the calculations, a time step of 36.667 µs is taken. A typical computational grid distribution at the end of compression is shown
in Fig. 1c. In the axial and radial directions, fine grids are taken near the walls/surfaces. Independence with respect to the size of time step and grid are ascertained by checking calculations on a finer grid with smaller time step size. At every time step, convergence is ascertained by monitoring residues. A much finer time step size is required to accurately resolve temperature and pressure rise during the rapid ignition event. However, the use of such a fine time step does not alter the solution before hot ignition or the value of ignition delay. This is further confirmed by conducting calculations near the ignition event by using a time step size smaller by two orders of magnitude. As such, the very fine resolution of solutions during the entire hot ignition process is deemed unnecessary for the present investigations. All calculations are performed for laminar flow conditions. By comparing the CFD calculations with the experimental measurements, it was shown [7] that laminar simulations adequately describe the temperature and velocity fields inside the current RCM, whereas turbulent simulations based on Reynoldsaverage Navier–Stokes equations smear out the features of temperature field and fail to predict multiple rotating vortices that are inferred to be present from the experimental measurements and are predicted by laminar flow calculations. 3.2. Zero-dimensional simulations In addition to the CFD simulations, zero-dimensional simulations are conducted using SENKIN [9]. As discussed in [14], the SENKIN calculations include the compression stroke and take the effect of heat loss during compression and postcompression period into account through an approach based on volume expansion, which uses the empirically determined heat loss parameters obtained from the nonreactive pressure history. The nonreactive pressure history can be used to back calculate the time-dependent “effective volume” of the core region. The heat loss model developed for the present RCM is briefly described in the following. To obtain the effective volume, an empirically determined parameter, Vadd , is added to the actual time-dependent geometric volume of the combustion chamber, Vg (t), to match the simulated pressure at the end of compression stroke, t = tTDC , with the experimental PC value. Here, Vadd simulates the effect of heat loss during compression. After the end of compression, t > tTDC , volume expansion is specified in terms of a polynomial fit, vp (t). The effective volume after compression is then taken as the product of the effective volume at TDC, Veff (tTDC ), and the fitted
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
421
volume expansion term, vp (t); namely, Veff (t) = Vg (t) + Vadd ,
t tTDC ,
Veff (t) = Veff (tTDC )vp (t), t > tTDC . Therefore, Vadd and the polynomial function vp (t) are the key parameters of the present volume expansion model. For simulating the RCM experiment of a given reactive mixture, a corresponding nonreactive experiment is first carried out with a mixture of the same values of heat capacity, initial pressure, and initial temperature as for the reactive mixture. As such, the parameters required for the present heat transfer model, including Vadd and vp (t), can be determined from the pressure history, P (t), of the nonreactive experiment. The empirically determined parameters are then used for simulating the corresponding reactive case. Numerical calculations are performed using the Sandia SENKIN code [9] in conjunction with CHEMKIN [17]. For modeling the RCM experiment, conditions specified in the SENKIN calculations are those of a closed adiabatic system with the timevarying volume prescribed by Veff (t). In this paper, the primary objective is to compare CFD calculations with the present zero-dimensional model to assess the applicability of such a zerodimensional approach. As done in the experiments, the pressure trace obtained from the corresponding nonreactive CFD simulation is used as the reference pressure trace for deriving the heat loss parameters that are used in the subsequent zero-dimensional simulation. The reactive CFD solutions are then compared with those obtained from the reactive zerodimensional simulation using SENKIN.
4. Results 4.1. Adiabatic constant volume calculations To ensure the predictability of the FLUENT package for reactive simulations, calculations are first conducted for a constant volume adiabatic system and compared with those for the solutions obtained from SENKIN. For these FLUENT calculations, a simple two-dimensional grid with adiabatic wall conditions is used and the reactive mixture is at rest initially. Fig. 3 shows such a comparison with initial pressure of 30 bar and initial temperature of 950 K using the stoichiometric H2 /O2 /N2 /Ar mixture specified earlier. Again, the reaction mechanism used is taken from Li et al. [19]. The FLUENT solutions closely follow the SENKIN solutions in terms of time variations of pressure, temperature, and species. Although a small discrepancy is noted when the time of ignition is approached, the predicted ignition delays of
Fig. 3. Comparison of solutions from FLUENT (line) with SENKIN (circles) for the constant volume adiabatic calculation using the mechanism of Li et al. [19]. Initial pressure = 30 bar, initial temperature = 950 K, and mixture molar composition: H2 /O2 /N2 /Ar = 12.5/6.25/18.125/63.125.
FLUENT (21.07 ms) and SENKIN (22.12 ms) are fairly close. Similar simulations are also conducted for other initial conditions of pressure and temperature, and FLUENT typically predicts an ignition delay 4–5% shorter than the SENKIN solution. Considering the differences in the solvers employed in FLUENT and SENKIN and the precision limitation typically in the CFD package, the FLUENT solution is deemed adequate, especially for the prediction of highly temperature-sensitive ignition delay. 4.2. Nonreactive RCM simulations FLUENT results for nonreactive compression in RCM are obtained here by simply turning off the chemistry during calculations, and are first compared with experimental pressure traces for nonreactive runs. Note that experimental nonreactive pressure traces are obtained by using a mixture of H2 /N2 /Ar with the same specific heat, initial pressure, and initial temperature as the H2 /O2 /N2 /Ar mixture. Through such a nonreactive comparison, we can assess whether the experimentally observed heat loss effect can be properly captured by the FLUENT simulation. Fig. 4 shows a comparison for P0 = 705.8 Torr and T0 = 298 K. The compression ratio and stroke length are 10.32 and 238.15 mm, respectively. Fig. 4
422
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
demonstrates that the simulation matches closely with the experiment except for some discrepancy near TDC. In particular, the FLUENT simulation predicts a slightly lower pressure near TDC. In simulation, the clearance volume at TDC is taken as slightly larger than the geometric volume to compensate for various dead volumes in experiment. The lower simulated pressure at TDC could be due to overcompensation for dead volumes. Simulation also shows a slightly lower heat loss near TDC than experiment, as reflected by lower rate of pressure drop in simulation. Nevertheless, the overall agreement is pretty good. Fig. 5 shows simulated temperature fields at varying times after TDC (tTDC = 33 ms) with fuel oxidation chemistry being suppressed during calculations. For the ease of representation, only the main reaction chamber apart from crevice is shown in Fig. 5. Furthermore, for clarity only a range of temperatures, as indicated in the legends, is shown. Also, the zone that appears white in the temperature field has temperatures lower than the lowest temperature in the legends. At this PC value of 29.08 bar, some effect
Fig. 4. Comparison of experimental (dotted line) and FLUENT-simulated (solid line) pressure traces. P0 = 705.8 Torr and T0 = 298 K.
of vortex on temperature distribution is seen in Fig. 5. From the FLUENT results, Fig. 6 further plots and compares the evolution of Tmax , Tmavmain , and Tmav , respectively representing the maximum instantaneous temperature, the mass averaged temperature of the main chamber (excluding crevice), and the mass averaged temperature of the entire chamber (including crevice). It is seen from Fig. 6 that near and after the TDC Tmavmain is lower than Tmax due to the effect of cold boundary layer and vortex. Hence, the difference between Tmavmain and Tmax indicates the extent of boundary layer. Moreover, Tmav is significantly lower than Tmavmain due to the crevice’s being at a much lower temperature. During compression, simulations also show that the mass fraction inside crevice increases and reaches 0.14 at TDC, which continues to increase gradually during the postcompression period due to the flow of mass into the crevice.
Fig. 6. Temperature profiles from the nonreactive FLUENT simulation. Tmax = maximum instantaneous temperature, Tmavmain = mass averaged temperature of the main reaction chamber (without crevice), Tmav = mass averaged temperature of the entire reaction chamber (including crevice). P0 = 705.8 Torr, T0 = 298 K, and PC = 29.08 bar. Mixture molar composition: H2 /O2 /N2 /Ar = 12.5/6.25/18.125/63.125.
Fig. 5. Comparison of simulated temperature fields within the main reaction chamber (excluding crevice) at varying times after the end of compression. P0 = 705.8 Torr, T0 = 298 K, and PC = 29.08 bar. The orientation of representation is identical to that in Fig. 1c. Mixture molar composition: H2 /O2 /N2 /Ar = 12.5/6.25/18.125/63.125. In this nonreactive FLUENT simulation, chemical reactions are not considered.
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Based on the above nonreactive FLUENT results, implications for the reactive cases are as follows. Because of the Arrhenius dependence of reaction rate on temperature, the higher temperature core region will experience higher reactivity than the colder regimes such as the boundary layer and crevice zone, leading to nonuniform heat release during reactions. In addition, the heat release within the main chamber could result in an enhanced flow of mass into the crevice. Such issues are not accounted for while using a zerodimensional simulation. The following reactive RCM simulations will help clarify the effect of nonuniform heat release on ignition delay.
Fig. 7. Time variations of Tmax , Tmavmain , and Tmav from the reactive and nonreactive FLUENT simulations. P0 = 705.8 Torr, T0 = 298 K, PC = 28.6 bar, and TC = 964.25 K.
423
4.3. Reactive RCM simulations It has been shown in our previous nonreactive RCM simulations [8] that temperature nonhomogeneity induced by boundary layer and vortex becomes more pronounced under conditions of low compressed pressure. Therefore, reactive RCM simulations are conducted over a range of operating pressures. In the following, two representative cases, PC ∼ 15 and 29 bar, are presented and discussed. 4.3.1. Case 1: PC ∼ 29 bar This simulation is conducted for initial conditions of P0 = 705.8 Torr and T0 = 298 K. The compression ratio, stroke length, and clearance are respectively 10.2, 237.87 mm, and 23.37 mm, yielding compressed conditions of PC = 28.6 bar and TC = 964.25 K. Fig. 7 shows the time variations of Tmax , Tmavmain , and Tmav for both reactive and nonreactive simulations. Half-way through the induction period, the temperature traces for the reactive case begin to deviate from the corresponding nonreactive traces, highlighting the importance of pre-ignition heat release. As time of ignition approaches, at time ∼59 ms, all temperatures for the reactive simulation in Fig. 7 show similar rapid rising patterns. By defining ignition delay as the time interval from the end of compression stroke to the event of ignition, it is ∼26 ms for the condition of Fig. 7. Fig. 8 shows and compares the temperature field and the distribution of OH mole fraction for various
Fig. 8. Temperature fields and OH mole fraction distributions within the main reaction chamber (excluding crevice) at varying times from the reactive FLUENT simulations. P0 = 705.8 Torr, T0 = 298 K, PC = 28.6 bar, and TC = 964.25 K. The orientation of representation is identical to that in Fig. 1c.
424
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Fig. 9. Mass fraction inside the crevice for the reactive and nonreactive FLUENT simulations. P0 = 705.8 Torr, T0 = 298 K, PC = 28.6 bar, and TC = 964.25 K.
times near the point of ignition. For the present case of hydrogen ignition, local ignition is seen to occur at the location of maximum temperature and OH mole fraction. Fig. 9 further presents a plot of the mass fraction that is contained in the crevice for both reactive and nonreactive cases. For the reactive case, due to heat release in the main reaction chamber, there is mass flow into the crevice, and hence mass fraction inside the crevice begins to increase in comparison to the nonreactive case. At the instant of ignition, there is significant mass flow into the crevice. Since our major concern is the effect of this mass flow on altering ignition delay, it is noted that prior to ignition the effect of this mass flow is much less as the mass fraction inside the crevice increases from 0.175 for the nonreactive case to 0.185 for the reactive case at time = 58.9 ms. Figs. 10a and 10b compare the results of SENKIN and FLUENT for the nonreactive and reactive simulations, respectively. For the FLUENT simulation, the maximum instantaneous temperature (Tmax ) is plotted. To obtain the pressure and temperature traces from SENKIN, the nonreactive pressure trace from the FLUENT calculations is first used to back calculate the “effective volume” as described previously. The effective volume is then used for the SENKIN simulation. It is noted that based on this approach not only the pressure but also the maximum instantaneous temperature is also exactly reproduced by the zerodimensional SENKIN solution, as shown in Fig. 10a. For the reactive case, it is seen from Fig. 10b that both calculations show the same features of pressure and temperature variations, and predict practically identical ignition delays; 25.94 ms from FLUENT and 26.564 ms from SENKIN. Therefore, despite the temperature nonuniformity and nonuniform heat release, the ignition event is adequately predicted by the zerodimensional model for this case of PC ∼ 29 bar and ignition delay of ∼26 ms. Since the effect of vortex
Fig. 10. Comparison of pressure and temperature traces from FLUENT (lines) and SENKIN (circles): (a) nonreactive simulation and (b) reactive simulation. For the FLUENT results, maximum instantaneous values are plotted. P0 = 705.8 Torr, T0 = 298 K, PC = 28.6 bar, and TC = 964.25 K.
becomes more pronounced when compressed pressure is lower [8], a case of PC ∼ 15 bar with a longer ignition delay is discussed next. 4.3.2. Case 2: PC ∼ 15 bar The simulation described in this section is conducted for P0 = 371.28 Torr, T0 = 298 K, compression ratio of 10.44, stroke length of 238.46 mm, and clearance of 22.81 mm, yielding compressed conditions of PC = 15.2 bar and TC = 967.6 K. The time variations of Tmax , Tmavmain , and Tmav for the reactive and nonreactive simulations are shown and compared in Fig. 11. With ignition occurring at time = 100.25 ms, ignition delay is 67.25 ms for this case. From time = 50 ms onward, the temperature traces for the reactive case begin to deviate from the nonreactive case. In contrast to the high pressure case (cf. Fig. 7), the features of temperature increase in Fig. 11 are somewhat different. Specifically, from time = 66 ms onward, while the maximum instantaneous temperature keeps on increasing, Tmavmain and Tmav continue to decrease due to significant effect of the cooler boundary layer and crevice zone. In addition, whereas the maximum instantaneous temperature shows a gradual rise in temperature before the event of ignition, the two mass averaged tempera-
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Fig. 11. Time variations of Tmax , Tmavmain , and Tmav from the reactive and nonreactive FLUENT simulations. P0 = 371.28 Torr, T0 = 298 K, PC = 15.2 bar, and TC = 967.6 K.
tures exhibit abrupt temperature rise only at the point of ignition. It is also noted that prior to ignition, the mass fraction inside the crevice increases from 0.202 for the nonreactive case to 0.207 for the reactive case. Fig. 12 shows the distributions of temperature and mole fractions of OH and H2 O2 at varying times. The temperature field for the nonreactive simulation is also included in Fig. 12a as a reference. For this relatively low PC condition, the effect of vortex is noted to be more significant than for Case 1. Comparing reactive and nonreactive temperature fields at time = 99 ms (66 ms after the end of compression), the nonuniform heat release is seen to increase temperature nonuniformity. As with the temperature field, there is similar nonuniformity in the mole fraction distributions of radicals. Figs. 12c and 12d also demonstrate that the peak mole fractions of OH and H2 O2 at TDC (time = 33 ms) are in the range of 10−12 and 10−8 , respectively, which are orders of magnitude smaller than those during the induction period and prior to ignition. As with Case 1, the nonreactive simulation using SENKIN for Case 2 also shows the same level of agreement as noted in Fig. 10a. Fig. 13a compares the pressure and temperature traces of FLUENT and SENKIN for the reactive case. It is seen that the effect of nonuniform heat release plays a role as FLUENT predicts an ignition delay slightly longer than the SENKIN calculations; 67.25 ms for FLUENT and 63 ms for SENKIN. Furthermore, in contrast to Case 1, different behavior of temperature and pressure rise is noted here. In the SENKIN solution, temperature and pressure increase rather gradually before the event of ignition. Whereas in the FLUENT simulation, although the maximum instantaneous temperature increases in the same way as SENKIN, the pressure does not show similar behavior and increases abruptly at the point of ignition. The behavior of abrupt pressure rise could be caused by nonuniform heat release
425
in that it lowers the rate of pressure rise, while the local instantaneous temperature may continue to increase due to local heat release, and eventually lead to rapid ignition. The above results for Case 2 demonstrate that although the ignition delay may still be predicted correctly within a few percent by the zero-dimensional simulation, the features of pressure rise can be different. Thus, caution needs to be exercised in conducting mechanism validation for pre-ignition heat release. Specifically, depending on the degree of nonuniformity, the measured pressure traces in RCM experiments may show a slower rate of rise in pressure near the point of ignition than that obtained by the zerodimensional modeling. Fig. 13b further compares the mole fraction profiles of some key radicals from the FLUENT and SENKIN simulations. For the FLUENT results, the maximum instantaneous mole fractions are plotted in Fig. 13b. Except for the discrepancy in ignition delay, the mole fraction profiles from FLUENT match with the SENKIN solutions closely. Again, the radical mole fractions at the end of the compression (time = 33 ms) are orders of magnitude smaller than those during the chemical induction event. 4.4. Comparison of simulated ignition delays The validity of the SENKIN simulation in predicting ignition delay is further examined by comparing to the FLUENT simulation for varying compressed conditions, covering the variation of ignition delay from a few ms to nearly 100 ms. The results are summarized and compared in Fig. 14. The solid line in Fig. 14 represents the situation of both the FLUENTand SENKIN-simulated ignition delays being equal. The pressures and temperatures indicated in Fig. 14 are the maximum instantaneous values at TDC for each case. The zero-dimensional model based on the approach of adiabatic volume expansion performs quite well for the present investigations of hydrogen ignition, and there is some error involved only when ignition delay is very long and compressed pressure is also low. 4.5. Implications for RCM species sampling Apart from the measurement of global combustion property as ignition delay, RCMs can also be used for acquiring species evolution profiles during the chemical induction period by incorporating a sampling apparatus. Such sampling measurements provide detailed experimental data for rigorous validation and development of chemical kinetic mechanisms. However, because the sampling procedure results in mixing of the reacting mixture from all different zones
426
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Fig. 12. Temperature fields and mole fraction distributions of OH and H2 O2 within the main reaction chamber (excluding crevice) at varying times. P0 = 371.28 Torr, T0 = 298 K, PC = 15.2 bar, and TC = 967.6 K. The orientation of representation is identical to that in Fig. 1c.
in an RCM, including the high-temperature core region, the boundary layer, and the crevice, such mixing can pose a difficulty in interpreting the measured species concentration as an average quantity. As such, it is unclear whether it is appropriate to directly com-
pare the measured species profiles with the simulated results using a zero-dimensional model. The present FLUENT simulations can then be used to assess the influence of such mixing on the RCM sampling experiments.
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Fig. 13. Comparison of the time-varying results from FLUENT (lines) and SENKIN (circles): (a) pressure and temperature; (b) mole fractions of key radicals. For the FLUENT results, maximum instantaneous values are plotted. P0 = 371.28 Torr, T0 = 298 K, PC = 15.2 bar, and TC = 967.6 K.
Fig. 14. Comparison of ignition delays obtained from FLUENT and SENKIN under various compressed conditions of PC and TC .
Fig. 15 compares the profiles of the mass-weighted average mass fractions of HO2 , H2 O2 , and H2 O obtained from FLUENT with the SENKIN results for two cases of PC = 15.2 and 32.52 bar. For both compressed pressures, ignition delays are in the range of 63 ms. For PC = 15.2 bar, although Fig. 13 shows that the maximum instantaneous mole fractions match closely between FLUENT and SENKIN, it is seen from Fig. 15a that the average mass fractions from FLUENT are generally lower than those from the
427
Fig. 15. Comparison of mass-weighted average mass fractions of HO2 , H2 O2 , and H2 O from FLUENT (triangles) and SENKIN (circles). (a) P0 = 371.28 Torr, T0 = 298 K, PC = 15.2 bar, and TC = 967.6 K; (b) P0 = 866.3 Torr, T0 = 298 K, PC = 32.52 bar, and TC = 939 K.
SENKIN solutions. Fig. 15a also demonstrates that the ratio of mass fractions of H2 O, HO2 or H2 O2 from SENKIN with the average mass fractions from FLUENT is ∼1.2 near TDC (time = 33 ms). This difference increases with postcompression time due to the increased influence of boundary layer and crevice, with the ratio becoming ∼3 near the ignition event. On the other hand, for PC = 32.52 bar shown in Fig. 15b, this ratio is ∼1.2 at TDC and increases slightly to become ∼1.5 prior to ignition. The extent of nonuniform species concentration is less at higher pressure due to the reduced effects of heat loss and vortex. These observations have significant implications for the species sampling measurements in an RCM. It is clear that RCM sampling experiments are expected to provide data of better quality for conditions of elevated pressures. Furthermore, when a zero-dimensional modeling is used to simulate RCM sampling experiments for the purpose of mechanism validation, proper interpretation of the measured species profiles is very essential. 5. Conclusions Because CFD simulations of RCM experiments with detailed chemistry can be computationally pro-
428
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
hibitive, for the purpose of validating detailed kinetic mechanisms of various fuels, the use of a zero-dimensional modeling approach is practically inevitable. In this work, FLUENT simulations including detailed chemistry are conducted for modeling hydrogen ignition in an RCM with a creviced piston in order to assess the validity of zero-dimensional modeling in ignition delay prediction. The key results of this computational study on hydrogen ignition are as follows: 1. Based on the comparison of FLUENT simulations and SENKIN calculations, it is demonstrated that the zero-dimensional model based on the approach of adiabatic volume expansion generally performs well in adequately capturing the time variations of pressure and maximum instantaneous temperature for the nonreactive cases and correctly predicting the hydrogen ignition delays for the reactive cases. 2. When the effect of temperature nonuniformity becomes significant, particularly for the conditions of low compressed pressures and long ignition delay times, the zero-dimensional model may predict hydrogen ignition delays somewhat shorter than the CFD modeling. However, under such conditions of longer ignition delays, the resulting uncertainty is of few percent. 3. For the hydrogen cases investigated, zerodimensional modeling along with the adiabatic volume expansion approach would be still considered acceptable in terms of mechanism validation for such global combustion property as ignition delay. 4. For species profile validation, on the other hand, there are issues involved with the use of zerodimensional model in that sampling measurements in an RCM are average species concentration levels at a given instant, which can be significantly different from the maximum instantaneous values typically calculated by the zero-dimensional simulations. Therefore, caution should be exercised in interpreting the measured RCM species data and comparing them with the results obtained using a zero-dimensional model. 5. Furthermore, the present FLUENT results indicate that the sampling experiments are expected to provide better data under elevated pressure conditions. Finally, although the conclusions outlined above are expected to be valid for the ignition cases of other hydrocarbons exhibiting single-stage ignition characteristics, they may not be applicable to the conditions of two-stage ignition phenomena. Further CFD studies for hydrocarbon two-stage ignition are required to ascertain the validity of the present zero-dimensional modeling approach for such cases.
Acknowledgments This work was supported by the National Aeronautics and Space Administration under Grant NNX07AB36Z, with the technical monitoring of Dr. Krishna P. Kundu.
References [1] J.F. Griffiths, Q. Jiao, A. Schreiber, J. Meyer, K.F. Knoche, W. Kardylewski, Combust. Flame 93 (1993) 303–315. [2] D. Lee, S. Hochgreb, Combust. Flame 114 (1998) 531– 545. [3] J. Clarkson, J.F. Griffiths, J.P. Macnamara, B.J. Whitaker, Combust. Flame 125 (2001) 1162–1175. [4] J.F. Griffiths, J.P. MacNamara, C. Mohamed, B.J. Whitaker, J. Pan, C.G.W. Sheppard, Faraday Discuss. 119 (2001) 287–303. [5] J. Würmel, J.M. Simmie, Combust. Flame 41 (2005) 417–430. [6] P. Park, J.C. Keck, SAE paper 9000027, 1990. [7] G. Mittal, C.J. Sung, Combust. Flame 145 (2006) 160– 180. [8] G. Mittal, M.P. Raju, C.J. Sung, K. Kundu, submitted for publication. [9] A.E. Lutz, R.J. Kee, J.A. Miller, SENKIN: A Fortran Program for Predicting Homogeneous Gas Phase Chemical Kinetics with Sensitivity Analysis, Report No. SAND 87-8248, Sandia National Laboratories, 1988. [10] C.K. Westbrook, W.J. Pitz, J.E. Boercker, H.J. Curran, J.F. Griffiths, C. Mohamed, M. Ribaucour, Proc. Combust. Inst. 29 (2002) 1311–1318. [11] S. Tanaka, F. Ayala, J.C. Keck, Combust. Flame 133 (2003) 467–481. [12] M. Ribaucour, R. Minetti, L.R. Sochet, H.J. Curran, W.J. Pitz, C.K. Westbrook, Proc. Combust. Inst. 28 (2000) 1671–1678. [13] G. Mittal, C.J. Sung, R.A. Yetter, Int. J. Chem. Kinet. 38 (2006) 516–529. [14] G. Mittal, C.J. Sung, Combust. Sci. Technol. 179 (3) (2007) 497–530. [15] G. Mittal, C.J. Sung, Combust. Flame 150 (2007) 355– 368. [16] FLUENT 6.1 Documentation, Fluent Inc. [17] R.J. Kee, F.M. Rupley, J.A. Miller, Chemkin-II: A Fortran Chemical Kinetics Package for the Analysis of Gas Phase Chemical Kinetics, Report No. SAND89-8009B, Sandia National Laboratories, 1989. [18] R.J. Kee, G.D. Lewis, J. Warnatz, M.E. Coltrin, J.A. Miller, A Fortran Computer Code Package for the Evaluation of Gas-Phase, Multicomponent Transport Properties, Report No. SAND 86-8246, Sandia National Laboratories, 1986. [19] J. Li, Z. Zhao, A. Kazakov, F.L. Dryer, Int. J. Chem. Kinet. 36 (2004) 566–575.
Computational fluid dynamics modeling of hydrogen ignition in a rapid compression machine Gaurav Mittal ∗ , Mandhapati P. Raju, Chih-Jen Sung Department of Mechanical and Aerospace Engineering, Case Western Reserve University, Cleveland, OH 44106, USA Received 3 March 2008; received in revised form 4 June 2008; accepted 9 June 2008
Abstract In modeling a rapid compression machine (RCM) experiment, a zero-dimensional code is commonly used along with an associated heat loss model. However, the applicability of such a zero-dimensional modeling needs to be assessed over a range of accessible experimental conditions. It is expected that when there exists significant influence of the multidimensional effects, including boundary layer, vortex roll-up, and nonuniform heat release, the zero-dimensional modeling may not be adequate. In this work, we simulate ignition of hydrogen in an RCM by employing computational fluid dynamics (CFD) studies with detailed chemistry. Through the comparison of CFD simulations with zero-dimensional results, the validity of a zero-dimensional modeling for simulating RCM experiments is assessed. Results show that the zero-dimensional modeling based on the approach of “adiabatic volume expansion” generally performs very well in adequately predicting the ignition delay of hydrogen, especially when a well-defined homogeneous core is retained within an RCM. As expected, the performance of this zero-dimensional modeling deteriorates with increasing temperature nonuniformity within the reaction chamber. Implications for the species sampling experiments in an RCM are further discussed. Proper interpretation of the measured species concentrations is emphasized and the validity of simulating RCM species sampling results with a zero-dimensional model is assessed. © 2008 The Combustion Institute. Published by Elsevier Inc. All rights reserved. Keywords: Rapid compression machine; Hydrogen; Ignition delay; CFD modeling
1. Introduction A rapid compression machine (RCM) typically simulates a single compression stroke of an engine and provides a direct measurement of ignition delay time based on the experimental pressure history. Recognizing the importance of the temperature distribution inside an RCM, many investigators (e.g., [1–5]) * Corresponding author.
E-mail address: [email protected] (G. Mittal).
have highlighted the effect of roll-up vortex caused by the piston movement during the compression stroke and the temperature nonhomogeneity induced by such fluid motion. On the basis of the suggestion of Park and Keck [6], Lee and Hochgreb [2] computationally showed that a properly designed piston crevice can suppress the vortex formation, leading to better definition of reacting core conditions by confining the cold gases to the wall regime. Recently, Mittal and Sung [7] conclusively demonstrated through experiments and simulations that the resulting temperature
0010-2180/$ – see front matter © 2008 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.combustflame.2008.06.006
418
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
nonhomogeneity inside an RCM is very severe when a flat piston is used, whereas a much better homogeneous environment is obtained with the use of an optimized creviced piston. Mittal et al. [8] further conducted a computational fluid dynamics (CFD) study using nonreactive mixtures to investigate the effect of physical and operating conditions on the performance of an RCM with creviced pistons. The work of Mittal et al. [8] emphasized the importance of assessing the creviced piston performance over the range of operating conditions to obtain reliable kinetic data using an RCM. In addition, the individual effects of compressed pressure, clearance volume, and stroke length on effective suppression of vortex are highlighted in [8]. Specifically, the influence of vortex was noted to increase at lower compressed pressures and by using longer stroke length and shorter clearance [8]. To validate chemical kinetic mechanisms with experimental data obtained from RCMs, proper modeling of RCM experiments is required. Typically, modeling of RCM experiments has been conducted by using a zero-dimensional code, such as SENKIN [9], while accounting for the effect of heat loss through empirical parameters based on the approach of Newtonian heat loss or the adiabatic volume expansion model [10–12]. Taking our earlier investigations [13– 15] as examples, for each reactive experiment, a corresponding nonreactive experiment is first conducted by compressing an inert mixture with the same specific heat and under the same operating conditions as the reactive mixture. The empirical heat loss parameters are then derived from the measured pressure history of the nonreactive mixture by using a specific heat loss model. Moreover, the heat loss modeling is considered adequate if the simulated and experimental pressure traces for the nonreactive case match. Subsequently, the modeling of the reactive experiment is conducted by using the same heat loss parameters. Although various heat loss modeling approaches (e.g., [10–12]) can be used to simulate and match the nonreactive pressure history, different approaches can lead to different temperature histories during the postcompression period [8,15]. The approach of Newtonian heat loss implicitly assumes uniform temperature field inside the reaction chamber and does not account for the fact that the boundary layer region is at a lower temperature and the core region is at a higher temperature. In contrast, the approach based on volume expansion assumes that there is no mixing between the cold boundary layer and the hot core region, and the only way the effect of near-wall heat loss penetrates the core region is through expansion of core caused by boundary layer cooling. Therefore, even though the geometric volume of reaction chamber re-
mains unchanged after the piston reaches the top dead center (TDC), the core region experiences an expansion that can be modeled as adiabatic volume expansion. For nonreactive experiments in particular, our CFD simulations have demonstrated the fallacies associated with using the Newtonian heat loss approach, which overcompensates for the heat loss effect, and the adequacy of the volume expansion approach in predicting the maximum temperature history when a creviced piston is used [7,8,15]. The validation of the modeling approach based on “adiabatic volume expansion,” however, has not been performed under reactive conditions. Considering that the phenomena inside an RCM are generally multidimensional in nature and that certain amounts of heat can be released prior to ignition for reactive mixtures, the adequacy of using a zero-dimensional model along with the approach of adiabatic volume expansion needs to be further scrutinized. Especially, the zero-dimensional modeling is based on the assumption of uniform heat release throughout the entire reaction chamber, and could overlook the effect of nonuniform heat release due to the presence of boundary layer and the piston crevice zone. Consequently, it is unclear a priori whether applying the heat loss parameters derived from a corresponding nonreactive experiment to model a reactive experiment is still valid when the extent of pre-ignition reactions become significant. This study aims to address the above-mentioned issues by conducting a CFD-based combustion modeling of a H2 /O2 /N2 /Ar mixture in an RCM. H2 was chosen as a representative fuel because its chemical kinetics is fairly well-established and the small size of the associated reaction mechanism facilitates CFD modeling with detailed chemistry. It is also noted that for the conditions above the extended second limit, the chemistry associated with the formation of HO2 and H2 O2 during the induction period can lead to significant amount of heat release prior to ignition. For such conditions, the influence of pre-ignition heat release on the hot ignition event will be systematically examined. Recognizing that the ignition characteristics of hydrocarbon fuels are different from hydrogen ignition, especially those exhibiting two-stage ignition phenomena, similar CFD studies on hydrocarbon ignition are warranted. In the following, we first describe the RCM setup to be simulated, including the configuration and the associated piston velocity profile. Specifications of the present numerical simulations are detailed next. Subsequently, the CFD results for nonreactive and reactive simulations are discussed, with special emphasis on assessing the validity of zero-dimensional simulations based on the adiabatic volume expansion approach of [13–15].
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
419
Fig. 1. Schematic of (a) the entire RCM, (b) creviced piston head configuration, and (c) typical computational grid distribution at the end of compression stroke. All length dimensions are in mm. In (b), A = 0.50, B = 4.0, C = 0.15, D = 20.0, and E = 1.50.
2. Experimental apparatus The specifications of the RCM we simulate here have been detailed in [7,13,14] and a schematic of this RCM is shown in Fig. 1a. Briefly, this RCM system consists of a driver piston, a reactor piston, a hydraulic motion control chamber, and a driving air tank. The machine is pneumatically driven and hydraulically stopped. Moreover, the driver cylinder has a bore of 127 mm and the reactor cylinder bore is 50.8 mm. This RCM also allows variations of stroke length and clearance volume. The reaction chamber is equipped with sensing devices for measuring pressure and temperature, gas inlet/outlet ports for preparing the reactant mixture, and quartz windows for optical access. Additionally, the machine incorporates an optimized creviced piston head design, as shown in Fig. 1b, to promote a homogeneous reaction zone. A homogeneous mixture of reactants is prepared manometrically inside a mixing tank equipped with a magnetic stirrer. Here we simulate RCM experiments for the stoichiometric mixture of H2 /O2 /N2 /Ar = 12.5/6.25/18.125/63.125 by mole. RCM ignition experiments using this mixture composition were performed previously over a range of compressed pressures and temperatures [13]. Because the mixture composition is fixed and the initial mixture temperature is kept at room temperature, the compressed temperature at TDC is varied through variation of compression ratio.
Based on the measured pressure traces, the procedure to deduce the velocity of the piston during the compression process is described as follows. This piston velocity profile is required by simulations. Because the compression process can be approximated as polytropic in nature, the relation of PC /P0 = (CR)n is used to determine the index of polytropic compression, n, which is taken as a constant. Here, PC is the measured pressure at the TDC, P0 is the initial pressure, and CR is the geometric compression ratio. Once the constant index of polytropic compression is calculated, the volume of the combustion chamber during compression, V (t), can be expressed as V (t) = V0 [P0 /P (t)]1/n , where V0 is the initial volume before compression begins and P (t) is the experimental pressure trace. Subsequently, the velocity of the piston during compression, Vel(t), can be obtained by Vel(t) = [dV (t)/dt]/(π D 2 /4), where D is the diameter of the reaction chamber and dV (t)/dt is the time rate of change of the chamber volume. To obtain dV (t)/dt, a piecewise continuous polynomial (including continuous first derivatives) is fitted to the time history of the calculated volume of the reaction chamber. This polynomial fit for V (t) is then used to determine the derivative dV (t)/dt. Based on this approach, the derived piston velocity is shown in Fig. 2. It is seen that after the initial rapid acceleration, the piston accelerates at a slower rate and reaches a peak velocity. When the TDC is approached, there is rapid but uniform deceleration. For simulations, the piston velocity is approximated by the four step piecewise
420
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Fig. 2. Piston velocity profiles deduced from experimental pressure history (dotted line) and the one used in simulations (solid line).
continuous linear profiles, which are also shown in Fig. 2.
3. Numerical specifications 3.1. Computational fluid dynamics simulations Simulations are conducted for an axisymmetric configuration using the FLUENT CFD package [16]. Computations are performed from the beginning of the compression stroke with a compression time of 33 ms, as illustrated in Fig. 2. Initially (at time = 0), the gas mixture at rest is specified with uniform temperature of 298 K and pressure. A fixed temperature of 298 K and no slip conditions are specified at the cylinder wall boundary and the piston surface. The temperature-dependent thermal conductivity and viscosity of the given gas mixture are calculated using the CHEMKIN subroutines [17, 18] and are specified as a temperature-dependent polynomial in the FLUENT calculations. The kinetic mechanism of H2 and the associated thermodynamic parameters are taken in CHEMKIN format from the model of Li et al. [19], consisting of 11 species and 19 reversible reactions. FLUENT simulations use the segregated implicit solver with the PISO (pressure-implicit split-operator) algorithm for pressure–velocity coupling, the PRESTO (pressureimplicit split-operator) scheme for pressure, and the second-order upwind discretization for density, momentum, and species. In the simulation, the piston starts from rest, and its motion during the compression stroke is given in a manner specified in Fig. 2. At the top dead center, the piston comes to rest and remains there for the subsequent time steps. Throughout the calculations, a time step of 36.667 µs is taken. A typical computational grid distribution at the end of compression is shown
in Fig. 1c. In the axial and radial directions, fine grids are taken near the walls/surfaces. Independence with respect to the size of time step and grid are ascertained by checking calculations on a finer grid with smaller time step size. At every time step, convergence is ascertained by monitoring residues. A much finer time step size is required to accurately resolve temperature and pressure rise during the rapid ignition event. However, the use of such a fine time step does not alter the solution before hot ignition or the value of ignition delay. This is further confirmed by conducting calculations near the ignition event by using a time step size smaller by two orders of magnitude. As such, the very fine resolution of solutions during the entire hot ignition process is deemed unnecessary for the present investigations. All calculations are performed for laminar flow conditions. By comparing the CFD calculations with the experimental measurements, it was shown [7] that laminar simulations adequately describe the temperature and velocity fields inside the current RCM, whereas turbulent simulations based on Reynoldsaverage Navier–Stokes equations smear out the features of temperature field and fail to predict multiple rotating vortices that are inferred to be present from the experimental measurements and are predicted by laminar flow calculations. 3.2. Zero-dimensional simulations In addition to the CFD simulations, zero-dimensional simulations are conducted using SENKIN [9]. As discussed in [14], the SENKIN calculations include the compression stroke and take the effect of heat loss during compression and postcompression period into account through an approach based on volume expansion, which uses the empirically determined heat loss parameters obtained from the nonreactive pressure history. The nonreactive pressure history can be used to back calculate the time-dependent “effective volume” of the core region. The heat loss model developed for the present RCM is briefly described in the following. To obtain the effective volume, an empirically determined parameter, Vadd , is added to the actual time-dependent geometric volume of the combustion chamber, Vg (t), to match the simulated pressure at the end of compression stroke, t = tTDC , with the experimental PC value. Here, Vadd simulates the effect of heat loss during compression. After the end of compression, t > tTDC , volume expansion is specified in terms of a polynomial fit, vp (t). The effective volume after compression is then taken as the product of the effective volume at TDC, Veff (tTDC ), and the fitted
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
421
volume expansion term, vp (t); namely, Veff (t) = Vg (t) + Vadd ,
t tTDC ,
Veff (t) = Veff (tTDC )vp (t), t > tTDC . Therefore, Vadd and the polynomial function vp (t) are the key parameters of the present volume expansion model. For simulating the RCM experiment of a given reactive mixture, a corresponding nonreactive experiment is first carried out with a mixture of the same values of heat capacity, initial pressure, and initial temperature as for the reactive mixture. As such, the parameters required for the present heat transfer model, including Vadd and vp (t), can be determined from the pressure history, P (t), of the nonreactive experiment. The empirically determined parameters are then used for simulating the corresponding reactive case. Numerical calculations are performed using the Sandia SENKIN code [9] in conjunction with CHEMKIN [17]. For modeling the RCM experiment, conditions specified in the SENKIN calculations are those of a closed adiabatic system with the timevarying volume prescribed by Veff (t). In this paper, the primary objective is to compare CFD calculations with the present zero-dimensional model to assess the applicability of such a zerodimensional approach. As done in the experiments, the pressure trace obtained from the corresponding nonreactive CFD simulation is used as the reference pressure trace for deriving the heat loss parameters that are used in the subsequent zero-dimensional simulation. The reactive CFD solutions are then compared with those obtained from the reactive zerodimensional simulation using SENKIN.
4. Results 4.1. Adiabatic constant volume calculations To ensure the predictability of the FLUENT package for reactive simulations, calculations are first conducted for a constant volume adiabatic system and compared with those for the solutions obtained from SENKIN. For these FLUENT calculations, a simple two-dimensional grid with adiabatic wall conditions is used and the reactive mixture is at rest initially. Fig. 3 shows such a comparison with initial pressure of 30 bar and initial temperature of 950 K using the stoichiometric H2 /O2 /N2 /Ar mixture specified earlier. Again, the reaction mechanism used is taken from Li et al. [19]. The FLUENT solutions closely follow the SENKIN solutions in terms of time variations of pressure, temperature, and species. Although a small discrepancy is noted when the time of ignition is approached, the predicted ignition delays of
Fig. 3. Comparison of solutions from FLUENT (line) with SENKIN (circles) for the constant volume adiabatic calculation using the mechanism of Li et al. [19]. Initial pressure = 30 bar, initial temperature = 950 K, and mixture molar composition: H2 /O2 /N2 /Ar = 12.5/6.25/18.125/63.125.
FLUENT (21.07 ms) and SENKIN (22.12 ms) are fairly close. Similar simulations are also conducted for other initial conditions of pressure and temperature, and FLUENT typically predicts an ignition delay 4–5% shorter than the SENKIN solution. Considering the differences in the solvers employed in FLUENT and SENKIN and the precision limitation typically in the CFD package, the FLUENT solution is deemed adequate, especially for the prediction of highly temperature-sensitive ignition delay. 4.2. Nonreactive RCM simulations FLUENT results for nonreactive compression in RCM are obtained here by simply turning off the chemistry during calculations, and are first compared with experimental pressure traces for nonreactive runs. Note that experimental nonreactive pressure traces are obtained by using a mixture of H2 /N2 /Ar with the same specific heat, initial pressure, and initial temperature as the H2 /O2 /N2 /Ar mixture. Through such a nonreactive comparison, we can assess whether the experimentally observed heat loss effect can be properly captured by the FLUENT simulation. Fig. 4 shows a comparison for P0 = 705.8 Torr and T0 = 298 K. The compression ratio and stroke length are 10.32 and 238.15 mm, respectively. Fig. 4
422
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
demonstrates that the simulation matches closely with the experiment except for some discrepancy near TDC. In particular, the FLUENT simulation predicts a slightly lower pressure near TDC. In simulation, the clearance volume at TDC is taken as slightly larger than the geometric volume to compensate for various dead volumes in experiment. The lower simulated pressure at TDC could be due to overcompensation for dead volumes. Simulation also shows a slightly lower heat loss near TDC than experiment, as reflected by lower rate of pressure drop in simulation. Nevertheless, the overall agreement is pretty good. Fig. 5 shows simulated temperature fields at varying times after TDC (tTDC = 33 ms) with fuel oxidation chemistry being suppressed during calculations. For the ease of representation, only the main reaction chamber apart from crevice is shown in Fig. 5. Furthermore, for clarity only a range of temperatures, as indicated in the legends, is shown. Also, the zone that appears white in the temperature field has temperatures lower than the lowest temperature in the legends. At this PC value of 29.08 bar, some effect
Fig. 4. Comparison of experimental (dotted line) and FLUENT-simulated (solid line) pressure traces. P0 = 705.8 Torr and T0 = 298 K.
of vortex on temperature distribution is seen in Fig. 5. From the FLUENT results, Fig. 6 further plots and compares the evolution of Tmax , Tmavmain , and Tmav , respectively representing the maximum instantaneous temperature, the mass averaged temperature of the main chamber (excluding crevice), and the mass averaged temperature of the entire chamber (including crevice). It is seen from Fig. 6 that near and after the TDC Tmavmain is lower than Tmax due to the effect of cold boundary layer and vortex. Hence, the difference between Tmavmain and Tmax indicates the extent of boundary layer. Moreover, Tmav is significantly lower than Tmavmain due to the crevice’s being at a much lower temperature. During compression, simulations also show that the mass fraction inside crevice increases and reaches 0.14 at TDC, which continues to increase gradually during the postcompression period due to the flow of mass into the crevice.
Fig. 6. Temperature profiles from the nonreactive FLUENT simulation. Tmax = maximum instantaneous temperature, Tmavmain = mass averaged temperature of the main reaction chamber (without crevice), Tmav = mass averaged temperature of the entire reaction chamber (including crevice). P0 = 705.8 Torr, T0 = 298 K, and PC = 29.08 bar. Mixture molar composition: H2 /O2 /N2 /Ar = 12.5/6.25/18.125/63.125.
Fig. 5. Comparison of simulated temperature fields within the main reaction chamber (excluding crevice) at varying times after the end of compression. P0 = 705.8 Torr, T0 = 298 K, and PC = 29.08 bar. The orientation of representation is identical to that in Fig. 1c. Mixture molar composition: H2 /O2 /N2 /Ar = 12.5/6.25/18.125/63.125. In this nonreactive FLUENT simulation, chemical reactions are not considered.
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Based on the above nonreactive FLUENT results, implications for the reactive cases are as follows. Because of the Arrhenius dependence of reaction rate on temperature, the higher temperature core region will experience higher reactivity than the colder regimes such as the boundary layer and crevice zone, leading to nonuniform heat release during reactions. In addition, the heat release within the main chamber could result in an enhanced flow of mass into the crevice. Such issues are not accounted for while using a zerodimensional simulation. The following reactive RCM simulations will help clarify the effect of nonuniform heat release on ignition delay.
Fig. 7. Time variations of Tmax , Tmavmain , and Tmav from the reactive and nonreactive FLUENT simulations. P0 = 705.8 Torr, T0 = 298 K, PC = 28.6 bar, and TC = 964.25 K.
423
4.3. Reactive RCM simulations It has been shown in our previous nonreactive RCM simulations [8] that temperature nonhomogeneity induced by boundary layer and vortex becomes more pronounced under conditions of low compressed pressure. Therefore, reactive RCM simulations are conducted over a range of operating pressures. In the following, two representative cases, PC ∼ 15 and 29 bar, are presented and discussed. 4.3.1. Case 1: PC ∼ 29 bar This simulation is conducted for initial conditions of P0 = 705.8 Torr and T0 = 298 K. The compression ratio, stroke length, and clearance are respectively 10.2, 237.87 mm, and 23.37 mm, yielding compressed conditions of PC = 28.6 bar and TC = 964.25 K. Fig. 7 shows the time variations of Tmax , Tmavmain , and Tmav for both reactive and nonreactive simulations. Half-way through the induction period, the temperature traces for the reactive case begin to deviate from the corresponding nonreactive traces, highlighting the importance of pre-ignition heat release. As time of ignition approaches, at time ∼59 ms, all temperatures for the reactive simulation in Fig. 7 show similar rapid rising patterns. By defining ignition delay as the time interval from the end of compression stroke to the event of ignition, it is ∼26 ms for the condition of Fig. 7. Fig. 8 shows and compares the temperature field and the distribution of OH mole fraction for various
Fig. 8. Temperature fields and OH mole fraction distributions within the main reaction chamber (excluding crevice) at varying times from the reactive FLUENT simulations. P0 = 705.8 Torr, T0 = 298 K, PC = 28.6 bar, and TC = 964.25 K. The orientation of representation is identical to that in Fig. 1c.
424
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Fig. 9. Mass fraction inside the crevice for the reactive and nonreactive FLUENT simulations. P0 = 705.8 Torr, T0 = 298 K, PC = 28.6 bar, and TC = 964.25 K.
times near the point of ignition. For the present case of hydrogen ignition, local ignition is seen to occur at the location of maximum temperature and OH mole fraction. Fig. 9 further presents a plot of the mass fraction that is contained in the crevice for both reactive and nonreactive cases. For the reactive case, due to heat release in the main reaction chamber, there is mass flow into the crevice, and hence mass fraction inside the crevice begins to increase in comparison to the nonreactive case. At the instant of ignition, there is significant mass flow into the crevice. Since our major concern is the effect of this mass flow on altering ignition delay, it is noted that prior to ignition the effect of this mass flow is much less as the mass fraction inside the crevice increases from 0.175 for the nonreactive case to 0.185 for the reactive case at time = 58.9 ms. Figs. 10a and 10b compare the results of SENKIN and FLUENT for the nonreactive and reactive simulations, respectively. For the FLUENT simulation, the maximum instantaneous temperature (Tmax ) is plotted. To obtain the pressure and temperature traces from SENKIN, the nonreactive pressure trace from the FLUENT calculations is first used to back calculate the “effective volume” as described previously. The effective volume is then used for the SENKIN simulation. It is noted that based on this approach not only the pressure but also the maximum instantaneous temperature is also exactly reproduced by the zerodimensional SENKIN solution, as shown in Fig. 10a. For the reactive case, it is seen from Fig. 10b that both calculations show the same features of pressure and temperature variations, and predict practically identical ignition delays; 25.94 ms from FLUENT and 26.564 ms from SENKIN. Therefore, despite the temperature nonuniformity and nonuniform heat release, the ignition event is adequately predicted by the zerodimensional model for this case of PC ∼ 29 bar and ignition delay of ∼26 ms. Since the effect of vortex
Fig. 10. Comparison of pressure and temperature traces from FLUENT (lines) and SENKIN (circles): (a) nonreactive simulation and (b) reactive simulation. For the FLUENT results, maximum instantaneous values are plotted. P0 = 705.8 Torr, T0 = 298 K, PC = 28.6 bar, and TC = 964.25 K.
becomes more pronounced when compressed pressure is lower [8], a case of PC ∼ 15 bar with a longer ignition delay is discussed next. 4.3.2. Case 2: PC ∼ 15 bar The simulation described in this section is conducted for P0 = 371.28 Torr, T0 = 298 K, compression ratio of 10.44, stroke length of 238.46 mm, and clearance of 22.81 mm, yielding compressed conditions of PC = 15.2 bar and TC = 967.6 K. The time variations of Tmax , Tmavmain , and Tmav for the reactive and nonreactive simulations are shown and compared in Fig. 11. With ignition occurring at time = 100.25 ms, ignition delay is 67.25 ms for this case. From time = 50 ms onward, the temperature traces for the reactive case begin to deviate from the nonreactive case. In contrast to the high pressure case (cf. Fig. 7), the features of temperature increase in Fig. 11 are somewhat different. Specifically, from time = 66 ms onward, while the maximum instantaneous temperature keeps on increasing, Tmavmain and Tmav continue to decrease due to significant effect of the cooler boundary layer and crevice zone. In addition, whereas the maximum instantaneous temperature shows a gradual rise in temperature before the event of ignition, the two mass averaged tempera-
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Fig. 11. Time variations of Tmax , Tmavmain , and Tmav from the reactive and nonreactive FLUENT simulations. P0 = 371.28 Torr, T0 = 298 K, PC = 15.2 bar, and TC = 967.6 K.
tures exhibit abrupt temperature rise only at the point of ignition. It is also noted that prior to ignition, the mass fraction inside the crevice increases from 0.202 for the nonreactive case to 0.207 for the reactive case. Fig. 12 shows the distributions of temperature and mole fractions of OH and H2 O2 at varying times. The temperature field for the nonreactive simulation is also included in Fig. 12a as a reference. For this relatively low PC condition, the effect of vortex is noted to be more significant than for Case 1. Comparing reactive and nonreactive temperature fields at time = 99 ms (66 ms after the end of compression), the nonuniform heat release is seen to increase temperature nonuniformity. As with the temperature field, there is similar nonuniformity in the mole fraction distributions of radicals. Figs. 12c and 12d also demonstrate that the peak mole fractions of OH and H2 O2 at TDC (time = 33 ms) are in the range of 10−12 and 10−8 , respectively, which are orders of magnitude smaller than those during the induction period and prior to ignition. As with Case 1, the nonreactive simulation using SENKIN for Case 2 also shows the same level of agreement as noted in Fig. 10a. Fig. 13a compares the pressure and temperature traces of FLUENT and SENKIN for the reactive case. It is seen that the effect of nonuniform heat release plays a role as FLUENT predicts an ignition delay slightly longer than the SENKIN calculations; 67.25 ms for FLUENT and 63 ms for SENKIN. Furthermore, in contrast to Case 1, different behavior of temperature and pressure rise is noted here. In the SENKIN solution, temperature and pressure increase rather gradually before the event of ignition. Whereas in the FLUENT simulation, although the maximum instantaneous temperature increases in the same way as SENKIN, the pressure does not show similar behavior and increases abruptly at the point of ignition. The behavior of abrupt pressure rise could be caused by nonuniform heat release
425
in that it lowers the rate of pressure rise, while the local instantaneous temperature may continue to increase due to local heat release, and eventually lead to rapid ignition. The above results for Case 2 demonstrate that although the ignition delay may still be predicted correctly within a few percent by the zero-dimensional simulation, the features of pressure rise can be different. Thus, caution needs to be exercised in conducting mechanism validation for pre-ignition heat release. Specifically, depending on the degree of nonuniformity, the measured pressure traces in RCM experiments may show a slower rate of rise in pressure near the point of ignition than that obtained by the zerodimensional modeling. Fig. 13b further compares the mole fraction profiles of some key radicals from the FLUENT and SENKIN simulations. For the FLUENT results, the maximum instantaneous mole fractions are plotted in Fig. 13b. Except for the discrepancy in ignition delay, the mole fraction profiles from FLUENT match with the SENKIN solutions closely. Again, the radical mole fractions at the end of the compression (time = 33 ms) are orders of magnitude smaller than those during the chemical induction event. 4.4. Comparison of simulated ignition delays The validity of the SENKIN simulation in predicting ignition delay is further examined by comparing to the FLUENT simulation for varying compressed conditions, covering the variation of ignition delay from a few ms to nearly 100 ms. The results are summarized and compared in Fig. 14. The solid line in Fig. 14 represents the situation of both the FLUENTand SENKIN-simulated ignition delays being equal. The pressures and temperatures indicated in Fig. 14 are the maximum instantaneous values at TDC for each case. The zero-dimensional model based on the approach of adiabatic volume expansion performs quite well for the present investigations of hydrogen ignition, and there is some error involved only when ignition delay is very long and compressed pressure is also low. 4.5. Implications for RCM species sampling Apart from the measurement of global combustion property as ignition delay, RCMs can also be used for acquiring species evolution profiles during the chemical induction period by incorporating a sampling apparatus. Such sampling measurements provide detailed experimental data for rigorous validation and development of chemical kinetic mechanisms. However, because the sampling procedure results in mixing of the reacting mixture from all different zones
426
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Fig. 12. Temperature fields and mole fraction distributions of OH and H2 O2 within the main reaction chamber (excluding crevice) at varying times. P0 = 371.28 Torr, T0 = 298 K, PC = 15.2 bar, and TC = 967.6 K. The orientation of representation is identical to that in Fig. 1c.
in an RCM, including the high-temperature core region, the boundary layer, and the crevice, such mixing can pose a difficulty in interpreting the measured species concentration as an average quantity. As such, it is unclear whether it is appropriate to directly com-
pare the measured species profiles with the simulated results using a zero-dimensional model. The present FLUENT simulations can then be used to assess the influence of such mixing on the RCM sampling experiments.
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
Fig. 13. Comparison of the time-varying results from FLUENT (lines) and SENKIN (circles): (a) pressure and temperature; (b) mole fractions of key radicals. For the FLUENT results, maximum instantaneous values are plotted. P0 = 371.28 Torr, T0 = 298 K, PC = 15.2 bar, and TC = 967.6 K.
Fig. 14. Comparison of ignition delays obtained from FLUENT and SENKIN under various compressed conditions of PC and TC .
Fig. 15 compares the profiles of the mass-weighted average mass fractions of HO2 , H2 O2 , and H2 O obtained from FLUENT with the SENKIN results for two cases of PC = 15.2 and 32.52 bar. For both compressed pressures, ignition delays are in the range of 63 ms. For PC = 15.2 bar, although Fig. 13 shows that the maximum instantaneous mole fractions match closely between FLUENT and SENKIN, it is seen from Fig. 15a that the average mass fractions from FLUENT are generally lower than those from the
427
Fig. 15. Comparison of mass-weighted average mass fractions of HO2 , H2 O2 , and H2 O from FLUENT (triangles) and SENKIN (circles). (a) P0 = 371.28 Torr, T0 = 298 K, PC = 15.2 bar, and TC = 967.6 K; (b) P0 = 866.3 Torr, T0 = 298 K, PC = 32.52 bar, and TC = 939 K.
SENKIN solutions. Fig. 15a also demonstrates that the ratio of mass fractions of H2 O, HO2 or H2 O2 from SENKIN with the average mass fractions from FLUENT is ∼1.2 near TDC (time = 33 ms). This difference increases with postcompression time due to the increased influence of boundary layer and crevice, with the ratio becoming ∼3 near the ignition event. On the other hand, for PC = 32.52 bar shown in Fig. 15b, this ratio is ∼1.2 at TDC and increases slightly to become ∼1.5 prior to ignition. The extent of nonuniform species concentration is less at higher pressure due to the reduced effects of heat loss and vortex. These observations have significant implications for the species sampling measurements in an RCM. It is clear that RCM sampling experiments are expected to provide data of better quality for conditions of elevated pressures. Furthermore, when a zero-dimensional modeling is used to simulate RCM sampling experiments for the purpose of mechanism validation, proper interpretation of the measured species profiles is very essential. 5. Conclusions Because CFD simulations of RCM experiments with detailed chemistry can be computationally pro-
428
G. Mittal et al. / Combustion and Flame 155 (2008) 417–428
hibitive, for the purpose of validating detailed kinetic mechanisms of various fuels, the use of a zero-dimensional modeling approach is practically inevitable. In this work, FLUENT simulations including detailed chemistry are conducted for modeling hydrogen ignition in an RCM with a creviced piston in order to assess the validity of zero-dimensional modeling in ignition delay prediction. The key results of this computational study on hydrogen ignition are as follows: 1. Based on the comparison of FLUENT simulations and SENKIN calculations, it is demonstrated that the zero-dimensional model based on the approach of adiabatic volume expansion generally performs well in adequately capturing the time variations of pressure and maximum instantaneous temperature for the nonreactive cases and correctly predicting the hydrogen ignition delays for the reactive cases. 2. When the effect of temperature nonuniformity becomes significant, particularly for the conditions of low compressed pressures and long ignition delay times, the zero-dimensional model may predict hydrogen ignition delays somewhat shorter than the CFD modeling. However, under such conditions of longer ignition delays, the resulting uncertainty is of few percent. 3. For the hydrogen cases investigated, zerodimensional modeling along with the adiabatic volume expansion approach would be still considered acceptable in terms of mechanism validation for such global combustion property as ignition delay. 4. For species profile validation, on the other hand, there are issues involved with the use of zerodimensional model in that sampling measurements in an RCM are average species concentration levels at a given instant, which can be significantly different from the maximum instantaneous values typically calculated by the zero-dimensional simulations. Therefore, caution should be exercised in interpreting the measured RCM species data and comparing them with the results obtained using a zero-dimensional model. 5. Furthermore, the present FLUENT results indicate that the sampling experiments are expected to provide better data under elevated pressure conditions. Finally, although the conclusions outlined above are expected to be valid for the ignition cases of other hydrocarbons exhibiting single-stage ignition characteristics, they may not be applicable to the conditions of two-stage ignition phenomena. Further CFD studies for hydrocarbon two-stage ignition are required to ascertain the validity of the present zero-dimensional modeling approach for such cases.
Acknowledgments This work was supported by the National Aeronautics and Space Administration under Grant NNX07AB36Z, with the technical monitoring of Dr. Krishna P. Kundu.
References [1] J.F. Griffiths, Q. Jiao, A. Schreiber, J. Meyer, K.F. Knoche, W. Kardylewski, Combust. Flame 93 (1993) 303–315. [2] D. Lee, S. Hochgreb, Combust. Flame 114 (1998) 531– 545. [3] J. Clarkson, J.F. Griffiths, J.P. Macnamara, B.J. Whitaker, Combust. Flame 125 (2001) 1162–1175. [4] J.F. Griffiths, J.P. MacNamara, C. Mohamed, B.J. Whitaker, J. Pan, C.G.W. Sheppard, Faraday Discuss. 119 (2001) 287–303. [5] J. Würmel, J.M. Simmie, Combust. Flame 41 (2005) 417–430. [6] P. Park, J.C. Keck, SAE paper 9000027, 1990. [7] G. Mittal, C.J. Sung, Combust. Flame 145 (2006) 160– 180. [8] G. Mittal, M.P. Raju, C.J. Sung, K. Kundu, submitted for publication. [9] A.E. Lutz, R.J. Kee, J.A. Miller, SENKIN: A Fortran Program for Predicting Homogeneous Gas Phase Chemical Kinetics with Sensitivity Analysis, Report No. SAND 87-8248, Sandia National Laboratories, 1988. [10] C.K. Westbrook, W.J. Pitz, J.E. Boercker, H.J. Curran, J.F. Griffiths, C. Mohamed, M. Ribaucour, Proc. Combust. Inst. 29 (2002) 1311–1318. [11] S. Tanaka, F. Ayala, J.C. Keck, Combust. Flame 133 (2003) 467–481. [12] M. Ribaucour, R. Minetti, L.R. Sochet, H.J. Curran, W.J. Pitz, C.K. Westbrook, Proc. Combust. Inst. 28 (2000) 1671–1678. [13] G. Mittal, C.J. Sung, R.A. Yetter, Int. J. Chem. Kinet. 38 (2006) 516–529. [14] G. Mittal, C.J. Sung, Combust. Sci. Technol. 179 (3) (2007) 497–530. [15] G. Mittal, C.J. Sung, Combust. Flame 150 (2007) 355– 368. [16] FLUENT 6.1 Documentation, Fluent Inc. [17] R.J. Kee, F.M. Rupley, J.A. Miller, Chemkin-II: A Fortran Chemical Kinetics Package for the Analysis of Gas Phase Chemical Kinetics, Report No. SAND89-8009B, Sandia National Laboratories, 1989. [18] R.J. Kee, G.D. Lewis, J. Warnatz, M.E. Coltrin, J.A. Miller, A Fortran Computer Code Package for the Evaluation of Gas-Phase, Multicomponent Transport Properties, Report No. SAND 86-8246, Sandia National Laboratories, 1986. [19] J. Li, Z. Zhao, A. Kazakov, F.L. Dryer, Int. J. Chem. Kinet. 36 (2004) 566–575.