- Email: [email protected]
Predicting chemical species in spark-ignition engines
Energy 29 (2004) 449–465 www.elsevier.com/locate/energy
Predicting chemical species in spark-ignition engines Emel Evren Selamet a,, Ahmet Selamet a, James M. Novak b a
Department of Mechanical Engineering, The Ohio State University, Columbus, OH 43210, USA b Ford Motor Company, Powertrain Operations, Dearborn, MI 48121, USA
Abstract The unsteady motion of chemical species in the intake and exhaust ducts of spark-ignition internal combustion (I.C.) engines is studied numerically by employing a finite-difference-based engine simulation code. The time-dependent mass fractions of six dominant products (CO2, H2O, CO, H2, O2, N2) combined with five additional minor species (H, NO, O, OH, N) are computed in the combustion chamber and tracked throughout the breathing system for wide-open-throttle as well as part-load operating conditions. The property calculations employ the NASA database to determine the composition. The effect of the number of product species on the engine performance and dynamic quantities, including pressure and temperature in the breathing system, is also investigated. # 2003 Elsevier Ltd. All rights reserved.
1. Introduction The ability to accurately predict the temporal and spatial distribution of chemical species in the breathing system of spark-ignition (SI) engines is critical for an improved engine design in terms of performance and reduced pollutant and noise emissions. The relevance of such information is described as applied specifically to the induction, in-cylinder, and exhaust systems. In connection with the time-resolved measurements of hydrocarbon concentrations along the intake port with a fast-response flame ionization detector, Cheng et al. [1] discuss the two typical backflows of cylinder contents into the intake ports and manifold: (1) overlap backflow due to significant vacuum generated between the throttle body and the intake valves during the valve overlap period of part-load operation, and (2) displacement backflow at low speeds due to the piston moving some of the cylinder contents, including fresh charge and residuals back into
Corresponding author. Tel.: +1-614-247-7298; fax: +1-614-688-4111. E-mail address: [email protected] (E.E. Selamet).
0360-5442/$ - see front matter# 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.energy.2003.10.006
450
E.E. Selamet et al. / Energy 29 (2004) 449–465
the induction system at the beginning of compression stroke. Since the products in the former category need to be pushed back first into the cylinder rather than the fresh mixture and the cylinder loses in the latter category part of its fresh charge back into the induction system, both backflows are detrimental to the volumetric efficiency, and therefore to the torque of the engine. However, Cheng and coworkers further discuss the importance of the backflow on the mixture preparation in SI engines with port fuel injection, particularly at part-load conditions. The hot combustion products due to either type of backflow modify the thermal environment in the intake ports by elevating gas temperature and assist the fuel evaporation. Thus, the ability to determine species distribution is important for the induction system. Clearly, the source of indicated work, therefore power, transferred from the gases in the combustion chamber to the piston is the enthalpy of combustion. This energy is released over the burn duration starting late in the compression stroke following the burn delay and ending not too far from the top-center into the expansion stroke. As expected from the first law of thermodynamics, the released energy is a function of changing mixture composition, which gradually transforms during the combustion from fresh mixture and residuals into products. In addition to determining the indicated power, the species concentrations, as well as their temperature and pressure, will dictate the engine-out emissions of various pollutants, including carbon monoxide and nitrogen oxides (NOx). The behavior of the catalytic converter in the exhaust system depends on the instantaneous mixture composition incident on the coated substrate. External or internal exhaust gas recirculation (EGR) is a contemporary technique used to reduce the peak flame temperatures in the combustion chamber, thereby significantly reducing the amount of NOx [2–4]. Thus, to investigate either the light-off and conversion behavior of catalytic converters, and/or the impact of EGR on engine performance, combustion stability, and emissions, the chemical species need to be known. Furthermore, since the instantaneous local values of gas temperature will be influenced by the composition, the speed of sound, which is proportional to the square-root of temperature, will also be affected, influencing in turn, the overall speed of wave propagation, the tuning of wave dynamics, and noise propagation and cancellation in the breathing system of engines. The time-domain approach to treat the balance equations (with the exception of species equations) in the breathing system of engines was first proposed by Benson et al. [5] who employed the method of characteristics. Since then, a number of more efficient numerical techniques have been developed, including the finite difference approach of Chapman et al. [6] and CE–SE (conservation element–solution element) method applied by Onorati et al. [7,8]. The transport of chemical species with time domain techniques has been studied, for example, by Baruah [9], Pearson and Winterbone [10], Zhang and Widener [11], and Onorati et al. [7,8]. By employing an extension of the finite-difference method of Chapman et al. [6], the objective of the present study is to incorporate the species equations into the balance equations for unsteady intake and exhaust flows and investigate the temporal and spatial variation of species in the breathing system of engines operating under both full- and part-load conditions. Following Introduction, Section 2 briefly describes the balance equations and the numerical approach used primarily for ducts, and Section 3 performs equilibrium calculations for the combustion of isooctane-air. Section 4 predicts the species distribution in the breathing system of a
E.E. Selamet et al. / Energy 29 (2004) 449–465
451
contemporary six-cylinder engine, and discusses the findings. Section 5 concludes the study with final remarks.
2. Computational approach For one-dimensional, unsteady, and compressible flow in ducts of variable cross-section with neglected axial conduction, the balance equations for mass, momentum, and energy may be written as follows. Mass: @ @ ðqAÞ ¼ ðqAUÞ; @t @x
(1)
momentum: @ @ @ ðqAUÞ ¼ ðqAU 2 Þ ðpAÞ sw P; @t @x @x energy: @ @ @ðUAÞ ðqAeÞ ¼ ðqAUeÞ p þ sw PU qP; @t @x @x
(2)
(3)
where q is the density, A the cross-sectional area, U the velocity, p the pressure, P the perimeter, e the specific internal energy, 1 fF qU 2 2 is the wall shear stress, and sw ¼
(4)
q ¼ ht ðT T1 Þ
(5)
is the heat transfer from control volume to the surroundings, fF and ht being the Fanning friction and heat transfer coefficients, respectively, and T the temperature. The model implements the Reynolds analogy between these two coefficients (see, for example, Arpaci et al. [12]). The species in ducts are determined by solving the species continuity equations: @ @ j ¼ 1; 2; :::; N 1 (6) qAxj ¼ qAUxj ; @t @x combined with N X xj ¼ 1; j¼1
where xj is the mass fraction of species j and N the number of species considered. The diffusion and chemical reactions are neglected in Eq. (6). Introducing the equation of state, p ¼ qRT, R being the gas constant, and the internal energy as a function of temperature, pressure, and species mass fractions, the energy equation may be
452
E.E. Selamet et al. / Energy 29 (2004) 449–465
expressed in terms of temperature, pressure, gas constants, species mass fractions, and their derivatives as illustrated in Appendix A based on an approach by Olikara and Borman [13]. Eq. (A.4) in Appendix A yields the gas temperature both in ducts and combustion chamber from computed internal energy with changing composition. A single-zone model is used for the combustion of homogeneous mixture in the cylinder along with a simple sin2 burn curve [6]. The mass fraction burned is then determined by: dxb n ð h hb Þ 2 sin p ; ¼ 12 Dh Dh dt where n [rpm] is the engine speed, h the instantaneous crank angle, hb the crank angle at the initiation of burn, and Dh the burn duration. Properties such as specific heats and enthalpies for species are obtained from NASA database [14]. The property of the mixture is then determined by: propmix ¼
N X
xj propj
j¼1
to account for the effect of gas composition variation in combustion chambers and ducts. Eqs. (1)–(3) and (6) are discretized by employing the finite difference method of Chapman et al. [6] discussed by Selamet et al. [15], which may be referred to for a more detailed description. The technique is explicit; therefore, the selection of the time step and the grid size needs to satisfy the Courant, Friedrichs, and Lewy (CFL) stability criterion. The staggered mesh used in the discretization divides a duct into cells (differential control volumes) with vector quantities located at nodal points, and scalar quantities at cell midpoints. At locations where the flow is smooth, the computational approach is based on a second order numerical scheme. However, in regions where solution changes rapidly to cause nonphysical oscillations in a higher-order algorithm, such as a shock wavefront, the technique reverts to a first-order approximation to suppress this spurious behavior.
3. Equilibrium calculation for products The computation of equilibrium composition is based on the Olikara and Borman [13] approach. Similar to their work (with the exception of Ar which is excluded from the present study), a total of 11 product species are considered: CO2, H2O, CO, H2, O2, N2, H, NO, O, OH, and N. The equilibrium composition depends on the fuel, oxidizer, equivalence ratio, temperature, and pressure. Fig. 1 compares the product equilibrium mass fractions for rich (fuel/air equivalence ratio, / ¼ 1:1, a typical value for wide-open-throttle (WOT) operation) combustion of isooctane (C8H18) and air for two selected pressures, p ¼ 1 and p ¼ 4 MPa. The reason for the choice of these pressures will become clear in the next section. The effect of pressure becomes noticeable above 2000 K: xCO2 , xH2 O , and xN2 increase with pressure in this range in contrast to the rest of the species. At higher temperatures (T > 2500 K), xCO2 and xH2 O decrease with increasing T, while their behavior differs at lower temperatures. The mass fractions of O2, H, NO, O, OH, and N increase with increasing T for T > 2250 K including xCO which exhibits
E.E. Selamet et al. / Energy 29 (2004) 449–465
453
Fig. 1. Mass fraction of products for isooctane þ air combustion with / ¼ 1:1 as a function of temperature at P ¼ 1 and P ¼ 4 MPa.
454
E.E. Selamet et al. / Energy 29 (2004) 449–465
a more gradual increase at lower temperatures, while xH2 goes through a minimum near T ¼ 2400 K. The reader may be referred to Ibrahim et al. [16] for the product equilibrium composition resulting from the combustion of other fuels such as ethanol–octane and methanol– octane with air. While the foregoing computations include CO and NO in the equilibrium composition, the actual amount of these pollutants in internal combustion engines may significantly deviate from these estimates, since the available time for physical process (burn duration) becomes comparable or smaller than the time needed to reach chemical equilibrium. The models then have to incorporate the time-dependency of CO and NO reactions, and therefore the chemical kinetics. The detailed kinetics mechanism has been excluded from the present study for simplicity, although the simulation package has the ability to account for such phenomenon. For example, Miller et al. [17] presented a super-extended Zeldovich mechanism including 13 species and 67 kinetics reactions predominantly based on the work of others, including the classical study by Miller and Bowman [18]. Such predictive capability can then be used to investigate the tradeoffs between improved fuel economy and reduced NOx emissions. It should be stressed that, with the equilibrium arguments employed in the present study, the CO amounts are expected to be underestimated primarily due to the time-dependency of CO-freezing.
4. Results with engine simulation The computations are performed for a Ford 3.0L V6 four-stroke spark-ignition engine operating first wide-open-throttle at a mid-speed of n ¼ 3000 rpm burning isooctane and air mixture with an equivalence ratio of / ¼ 1:1 (fuel-rich). The basic features of this engine are summarized in Table 1 ([19]), where valve events are given for cylinder #1 with reference to the top-center (TC) at the end of compression stroke (ATC ¼ After TC). A schematic representation of the breathing system is given in Fig. 2 where large circles with C# designate cylinders (C1–C6), small circles the junctions (1–6, 9–13), and the two squares (8 and 14) the ambient junctions (interfaces). The ducts connecting these cylinders and junctions are also numbered from 1 to 24, with 1–13 designating the intake side depicted by solid lines and 14–24 the exhaust side by dashed lines. Junctions 1–8 represent the induction system, while 9–14 the exhaust system. The Table 1 Engine specifications Number of cylinders Bore Stroke Compression ratio Firing order Number of intake or exhaust valves per cylinder Intake valve open Intake valve close Exhaust valve open Exhaust valve close
6 8.9 cm 8.0 cm 9.25 142536 1 334.5 ATC 622.5 ATC 105.5 ATC 393.5 ATC
E.E. Selamet et al. / Energy 29 (2004) 449–465
455
Fig. 2. A schematic representation of the breathing system of a Ford 3.0L V6 engine.
throttle body is located upstream of junction #1. The results in the present work will be illustrated in terms of cylinder #1 (designated by C1 in Fig. 2) and the connecting intake (#1) and exhaust (#14) ducts. The products are initially assumed to be CO2, H2O, CO, H2, O2, N2. Fig. 3 depicts the composition of gas mixture with the exception of xN2 within cylinder #1 throughout the engine cycle: the mass fractions of CO2, H2O, CO increase as a result of the combustion followed by
456
E.E. Selamet et al. / Energy 29 (2004) 449–465
Fig. 3. Mass fraction of in-cylinder species as a function of crank angle.
slight increase in xCO2 and comparable decreases in xCO and xH2 O as the temperature is reduced during the expansion and exhaust strokes. The behavior of mass fractions of these three species with decreasing temperature resembles that of Fig. 1(a–c). After the intake valves are opened, the mass fractions of the products decrease and xO2 increases rapidly due to incoming fresh air. The fuel is injected into the cylinder after intake valves are closed. Due to its rather small magnitude, the mass fraction of H2 is excluded from Fig. 3 and the remainder of the study. Fig. 4 shows the exhaust and intake mass flow rates at the port. The negative and positive values represent flow into and out of the cylinder, respectively. Fig. 5 depicts the species mass fractions and temperature in the first control volume outside the intake valve (designated by I1 in Fig. 2) as a function of crank angle throughout the engine cycle. The results in Fig. 5 can be interpreted better in view of Fig. 4. The product species mass fractions and the temperature increase during the overlap period (when both intake and exhaust valves are open), which lasts, in view of Table 1, 393:5 334:5 ¼ 59 crank angle degrees (CAD) for the present engine. At the beginning of the compression stroke and before the intake valve is closed, the displacement backflow in Fig. 4 (near h ¼ 540 CAD) brings residuals from the cylinder back into the intake as observed in Fig. 5. Fig. 6 shows the species mass fractions as a function of location in duct #1 (recall Fig. 2) at three different times, at h ¼ 360, 390, and 420 CAD. The mass fractions of
Fig. 4. Mass flow rates at the ports as a function of crank angle.
E.E. Selamet et al. / Energy 29 (2004) 449–465
457
Fig. 5. Temperature and mass fraction of species in intake port as a function of crank angle.
CO2, H2O, and CO decrease with distance from the inlet port. The peak product mass fractions are observed at about the middle of the overlap period (approximately h ¼ 360 CAD). Fig. 7(a,b) shows the temperature and mass fractions as a function of crank angle in the first control volume immediately downstream of the exhaust valve, designated by E1 in Fig. 2. The
Fig. 6. Mass fraction of species in intake as a function of distance from intake valve.
458
E.E. Selamet et al. / Energy 29 (2004) 449–465
Fig. 7. Temperature and mass fraction of species in exhaust port as a function of crank angle at location E1.
temperature increases sharply when the exhaust valve is opened and xH2 O and xCO increase slightly due to water–gas reaction resulting in a comparable decrease in xCO2 . Due to its negligible variation, xN2 was excluded from this figure. The trends are similar with the exception of time lag further downstream at location E2. Similar calculations have also been performed with an alternative fuel, methane, for the same equivalence ratio of / ¼ 1:1 and speed of n ¼ 3000 rpm. Figs. 8 and 9 compare the species mass fractions in the combustion chamber and the exhaust port, respectively. As equilibrium arguments at constant T and p would illustrate, methane yields higher H2O and lower CO2, an
Fig. 8. Mass fraction of in-cylinder species as a function of crank angle for isooctane and methane.
E.E. Selamet et al. / Energy 29 (2004) 449–465
459
Fig. 9. Mass fraction of species in exhaust port as a function of crank angle for isooctane and methane.
expected result due to higher H/C ratio. The rest of the observations with methane will not be discussed here however, primarily because of the additional influence of gaseous fuel on the air breathing efficiency in the induction system. Hereafter, the study continues the presentation in terms of isooctane alone, while the current structure also allows the use of a number of other fuels, including methanol, hydrogen, indolene, and European gasoline. The species mass fractions and temperature are examined next for the same engine and speed, now operating at part-load. The degree of throttling is illustrated in terms of the variation in pressure in junction #1 (which is part of the intake manifold) from WOT to part-load operation in Fig. 10. Fig. 11 compares the resulting predictions of in-cylinder pressure versus volume for WOT and part-load. Note that the peak pressures observed in this figure have been rounded to 1 and 4 MPa and used in the equilibrium calculations of Fig. 1 of the preceding section. The restricted air flow in throttled operation and the resulting lower trapped mass and reduced energy release lead to significantly lower pressure during the power loop. As expected, due to the increasing vacuum in the intake manifold as demonstrated in Fig. 10, the amount of backflow increases during the overlap period, as shown in Fig. 12. The behavior in Fig. 12 can be interpreted better in view of Figs. 13 and 14 which compare the in-cylinder pressure, respectively, with the intake (location I1) and exhaust (location E1) during the corresponding valve-
Fig. 10. Pressure in intake junction 1 as a function of crank angle.
460
E.E. Selamet et al. / Energy 29 (2004) 449–465
Fig. 11. In-cylinder pressure versus volume.
Fig. 12. Mass flow rates at ports as a function of crank angle.
Fig. 13. In-cylinder and intake port pressure versus crank angle.
open periods (Recall Table 1). For example, overlap and displacement backflows can be readily anticipated from Fig. 13, which shows higher cylinder pressure than the intake runner at the beginning and end of the intake process. Figs. 15 and 16 compare, respectively, the species mass fractions and temperature in the intake port at location I1 for part-load and wide-open throttle
E.E. Selamet et al. / Energy 29 (2004) 449–465
461
Fig. 14. In-cylinder and exhaust port pressure versus crank angle.
Fig. 15. Mass fraction of species in intake port versus crank angle.
Fig. 16. Temperature in intake port versus crank angle.
operating conditions. The product mass fractions and temperature are significantly higher during the overlap with the throttled engine. The resulting spatial distribution of species for three selected CADs is given in Fig. 17. Note that with the throttled operation products can penetrate
462
E.E. Selamet et al. / Energy 29 (2004) 449–465
Fig. 17. Mass fraction of species in intake as a function of distance from intake valve for throttled engine.
further upstream in the primary runners towards the plenum. Such a distribution will further be modified in the presence of charge motion control valves in the intake primary runners, which will influence the idle stability and exhaust gas recirculation tolerance of the engine. The temperature in the exhaust port at location E1 is about 15% lower during the early stages of blowdown (Fig. 18), but tends to stay higher in rest of the cycle. Next, in addition to the foregoing six dominant species in the combustion products, five minor species H, NO, O, OH, and N are included in the computations (for an example to a further increase to 20 product species, the reader is referred to Zacharias [20]). With products rep-
Fig. 18. Temperature in exhaust port versus crank angle.
E.E. Selamet et al. / Energy 29 (2004) 449–465
463
resented by eleven species, the volumetric efficiency and torque predictions are rather close to the case investigated already with six species. The peak in-cylinder temperature decreases only by 0.7%. The effect of additional five minor species on the exhaust temperature and pressure is found to be negligible.
5. Concluding remarks The temporal and spatial variation of mass fractions of chemical species in the breathing system of SI engines is predicted by employing a finite-difference-based engine simulation code. The temporal variation of mass fractions of six dominant products (CO2, H2O, CO, H2, O2, N2), combined with five additional minor species (H, NO, O, OH, N), is computed in the combustion chamber, as well as in rest of the breathing system for wide-open-throttle and partload operating conditions. Primarily, isooctane is used as the fuel in engine simulation. The effect of number of product species is also quantified on the engine performance and dynamic quantities, including pressure and temperature in the breathing system. In addition to a number of reasons described in Introduction motivating the present study, the ability to account for species will further promote the incorporation of detailed chemical kinetics for the prediction of major pollutants under different operating conditions, including cold-start and idle, thereby providing an analytical design tool to assess the impact on environment. Such time-dependency in chemical kinetics will lead to more realistic predictions of particularly CO and NOx that are known to freeze during the expansion stroke rather than reach an equilibrium. Future refinements should also include more realistic fuel injection physics and incorporate an unburned hydrocarbon mechanism. The resulting capability is then expected to facilitate the evaluation of contemporary concepts such as variable cam timing, variable valve timing, charge motion control valves, and throttless and camless engines.
Acknowledgements This work is supported by FORD Motor Company, Dearborn, MI, USA.
Appendix A. Temperature versus specific energy Let the specific energy, e, be a function of T, p and xi (mass fraction of species): e ¼ f ðT; p; xi Þ: The derivative of e with respect to time then becomes X @e de @e _ @e ¼ p_ þ x_ i ; Tþ dt @T @p @xi i
(A.1)
(A.2)
@e _ where ð_ Þ implies d( )/dt. Note that the last term in Eq. (A.2) is replaced by @/ / in the Olikara and Borman [13], where / is the equivalence ratio. The derivative of total internal energy with
464
E.E. Selamet et al. / Energy 29 (2004) 449–465
respect to time may readily be written in terms of specific internal energy as: dE de dm ¼m þe : dt dt dt
(A.3)
Combining Eqs. (3), (A.2) and (A.3) along with the derivative p_ from ideal gas relationship ð p ¼ qRT Þ yields: P _ ðdE=dtÞ=m B m =m Ci x_ i þ DV_ =V T_ ¼ (A.4) @e=@T þ DF =T for the time derivative of temperature. Here B ¼ e þ D; Ci ¼
@e D @R þ ; @xi R @xi
@e =½1 ðp=RÞð@R=@pÞ ; @p F ¼ 1 þ ðT=RÞð@R=@TÞ; D¼p
where R is the gas constant.
References [1] Cheng C, Cheng WK, Heywood JB, Maroteaux D, Collings N. Intake port phenomena in a spark ignition engine at part load. SAE Paper 912401. Society of Automotive Engineers, 1991. [2] Blumberg PN, Kummer JT. Predictions of NO formation in spark-ignited engines—an analysis of methods of control. Combust Sci Tech 1971;4:73–95. [3] Heywood JB, Higgins JM, Watts PA, Tabaczynski RJ. Development and use of a cycle simulation to predict SI engine efficiency and NOx emissions. SAE Paper 790291. Society of Automotive Engineers, 1979. [4] Heywood JB. Internal combustion engine fundamentals. New York: McGraw-Hill; 1988. [5] Benson RS, Garg RD, Woollatt DA. Numerical solution of unsteady flow problems. Int J Mech Sci 1964;6(1):117–44. [6] Chapman M, Novak JM, Stein RA. Numerical modeling of inlet and exhaust flows in multi-cylinder internal combustion engines. In: Uzkan Teoman, editor. Flows in internal combustion engines. ASME winter annual meeting. Phoenix, Arizona: American Society of Mechanical Engineers; 1982. [7] Onorati A, Ferrari G. Modeling of 1-D unsteady flows in I.C. engine pipe systems: numerical methods and transport of chemical species. SAE Paper 980782. Society of Automotive Engineers, 1998. [8] Onorati A, Ferrari G, D’errico G. Fluid dynamic modeling of the gas flow with chemical specie transport through the exhaust manifold of a four cylinder SI engine. SAE Paper 1999-01-0557. SAE SP-1451 SI engine modeling. Society of Automotive Engineers, 1999. [9] Baruah PC. Numerical solution of non-steady flows with variable specific heats and chemical reactions (including use of emission control devices). In: Horlock JH, Winterbone DE, editors. The thermodynamics and gas dynamics of internal-combustion engines, vol. II. New York: Oxford University Press; 1986, p. 869–902. [10] Pearson RJ, Winterbone DE. Calculating the effects of variations in composition on wave propagation in gases. Int J Mech Sci 1993;35(6):517–37. [11] Zhang H, Widener SK. An integrated engine cycle simulation model with species tracking in piping system. SAE Paper 960077. Society of Automotive Engineers, 1996. [12] Arpaci VS, Selamet A, Kao S. Introduction to heat transfer, 1st ed. 2nd print. New Jersey: Prentice-Hall; 2000.
E.E. Selamet et al. / Energy 29 (2004) 449–465
465
[13] Olikara C, Borman GL. A computer program for calculating properties of equilibrium combustion products with some applications to I.C. engines. SAE Paper 750468. Society of Automotive Engineers, 1975. [14] McBride BJ, Gordon S. Computer program for calculation of complex chemical equilibrium compositions and applications. II. Users manual and program description, NASA reference publication 1311; Also, PAC91 by McBride BJ: Properties and coefficients 1991, Cosmic Program # LEW-15779, 1992. [15] Selamet A, Dickey NS, Novak JM. The Herschel-Quincke tube: a theoretical, computational, and experimental investigation. J Acoust Soc Am 1994;96(5):3177–85. [16] Ibrahim AM, Abdel-Rahim YM, Ibrahim IG. Thermochemical properties of combustion gases of ethanol-octane and methanol-octane blends. SAE Paper 911311. Society of Automotive Engineers, 1991. [17] Miller R, Davis G, Lavoie G, Newman C, Gardner T. A super-extended Zel’dovich mechanism for NOx modeling and engine calibration. SAE Paper 980781. Society of Automotive Engineers, 1998. [18] Miller JA, Bowman CT. Mechanism and modeling of nitrogen chemistry in combustion. Prog Energy Combust Sci 1989;15:287–338. [19] Selamet A, Kothamasu V, Novak JM, Kach RA. Experimental investigation of in-duct insertion loss of catalysts in internal combustion engines. Appl Acoust 2000;60:451–87. [20] Zacharias F. Analytical representation of the thermodynamic properties of combustion gases. SAE Paper 670930. Society of Automotive Engineers, 1967.
Predicting chemical species in spark-ignition engines Emel Evren Selamet a,, Ahmet Selamet a, James M. Novak b a
Department of Mechanical Engineering, The Ohio State University, Columbus, OH 43210, USA b Ford Motor Company, Powertrain Operations, Dearborn, MI 48121, USA
Abstract The unsteady motion of chemical species in the intake and exhaust ducts of spark-ignition internal combustion (I.C.) engines is studied numerically by employing a finite-difference-based engine simulation code. The time-dependent mass fractions of six dominant products (CO2, H2O, CO, H2, O2, N2) combined with five additional minor species (H, NO, O, OH, N) are computed in the combustion chamber and tracked throughout the breathing system for wide-open-throttle as well as part-load operating conditions. The property calculations employ the NASA database to determine the composition. The effect of the number of product species on the engine performance and dynamic quantities, including pressure and temperature in the breathing system, is also investigated. # 2003 Elsevier Ltd. All rights reserved.
1. Introduction The ability to accurately predict the temporal and spatial distribution of chemical species in the breathing system of spark-ignition (SI) engines is critical for an improved engine design in terms of performance and reduced pollutant and noise emissions. The relevance of such information is described as applied specifically to the induction, in-cylinder, and exhaust systems. In connection with the time-resolved measurements of hydrocarbon concentrations along the intake port with a fast-response flame ionization detector, Cheng et al. [1] discuss the two typical backflows of cylinder contents into the intake ports and manifold: (1) overlap backflow due to significant vacuum generated between the throttle body and the intake valves during the valve overlap period of part-load operation, and (2) displacement backflow at low speeds due to the piston moving some of the cylinder contents, including fresh charge and residuals back into
Corresponding author. Tel.: +1-614-247-7298; fax: +1-614-688-4111. E-mail address: [email protected] (E.E. Selamet).
0360-5442/$ - see front matter# 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.energy.2003.10.006
450
E.E. Selamet et al. / Energy 29 (2004) 449–465
the induction system at the beginning of compression stroke. Since the products in the former category need to be pushed back first into the cylinder rather than the fresh mixture and the cylinder loses in the latter category part of its fresh charge back into the induction system, both backflows are detrimental to the volumetric efficiency, and therefore to the torque of the engine. However, Cheng and coworkers further discuss the importance of the backflow on the mixture preparation in SI engines with port fuel injection, particularly at part-load conditions. The hot combustion products due to either type of backflow modify the thermal environment in the intake ports by elevating gas temperature and assist the fuel evaporation. Thus, the ability to determine species distribution is important for the induction system. Clearly, the source of indicated work, therefore power, transferred from the gases in the combustion chamber to the piston is the enthalpy of combustion. This energy is released over the burn duration starting late in the compression stroke following the burn delay and ending not too far from the top-center into the expansion stroke. As expected from the first law of thermodynamics, the released energy is a function of changing mixture composition, which gradually transforms during the combustion from fresh mixture and residuals into products. In addition to determining the indicated power, the species concentrations, as well as their temperature and pressure, will dictate the engine-out emissions of various pollutants, including carbon monoxide and nitrogen oxides (NOx). The behavior of the catalytic converter in the exhaust system depends on the instantaneous mixture composition incident on the coated substrate. External or internal exhaust gas recirculation (EGR) is a contemporary technique used to reduce the peak flame temperatures in the combustion chamber, thereby significantly reducing the amount of NOx [2–4]. Thus, to investigate either the light-off and conversion behavior of catalytic converters, and/or the impact of EGR on engine performance, combustion stability, and emissions, the chemical species need to be known. Furthermore, since the instantaneous local values of gas temperature will be influenced by the composition, the speed of sound, which is proportional to the square-root of temperature, will also be affected, influencing in turn, the overall speed of wave propagation, the tuning of wave dynamics, and noise propagation and cancellation in the breathing system of engines. The time-domain approach to treat the balance equations (with the exception of species equations) in the breathing system of engines was first proposed by Benson et al. [5] who employed the method of characteristics. Since then, a number of more efficient numerical techniques have been developed, including the finite difference approach of Chapman et al. [6] and CE–SE (conservation element–solution element) method applied by Onorati et al. [7,8]. The transport of chemical species with time domain techniques has been studied, for example, by Baruah [9], Pearson and Winterbone [10], Zhang and Widener [11], and Onorati et al. [7,8]. By employing an extension of the finite-difference method of Chapman et al. [6], the objective of the present study is to incorporate the species equations into the balance equations for unsteady intake and exhaust flows and investigate the temporal and spatial variation of species in the breathing system of engines operating under both full- and part-load conditions. Following Introduction, Section 2 briefly describes the balance equations and the numerical approach used primarily for ducts, and Section 3 performs equilibrium calculations for the combustion of isooctane-air. Section 4 predicts the species distribution in the breathing system of a
E.E. Selamet et al. / Energy 29 (2004) 449–465
451
contemporary six-cylinder engine, and discusses the findings. Section 5 concludes the study with final remarks.
2. Computational approach For one-dimensional, unsteady, and compressible flow in ducts of variable cross-section with neglected axial conduction, the balance equations for mass, momentum, and energy may be written as follows. Mass: @ @ ðqAÞ ¼ ðqAUÞ; @t @x
(1)
momentum: @ @ @ ðqAUÞ ¼ ðqAU 2 Þ ðpAÞ sw P; @t @x @x energy: @ @ @ðUAÞ ðqAeÞ ¼ ðqAUeÞ p þ sw PU qP; @t @x @x
(2)
(3)
where q is the density, A the cross-sectional area, U the velocity, p the pressure, P the perimeter, e the specific internal energy, 1 fF qU 2 2 is the wall shear stress, and sw ¼
(4)
q ¼ ht ðT T1 Þ
(5)
is the heat transfer from control volume to the surroundings, fF and ht being the Fanning friction and heat transfer coefficients, respectively, and T the temperature. The model implements the Reynolds analogy between these two coefficients (see, for example, Arpaci et al. [12]). The species in ducts are determined by solving the species continuity equations: @ @ j ¼ 1; 2; :::; N 1 (6) qAxj ¼ qAUxj ; @t @x combined with N X xj ¼ 1; j¼1
where xj is the mass fraction of species j and N the number of species considered. The diffusion and chemical reactions are neglected in Eq. (6). Introducing the equation of state, p ¼ qRT, R being the gas constant, and the internal energy as a function of temperature, pressure, and species mass fractions, the energy equation may be
452
E.E. Selamet et al. / Energy 29 (2004) 449–465
expressed in terms of temperature, pressure, gas constants, species mass fractions, and their derivatives as illustrated in Appendix A based on an approach by Olikara and Borman [13]. Eq. (A.4) in Appendix A yields the gas temperature both in ducts and combustion chamber from computed internal energy with changing composition. A single-zone model is used for the combustion of homogeneous mixture in the cylinder along with a simple sin2 burn curve [6]. The mass fraction burned is then determined by: dxb n ð h hb Þ 2 sin p ; ¼ 12 Dh Dh dt where n [rpm] is the engine speed, h the instantaneous crank angle, hb the crank angle at the initiation of burn, and Dh the burn duration. Properties such as specific heats and enthalpies for species are obtained from NASA database [14]. The property of the mixture is then determined by: propmix ¼
N X
xj propj
j¼1
to account for the effect of gas composition variation in combustion chambers and ducts. Eqs. (1)–(3) and (6) are discretized by employing the finite difference method of Chapman et al. [6] discussed by Selamet et al. [15], which may be referred to for a more detailed description. The technique is explicit; therefore, the selection of the time step and the grid size needs to satisfy the Courant, Friedrichs, and Lewy (CFL) stability criterion. The staggered mesh used in the discretization divides a duct into cells (differential control volumes) with vector quantities located at nodal points, and scalar quantities at cell midpoints. At locations where the flow is smooth, the computational approach is based on a second order numerical scheme. However, in regions where solution changes rapidly to cause nonphysical oscillations in a higher-order algorithm, such as a shock wavefront, the technique reverts to a first-order approximation to suppress this spurious behavior.
3. Equilibrium calculation for products The computation of equilibrium composition is based on the Olikara and Borman [13] approach. Similar to their work (with the exception of Ar which is excluded from the present study), a total of 11 product species are considered: CO2, H2O, CO, H2, O2, N2, H, NO, O, OH, and N. The equilibrium composition depends on the fuel, oxidizer, equivalence ratio, temperature, and pressure. Fig. 1 compares the product equilibrium mass fractions for rich (fuel/air equivalence ratio, / ¼ 1:1, a typical value for wide-open-throttle (WOT) operation) combustion of isooctane (C8H18) and air for two selected pressures, p ¼ 1 and p ¼ 4 MPa. The reason for the choice of these pressures will become clear in the next section. The effect of pressure becomes noticeable above 2000 K: xCO2 , xH2 O , and xN2 increase with pressure in this range in contrast to the rest of the species. At higher temperatures (T > 2500 K), xCO2 and xH2 O decrease with increasing T, while their behavior differs at lower temperatures. The mass fractions of O2, H, NO, O, OH, and N increase with increasing T for T > 2250 K including xCO which exhibits
E.E. Selamet et al. / Energy 29 (2004) 449–465
453
Fig. 1. Mass fraction of products for isooctane þ air combustion with / ¼ 1:1 as a function of temperature at P ¼ 1 and P ¼ 4 MPa.
454
E.E. Selamet et al. / Energy 29 (2004) 449–465
a more gradual increase at lower temperatures, while xH2 goes through a minimum near T ¼ 2400 K. The reader may be referred to Ibrahim et al. [16] for the product equilibrium composition resulting from the combustion of other fuels such as ethanol–octane and methanol– octane with air. While the foregoing computations include CO and NO in the equilibrium composition, the actual amount of these pollutants in internal combustion engines may significantly deviate from these estimates, since the available time for physical process (burn duration) becomes comparable or smaller than the time needed to reach chemical equilibrium. The models then have to incorporate the time-dependency of CO and NO reactions, and therefore the chemical kinetics. The detailed kinetics mechanism has been excluded from the present study for simplicity, although the simulation package has the ability to account for such phenomenon. For example, Miller et al. [17] presented a super-extended Zeldovich mechanism including 13 species and 67 kinetics reactions predominantly based on the work of others, including the classical study by Miller and Bowman [18]. Such predictive capability can then be used to investigate the tradeoffs between improved fuel economy and reduced NOx emissions. It should be stressed that, with the equilibrium arguments employed in the present study, the CO amounts are expected to be underestimated primarily due to the time-dependency of CO-freezing.
4. Results with engine simulation The computations are performed for a Ford 3.0L V6 four-stroke spark-ignition engine operating first wide-open-throttle at a mid-speed of n ¼ 3000 rpm burning isooctane and air mixture with an equivalence ratio of / ¼ 1:1 (fuel-rich). The basic features of this engine are summarized in Table 1 ([19]), where valve events are given for cylinder #1 with reference to the top-center (TC) at the end of compression stroke (ATC ¼ After TC). A schematic representation of the breathing system is given in Fig. 2 where large circles with C# designate cylinders (C1–C6), small circles the junctions (1–6, 9–13), and the two squares (8 and 14) the ambient junctions (interfaces). The ducts connecting these cylinders and junctions are also numbered from 1 to 24, with 1–13 designating the intake side depicted by solid lines and 14–24 the exhaust side by dashed lines. Junctions 1–8 represent the induction system, while 9–14 the exhaust system. The Table 1 Engine specifications Number of cylinders Bore Stroke Compression ratio Firing order Number of intake or exhaust valves per cylinder Intake valve open Intake valve close Exhaust valve open Exhaust valve close
6 8.9 cm 8.0 cm 9.25 142536 1 334.5 ATC 622.5 ATC 105.5 ATC 393.5 ATC
E.E. Selamet et al. / Energy 29 (2004) 449–465
455
Fig. 2. A schematic representation of the breathing system of a Ford 3.0L V6 engine.
throttle body is located upstream of junction #1. The results in the present work will be illustrated in terms of cylinder #1 (designated by C1 in Fig. 2) and the connecting intake (#1) and exhaust (#14) ducts. The products are initially assumed to be CO2, H2O, CO, H2, O2, N2. Fig. 3 depicts the composition of gas mixture with the exception of xN2 within cylinder #1 throughout the engine cycle: the mass fractions of CO2, H2O, CO increase as a result of the combustion followed by
456
E.E. Selamet et al. / Energy 29 (2004) 449–465
Fig. 3. Mass fraction of in-cylinder species as a function of crank angle.
slight increase in xCO2 and comparable decreases in xCO and xH2 O as the temperature is reduced during the expansion and exhaust strokes. The behavior of mass fractions of these three species with decreasing temperature resembles that of Fig. 1(a–c). After the intake valves are opened, the mass fractions of the products decrease and xO2 increases rapidly due to incoming fresh air. The fuel is injected into the cylinder after intake valves are closed. Due to its rather small magnitude, the mass fraction of H2 is excluded from Fig. 3 and the remainder of the study. Fig. 4 shows the exhaust and intake mass flow rates at the port. The negative and positive values represent flow into and out of the cylinder, respectively. Fig. 5 depicts the species mass fractions and temperature in the first control volume outside the intake valve (designated by I1 in Fig. 2) as a function of crank angle throughout the engine cycle. The results in Fig. 5 can be interpreted better in view of Fig. 4. The product species mass fractions and the temperature increase during the overlap period (when both intake and exhaust valves are open), which lasts, in view of Table 1, 393:5 334:5 ¼ 59 crank angle degrees (CAD) for the present engine. At the beginning of the compression stroke and before the intake valve is closed, the displacement backflow in Fig. 4 (near h ¼ 540 CAD) brings residuals from the cylinder back into the intake as observed in Fig. 5. Fig. 6 shows the species mass fractions as a function of location in duct #1 (recall Fig. 2) at three different times, at h ¼ 360, 390, and 420 CAD. The mass fractions of
Fig. 4. Mass flow rates at the ports as a function of crank angle.
E.E. Selamet et al. / Energy 29 (2004) 449–465
457
Fig. 5. Temperature and mass fraction of species in intake port as a function of crank angle.
CO2, H2O, and CO decrease with distance from the inlet port. The peak product mass fractions are observed at about the middle of the overlap period (approximately h ¼ 360 CAD). Fig. 7(a,b) shows the temperature and mass fractions as a function of crank angle in the first control volume immediately downstream of the exhaust valve, designated by E1 in Fig. 2. The
Fig. 6. Mass fraction of species in intake as a function of distance from intake valve.
458
E.E. Selamet et al. / Energy 29 (2004) 449–465
Fig. 7. Temperature and mass fraction of species in exhaust port as a function of crank angle at location E1.
temperature increases sharply when the exhaust valve is opened and xH2 O and xCO increase slightly due to water–gas reaction resulting in a comparable decrease in xCO2 . Due to its negligible variation, xN2 was excluded from this figure. The trends are similar with the exception of time lag further downstream at location E2. Similar calculations have also been performed with an alternative fuel, methane, for the same equivalence ratio of / ¼ 1:1 and speed of n ¼ 3000 rpm. Figs. 8 and 9 compare the species mass fractions in the combustion chamber and the exhaust port, respectively. As equilibrium arguments at constant T and p would illustrate, methane yields higher H2O and lower CO2, an
Fig. 8. Mass fraction of in-cylinder species as a function of crank angle for isooctane and methane.
E.E. Selamet et al. / Energy 29 (2004) 449–465
459
Fig. 9. Mass fraction of species in exhaust port as a function of crank angle for isooctane and methane.
expected result due to higher H/C ratio. The rest of the observations with methane will not be discussed here however, primarily because of the additional influence of gaseous fuel on the air breathing efficiency in the induction system. Hereafter, the study continues the presentation in terms of isooctane alone, while the current structure also allows the use of a number of other fuels, including methanol, hydrogen, indolene, and European gasoline. The species mass fractions and temperature are examined next for the same engine and speed, now operating at part-load. The degree of throttling is illustrated in terms of the variation in pressure in junction #1 (which is part of the intake manifold) from WOT to part-load operation in Fig. 10. Fig. 11 compares the resulting predictions of in-cylinder pressure versus volume for WOT and part-load. Note that the peak pressures observed in this figure have been rounded to 1 and 4 MPa and used in the equilibrium calculations of Fig. 1 of the preceding section. The restricted air flow in throttled operation and the resulting lower trapped mass and reduced energy release lead to significantly lower pressure during the power loop. As expected, due to the increasing vacuum in the intake manifold as demonstrated in Fig. 10, the amount of backflow increases during the overlap period, as shown in Fig. 12. The behavior in Fig. 12 can be interpreted better in view of Figs. 13 and 14 which compare the in-cylinder pressure, respectively, with the intake (location I1) and exhaust (location E1) during the corresponding valve-
Fig. 10. Pressure in intake junction 1 as a function of crank angle.
460
E.E. Selamet et al. / Energy 29 (2004) 449–465
Fig. 11. In-cylinder pressure versus volume.
Fig. 12. Mass flow rates at ports as a function of crank angle.
Fig. 13. In-cylinder and intake port pressure versus crank angle.
open periods (Recall Table 1). For example, overlap and displacement backflows can be readily anticipated from Fig. 13, which shows higher cylinder pressure than the intake runner at the beginning and end of the intake process. Figs. 15 and 16 compare, respectively, the species mass fractions and temperature in the intake port at location I1 for part-load and wide-open throttle
E.E. Selamet et al. / Energy 29 (2004) 449–465
461
Fig. 14. In-cylinder and exhaust port pressure versus crank angle.
Fig. 15. Mass fraction of species in intake port versus crank angle.
Fig. 16. Temperature in intake port versus crank angle.
operating conditions. The product mass fractions and temperature are significantly higher during the overlap with the throttled engine. The resulting spatial distribution of species for three selected CADs is given in Fig. 17. Note that with the throttled operation products can penetrate
462
E.E. Selamet et al. / Energy 29 (2004) 449–465
Fig. 17. Mass fraction of species in intake as a function of distance from intake valve for throttled engine.
further upstream in the primary runners towards the plenum. Such a distribution will further be modified in the presence of charge motion control valves in the intake primary runners, which will influence the idle stability and exhaust gas recirculation tolerance of the engine. The temperature in the exhaust port at location E1 is about 15% lower during the early stages of blowdown (Fig. 18), but tends to stay higher in rest of the cycle. Next, in addition to the foregoing six dominant species in the combustion products, five minor species H, NO, O, OH, and N are included in the computations (for an example to a further increase to 20 product species, the reader is referred to Zacharias [20]). With products rep-
Fig. 18. Temperature in exhaust port versus crank angle.
E.E. Selamet et al. / Energy 29 (2004) 449–465
463
resented by eleven species, the volumetric efficiency and torque predictions are rather close to the case investigated already with six species. The peak in-cylinder temperature decreases only by 0.7%. The effect of additional five minor species on the exhaust temperature and pressure is found to be negligible.
5. Concluding remarks The temporal and spatial variation of mass fractions of chemical species in the breathing system of SI engines is predicted by employing a finite-difference-based engine simulation code. The temporal variation of mass fractions of six dominant products (CO2, H2O, CO, H2, O2, N2), combined with five additional minor species (H, NO, O, OH, N), is computed in the combustion chamber, as well as in rest of the breathing system for wide-open-throttle and partload operating conditions. Primarily, isooctane is used as the fuel in engine simulation. The effect of number of product species is also quantified on the engine performance and dynamic quantities, including pressure and temperature in the breathing system. In addition to a number of reasons described in Introduction motivating the present study, the ability to account for species will further promote the incorporation of detailed chemical kinetics for the prediction of major pollutants under different operating conditions, including cold-start and idle, thereby providing an analytical design tool to assess the impact on environment. Such time-dependency in chemical kinetics will lead to more realistic predictions of particularly CO and NOx that are known to freeze during the expansion stroke rather than reach an equilibrium. Future refinements should also include more realistic fuel injection physics and incorporate an unburned hydrocarbon mechanism. The resulting capability is then expected to facilitate the evaluation of contemporary concepts such as variable cam timing, variable valve timing, charge motion control valves, and throttless and camless engines.
Acknowledgements This work is supported by FORD Motor Company, Dearborn, MI, USA.
Appendix A. Temperature versus specific energy Let the specific energy, e, be a function of T, p and xi (mass fraction of species): e ¼ f ðT; p; xi Þ: The derivative of e with respect to time then becomes X @e de @e _ @e ¼ p_ þ x_ i ; Tþ dt @T @p @xi i
(A.1)
(A.2)
@e _ where ð_ Þ implies d( )/dt. Note that the last term in Eq. (A.2) is replaced by @/ / in the Olikara and Borman [13], where / is the equivalence ratio. The derivative of total internal energy with
464
E.E. Selamet et al. / Energy 29 (2004) 449–465
respect to time may readily be written in terms of specific internal energy as: dE de dm ¼m þe : dt dt dt
(A.3)
Combining Eqs. (3), (A.2) and (A.3) along with the derivative p_ from ideal gas relationship ð p ¼ qRT Þ yields: P _ ðdE=dtÞ=m B m =m Ci x_ i þ DV_ =V T_ ¼ (A.4) @e=@T þ DF =T for the time derivative of temperature. Here B ¼ e þ D; Ci ¼
@e D @R þ ; @xi R @xi
@e =½1 ðp=RÞð@R=@pÞ ; @p F ¼ 1 þ ðT=RÞð@R=@TÞ; D¼p
where R is the gas constant.
References [1] Cheng C, Cheng WK, Heywood JB, Maroteaux D, Collings N. Intake port phenomena in a spark ignition engine at part load. SAE Paper 912401. Society of Automotive Engineers, 1991. [2] Blumberg PN, Kummer JT. Predictions of NO formation in spark-ignited engines—an analysis of methods of control. Combust Sci Tech 1971;4:73–95. [3] Heywood JB, Higgins JM, Watts PA, Tabaczynski RJ. Development and use of a cycle simulation to predict SI engine efficiency and NOx emissions. SAE Paper 790291. Society of Automotive Engineers, 1979. [4] Heywood JB. Internal combustion engine fundamentals. New York: McGraw-Hill; 1988. [5] Benson RS, Garg RD, Woollatt DA. Numerical solution of unsteady flow problems. Int J Mech Sci 1964;6(1):117–44. [6] Chapman M, Novak JM, Stein RA. Numerical modeling of inlet and exhaust flows in multi-cylinder internal combustion engines. In: Uzkan Teoman, editor. Flows in internal combustion engines. ASME winter annual meeting. Phoenix, Arizona: American Society of Mechanical Engineers; 1982. [7] Onorati A, Ferrari G. Modeling of 1-D unsteady flows in I.C. engine pipe systems: numerical methods and transport of chemical species. SAE Paper 980782. Society of Automotive Engineers, 1998. [8] Onorati A, Ferrari G, D’errico G. Fluid dynamic modeling of the gas flow with chemical specie transport through the exhaust manifold of a four cylinder SI engine. SAE Paper 1999-01-0557. SAE SP-1451 SI engine modeling. Society of Automotive Engineers, 1999. [9] Baruah PC. Numerical solution of non-steady flows with variable specific heats and chemical reactions (including use of emission control devices). In: Horlock JH, Winterbone DE, editors. The thermodynamics and gas dynamics of internal-combustion engines, vol. II. New York: Oxford University Press; 1986, p. 869–902. [10] Pearson RJ, Winterbone DE. Calculating the effects of variations in composition on wave propagation in gases. Int J Mech Sci 1993;35(6):517–37. [11] Zhang H, Widener SK. An integrated engine cycle simulation model with species tracking in piping system. SAE Paper 960077. Society of Automotive Engineers, 1996. [12] Arpaci VS, Selamet A, Kao S. Introduction to heat transfer, 1st ed. 2nd print. New Jersey: Prentice-Hall; 2000.
E.E. Selamet et al. / Energy 29 (2004) 449–465
465
[13] Olikara C, Borman GL. A computer program for calculating properties of equilibrium combustion products with some applications to I.C. engines. SAE Paper 750468. Society of Automotive Engineers, 1975. [14] McBride BJ, Gordon S. Computer program for calculation of complex chemical equilibrium compositions and applications. II. Users manual and program description, NASA reference publication 1311; Also, PAC91 by McBride BJ: Properties and coefficients 1991, Cosmic Program # LEW-15779, 1992. [15] Selamet A, Dickey NS, Novak JM. The Herschel-Quincke tube: a theoretical, computational, and experimental investigation. J Acoust Soc Am 1994;96(5):3177–85. [16] Ibrahim AM, Abdel-Rahim YM, Ibrahim IG. Thermochemical properties of combustion gases of ethanol-octane and methanol-octane blends. SAE Paper 911311. Society of Automotive Engineers, 1991. [17] Miller R, Davis G, Lavoie G, Newman C, Gardner T. A super-extended Zel’dovich mechanism for NOx modeling and engine calibration. SAE Paper 980781. Society of Automotive Engineers, 1998. [18] Miller JA, Bowman CT. Mechanism and modeling of nitrogen chemistry in combustion. Prog Energy Combust Sci 1989;15:287–338. [19] Selamet A, Kothamasu V, Novak JM, Kach RA. Experimental investigation of in-duct insertion loss of catalysts in internal combustion engines. Appl Acoust 2000;60:451–87. [20] Zacharias F. Analytical representation of the thermodynamic properties of combustion gases. SAE Paper 670930. Society of Automotive Engineers, 1967.