- Email: [email protected]
Numerical simulation of spark ignition including ionization
Proceedings of the Combustion Institute, Volume 28, 2000/pp. 1177–1185
NUMERICAL SIMULATION OF SPARK IGNITION INCLUDING IONIZATION ¨ RGEN WARNATZ1 and ULRICH MAAS2 MAREN THIELE,1 STEFAN SELLE,1 UWE RIEDEL,1 JU 1
Interdisziplina¨res Zentrum fu¨r Wissenschaftliches Rechnen Universita¨t Heidelberg Im Neuenheimer Feld 368, D-69120 Heidelberg, Germany 2 Institut fu¨r Technische Verbrennung Universita¨t Stuttgart Paffenwaldring 12, D-70569 Stuttgart, Germany
A detailed understanding of the processes associated with spark ignition, as a first step during combustion, is of great importance for clean operation of spark ignition engines. In the past 10 years, a growing concern for environmental protection, including low emission of pollutants, has increased the interest in the numerical simulation of ignition phenomena to guarantee successful flame kernel development even for lean mixtures. However, the process of spark ignition in a combustible mixture is not yet fully understood. The use of detailed reaction mechanisms, combined with electrodynamical modeling of the spark, is necessary to optimize spark ignition for lean mixtures. This work presents the simulation of the coupling of flow, chemical reactions, and transport with discharge processes including ionization in order to investigate the development of a stable flame kernel initiated by an electrical spark in methane/air mixtures. A transport model taking into account the interactions of charged particles has been incorporated in the flow model. This model is based on the Chapman-Enskog theory with an extension for polyatomic gases and considers resonant charge transfer and ambipolar diffusion for the computation of the transport coefficients. A two-dimensional code to simulate the early stages of flame development, shortly after the breakdown discharge, has been developed. The modeling includes an equation for the electrical field. The spark plasma channel left behind by the breakdown is incorporated into the initial conditions. Due to the fast expansion of the plasma channel, a complicated flowfield develops after the emission of a shock wave by the expanding channel. The second phase, that is, the development of a propagating flame and the flame kernel expansion, can last up to several milliseconds and is dominated by diffusive processes and chemical reactions.
Introduction One of the major tasks in combustion technology is the reduction of pollutant emission, such as NOx, CO, or soot [1]. Designing modern internal combustion engines includes the optimization of the cylinder and piston geometry as well as exploiting new technologies such as direct injection and lean combustion. By approaching the lean ignition limit, the formation of pollutants can be reduced, but some cycles do not burn completely or do not ignite at all. As the first step in a combustion process, ignition has great impact on the development of the flame kernel. The discharge in the small gap between the electrodes depends on many parameters, such as spark energy, gas composition, heat losses, the flowfield, and many more. In the first phase of the ignition process, the properties of the combustible mixture are dominated by the shock wave, which is emitted during the expansion of the hot plasma channel. It induces a characteristic flowfield and a growing hot channel, from
which ignition is initiated. In a complete discharge process, only a part of the deposited energy is available for the initiation of the ignition, whereas the other part is lost due to heat conduction processes, the expanding blast wave, and radiation [2]. Experiments (e.g., Ref. [1]), show that rapid deposition of energy during the first phase reduces the induction period. The influence of different temperature profiles was simulated by Ref. [3], but the work excluded detailed chemical kinetics. Several authors have investigated the influence of spark energy and power on the structure of the kernel [3–7]. The geometry of the electrode (width, shape, and gap size) has a strong impact on the development of the flame kernel itself, as shown experimentally in Ref. [8] or numerically in Refs. [9,10], but also on its shape [11,12]. The electrodes account also for heat losses through conduction [13,14], which might even cause quenching. Not only the electrical and geometrical parameters influence the development of the flame kernel, but also the flowfield, governed by the shock
1177
1178
COMBUSTION IN ENGINES
Fig. 1. Computational region, electrode geometry, location of the plasma channel, and boundary treatment.
wave [15–17]. Shortly after breakdown, chemical reactions must produce enough energy to overcome heat losses; the kernel has to grow beyond a critical size or quenching will occur [8,18]. Experimental studies divide the spark ignition process into three phases: breakdown, arc, and glowdischarge, which all have particular electrical properties [2]. Theoretical studies sometimes separate the overall process according to numerical aspects in order to simplify the mathematical treatment. The heating of the plasma channel is sometimes considered as a constant volume process followed by the shock-determined expansion phase [19,20]. Then, as soon as the rapid movement has slowed down, the diffusive phase starts. In this work, we perform numerical simulations of the spark ignition process, including the calculation of the electric field, detailed chemical reactions, and transport models of neutral as well as ionized species, taking into account the coupling of the gasphase processes with the discharge via the electric conductivity. The interaction of these processes is very close. All of them have to be considered in order to obtain a broader understanding of the spark ignition process. Model Description The early development of the flame kernel initiated by an electrical spark is simulated numerically by solving the set of conservation equations for reacting flows, including detailed chemistry [21]. Because the early shape of the flame kernel is assumed to have a rotational symmetry, we use a cylindrical computational domain. This reduces the geometry
to a two-dimensional coordinate system, which is shown in Fig. 1. The electrodes are modeled as cylinders with 45 tips. Setting up an appropriate set of conservation equations for the reactive flow (see below), we assume that the electrodes have a uniform temperature distribution, that the influence of the magnetic field is negligible, and that the plasma is transparent to its own radiation. To describe the reactive flow, the complete compressible Navier-Stokes equations together with the energy and species mass conservation were employed [22]. In the energy conservation equation, the coupling with the electrodynamic equations is modeled by a source term q˙e, which accounts for energy deposition by the discharge due to resistive heating and is given by q˙e ⳱ r • |grad U|, where U is the potential and r the conductivity, and assuming that the plasma is electrically neutral [23]. The gradient of the potential can be obtained by a solution of Maxwell’s equations [23]. Assuming that the medium is isotropically conducting, and that the field attains a quasi-steady state, these equations can be reduced to the simple, time-independent equation div(ⳮr grad U) ⳱ 0
(1)
After the breakdown, the plasma channel is modeled as a cylinder of hot gases [2,17]. Therefore, during this phase one can assume that electrical current flows in the z direction only, that is, the radial component of the electric field inside the kernel is much smaller than the axial component. Outside the channel, the heating by the electric current is negligible due to vanishing conductivity. In this case, the Joule heat source can be calculated directly from the total electrical current I [16,17] by q˙e ⳱ r
I2
冢冮
2
Rc
0
冣
(2)
2prr(r, z)dr
As soon as the heated gas leaves the electrode gap, the radial component of the electric field becomes more and more important, because the flame kernel develops a spherical shape. The assumption of a negligible radial field component is therefore no longer justified, and equation 1 has to be solved to recover the electrical field. Initial and Boundary Conditions The calculated properties of the plasma channel shortly after breakdown [17,19] were incorporated as an initial profile at the beginning of the simulation. It was assumed that local thermal equilibrium exists, which is true a few nanoseconds after the breakdown [2]. The first set of computations in this work uses the temperature and pressure profiles calculated in Ref. [19] as input. There, the channel has cooled down to about 6000 K, with a pressure of
SPARK IGNITION INCLUDING IONIZATION
1179
Fig. 2. Validation of the results of our simulations with experiments [2] and simulations [17, 19]. Left, CH4/air, stoichiometric (this work, [2], [18]); right, CH4/air (this work), pure air ([12,18]). In both figures, the upper set of lines shows pressure development. The lower set of lines shows flame kernel development.
around 10 bar. The dimensions of the computational region are 2.5 ⳯ 2.5 mm. The electrodes’ diameter is 1 mm, and the gap size is 1 mm. From now on, this geometry is referred to as the first configuration. The second set of computations uses the channel properties measured in Ref. [2]. The maximal temperature at 60 ns after breakdown is about 35,000 K, with a corresponding pressure maximum of about 100 bar. The computational regions’ dimensions are 4 ⳯ 4 mm. The electrodes’ diameter is 2 mm, and the gap size is 2 mm. This geometry is referred to as the second configuration throughout the paper. Initial profiles for the species were calculated from the unburned gas according to the channel’s properties. The cold mixture surrounding the plasma channel has an initial temperature of 300 K and is at atmospheric pressure. The initial mixture is quiescent. The electrodes are incorporated as walls in the computational region (Fig. 1). The boundary r ⳱ 0 is a line of symmetry for all variables, including the scalar potential U. The same is true for z ⳱ Z, except for U. Because a reference potential has to be defined for U, the equipotential surface U ⳱ 0 is located at z ⳱ Z. The other two boundaries are treated as outer boundaries, where no-slip conditions are implemented for the velocities and vanishing normal derivatives for the other variables. For the symmetry lines, vanishing normal derivatives are implemented for all variables except the velocity component normal to the symmetry line, which is set to zero. Special boundary conditions are formulated at the electrode boundaries to account for heat transfer and for current input through the electrode at z ⳱ 0.
To account for heat losses, three different types of boundary conditions were implemented: no heat loss to electrodes (Neumann condition), T/n ⳱ 0; constant wall temperature, T ⳱ Tb ⳱ 300 K (Dirichlet); ¯ ⳮ Tb) ⳮ k grad T ⳱ 0. For or mixed conditions h(T the electric current input, a time-dependent profile was assumed for the total current I, which is either modeled as linearly decreasing [19] (second configuration) or as step function (first configuration):
冦
冧
6 A, t ⱕ 0.3 ls I1st(t) ⳱ 0.2 A, 0.3 ls ⱕ t ⱕ 100 ls 0, t 100 ls
冢
I2nd(t) ⳱ 80 mA max 1 ⳮ
冣
t , 0 tign
This leads to the boundary condition U/n ⳱ ⳮ I(t)/(rpr2ign) for the scalar potential. For the one-dimensional electric source model, U is not needed, since the source term can be computed directly from equation 2. Reaction Mechanism and Transport Properties A detailed reaction mechanism consisting of 35 species and 288 elementary reactions for methane/ air mixture was used [24]. The transport model is based on the theory of dilute gases [21]; details on the evaluation of the transport coefficients can be found in Ref. [25]. A comparison between three multicomponent formulations was carried out. The frequently used (1,ⳮ1) mixing rules [26] were tested against the method proposed in Ref. [27], which uses the first
1180
COMBUSTION IN ENGINES
Fig. 3. Traveling pressure wave for the first configuration (top) at 0.2 ls and 2 ls and for the second configuration (bottom) at 0.5 ls and 5 ls for a CH4/air mixture (U ⳱ 0.8). Ignition time, 100 ls. Vectors indicate velocities. Stream traces indicate flow and vortice development.
SPARK IGNITION INCLUDING IONIZATION
1181
Fig. 4. Development of the flame kernel at 0.3 ls, 2 ls, 10 ls, and 30 ls for a CH4/air mixture (U ⳱ 0.8). First configuration.
iterate in a conjugate gradient solution of the transport system, and against solving the complete transport linear system [26]. Since the solution of the complete linear system is very time consuming compared to the other two and the improvement is negligible if the first iterate is used instead of the (1,ⳮ1) rules, all computations were done using the approximations in Ref. [27]. The electrodynamical model must be completed by an equation to derive the electrical conductivity
r. Since apart from transport of ionized air very little data can be found in the literature, the transport models are based on the properties of ionized air [28]. As a first approximation for the electrical conductivity, a two-dimensional piecewise polynomial fit in P and T of degree two based on this data was chosen. This fit is limited to temperatures lower than 7000 K. For higher temperatures, shape-preserving cubic piecewise Hermite interpolation in T together with a spline interpolation in P was used.
1182
COMBUSTION IN ENGINES
Fig. 5. Profiles of some major species in a CH4/air mixture (U ⳱ 0.8) at t ⳱ 1 ls and t ⳱ 100 ls. Ignition time, 100 ls.
Numerical Solution The differential-algebraical equation system (DAE) presented above has to be solved numerically. The extremely differing time and length scales lead to very stiff systems. Discretization in space and time was treated separately by the method of lines. Space discretization was done with central finite differences using five-point stencils for all derivatives but mixed second derivatives, as in the dissipative terms, where nine-point stencils were employed. Since the pressure wave has an important influence on the flowfield in the beginning, it has to be treated accurately by the scheme. In order to cope with the steep gradients at the shock fronts, an artificial viscosity method was used, which broadens the wave but does not influence the flowfield. Implicit time integration methods were applied due to stiffness. The extrapolation solver LIMEX [29], designed for stiff DAEs, which includes step size and order control, was used for solving the timedependent system. The discretization scheme leads to block-nonadiagonal structures of the Jacobian and the iteration matrices [30]. Results Validation During the breakdown phase, the gas temperature in a very small channel (40 lm) rises to several thousands of degrees. Because the time is too short for the pressure to equilibrate, it rises immediately to several tens of mega-Pascals. A shock wave is
formed and the plasma channel starts expanding violently, causing the channel to cool down. During the first 5 ls, the flowfield is dominated by the processes induced by the shock wave. To validate our code, we compared the experimental results from Ref. [2] for a stoichiometric methane/air mixture with our computations for the same setting during this initial phase. For the first configuration, Fig. 2 (left) shows the pressure wave and the flame kernel development from Ref. [2] together with the simulations in Ref. [17] in comparison with our results. Fig. 2 (right) shows the same for the second configuration, together with the simulations in Ref. [19]. The pressure wave is well resolved in both cases, whereas in our computations the flame kernels radius is a little too big, but within tolerable agreement, because the actual speed of the pressure wave is consistent with the experiment. Lean Methane/Air Mixture Figure 3 shows a two-dimensional contour plot of the traveling pressure wave for a lean (U ⳱ 0.8) CH4/air mixture, an ambient pressure of 1•105 Pa and an initial temperature of 300 K for both configurations. The development of the shock wave, the beginning rarefaction wave, and the transition of the expanding channel from cylindrical to spherical configuration is clearly visible in both cases. The pressure maxima decrease very fast, as typical for cylindrical expanding shocks. Due to the very high initial temperatures in the plasma channel, the pressure wave is much stronger in the first case. The rotating vortices which develop near the electrodes’ corners as well as the beginning of recirculation caused by
SPARK IGNITION INCLUDING IONIZATION
1183
Fig. 6. Influence of ions on the temperature (M1, ----; M2, •• –) and on the temperature source term (M1, – – –, M2, ——) (upper part). Mole fraction of key species (lower part) along the symmetry line between the electrodes at 10 ls (left) and 25 ls (right). Symbols: O, ⵧ; OH, 䉭; CH3, 䉮; CHO, (no symbol); O2, 䊊; and NOⳭ, -----. Scale: 10 ls: xNOⳭ • 40 (M2); 25 ls: xCHO • 5, xNOⳭ • 5 • 104 (M2).
fresh gas flowing inward along the electrode walls are clearly visible in Fig. 3 (bottom right), for the second configuration. Due to increased friction, rotating vortices do not appear for the first configuration. The shock wave as well as the hot kernel undergoes a transition from cylindrical to spherical geometry (see Fig. 4). Due to the increase in axial velocity and the decrease in radial velocity, when the shock leaves the gap, the expansion slows down. After the initial phase, a flame kernel is formed. Fig. 4 shows four stages of its development at 0.3 ls, 2 ls, 10 ls, and 30 ls, starting out as cylinder, but changing into spherical shape. While the ignition source
still supplies energy, the temperature stays around 5000 to 6000 K inside the hot kernel. After switching off the ignition source at 100 ls, the maximal kernel temperature gradually decreases. The expansion velocity of the kernel also decreases. The speed of the kernel expansion approaches the burning velocity after the spark is turned off. Figure 5 shows the mass fractions of some species for the first configuration at t ⳱ 1 ls and at t ⳱ 100 ls in radial direction near the line of symmetry between the electrodes. At 1 ls, high radical concentrations are found at the transition region between the plasma channel and the surrounding unburned mixture, but also in the plasma channel itself. At 100
1184
COMBUSTION IN ENGINES
ls, high radical concentrations mark the location of the flame front. Species such as O, H, or CO can mainly be found inside the channel, where the temperature is still higher than 4000 K. The steep profiles of O2 and CH4 indicate the location of the flame front. The profiles show how the dissociation mechanism and the combustion reactions intersect. In the spark channel (T 3000 K), species, such as H2O, which are combustion products, are dissociated. The Influence of Ions To investigate the influence of ions in the igniting mixture, the second configuration was reexamined for the lean methane/air mixture and the same initial condition as discussed in the section above. The reaction scheme was extended by seven species (NⳭ 2 , NOⳭ, NⳭ, OⳭ, H3OⳭ, HCOⳭ, COⳭ). No ions of CH4 species were included because air is the main component of the mixture, and little data are available on ionization of CH4. The reaction scheme including the ions is referred to as M2; M1 is the scheme without ions. Figure 6 shows the profiles of temperature and chemical source term (upper part) as well as the profiles of several species (lower part) along the symmetry line between the electrodes after 10 ls and 25 ls. The ignition time in this example is 10 ls. It can be seen that in the early phase of the process the formation of ions (consuming energy) leads to an about 150 K lower maximum temperature. However, this energy used for the ionization is released later when the temperature decreases. After 25 ls (Fig. 6, upper right), the flame kernel of the calculation M2 is about 100 K hotter than in the M1 case. It can also be seen from the lower part of Fig. 6 that with scheme M2, significantly more radical species, such as O, are formed, and the concentration of reactants, such as O2, is less in the center region. To summarize, neglecting ionization in the simulation, a faster cooling of the inner kernel and a smaller activated volume were observed.
Conclusions A model for the coupling of a two-dimensional cylindrical reactive flow with electrodynamics was presented. This allows a study of the early phase of the development of a stable flame kernel, including the emission of the pressure wave by the expanding hot, conducting channel. The detailed description of the electric field reveals insight into the behavior of the discharge process. The computations show that the influence of the shock wave on all processes dominates during the first few microseconds (less than 5 ls) and is then followed by a diffusion-controlled phase during which a flame front develops at the rim of the plasma
channel. The main emphasis of the work, however, has been to present the results of a model which allows to couple the electrical process with chemical reactions, flow, and transport during spark ignition. Further calculations using this model will be able to provide a better understanding of the complex processes associated with spark ignition. Acknowledgments This work was financially supported by Deutsche Forschungsgemeinschaft, the Fonds der Chemischen Industrie, and the Max-Buchner-Forschungsstiftung. REFERENCES 1. Maly, R., and Vogel, M., Proc. Combust. Inst. 17:821– 831 (1978). 2. Maly, R., Fuel Economy in Road Vehicles Powered by Spark Ignition Engines, (J. C. Hilliard and G. S. Springer, eds.), Plenum Press, New York, 1984. 3. Bradley, D., and Lung, F. K.-K., Combust. Flame 47:71–93 (1987). 4. Sher, E., and Rafael, S., Proc. Combust. Inst. 19:251– 257 (1982). 5. Lim, M. T., Anderson, R. W., and Arpaci, V. S., Combust. Flame 69:303–313 (1987). 6. Maly, R., Proc. Combust. Inst. 18:1747–1754 (1981). 7. Pischinger, S., and Heywood, J. B., A Study of Flame Development and Engine Performance with Breakdown Ignition Systems in a Visualization Engine, SAE Paper 88-0518. 8. Xu, J., “Untersuchung der Funkenzu¨ndung mittels 20Laserinduzierter Fluoreszenz von Otf-Radikalen,” Dissertation, Institut fu¨r Technische Verbrennung, Universita¨t Stuttgart, Germany, 1995. 9. Ishii, K., Tsukamoto, T., Ujiie, Y., and Kono, M., Combust. Flame 91:153–164 (1992). 10. Kono, M., Ishii, K., Niu, K., Tsukamoto, T., and Ujiie, Y., Prog. Astronaut. Aeronaut. 131:55–70 (1991). 11. Pitt, P. L., Clements, R. M., and Topham, D. R., Combust. Sci. Technol. 78:289–314 (1991). 12. Kono, M., Niu, K., Tsukamoto, T., and Ujiie, Y., Proc. Combust. Inst. 22:1643–1649 (1988). 13. Pischinger, S., and Heywood, J. B., How Heat Losses to the Spark Plug Electrodes Affect Flame Kernel Development in an SI Engine, SAE report 90-0021. 14. Ko, Y., and Anderson, R. W., Electrode Heat Transfer During Spark Ignition, SAE report 89-2083. 15. Borghese, A., Diana, M., Moccia, V., and Tamai, R., Combust. Sci. Technol. 76:219–231 (1991). 16. Akram, M., AIAA J. 34(9):1835–1842 (1996). 17. Kravchik, T., Sher, E., and Heywood, J. B., Combust. Sci. Technol. 108:1–30 (1995). 18. Akindelle, O. O., Bradley, D., Mak, P. W., and McMahon, M., Combust. Flame 47:129–155 (1982). 19. Scha¨fer, M., “Der Zu¨ndfunke,” Dissertation, Institut fu¨r Physikalische Elektronik, Universita¨t Stuttgart, Germany, 1997.
SPARK IGNITION INCLUDING IONIZATION 20. Sher, E., Ben-Ya’ish, J., and Kravchik, T., Combust. Flame 89:186–194 (1992). 21. Hirschfelder, J. O., and Curtiss, C. F., Proc. Combust. Inst. 3:121–127 (1949). 22. Maas, U., and Warnatz, J., IMPACT Comput. Sci. Eng. 1:394–420 (1989). 23. Jackson, J. D., Classical Electrodynamics, John Wiley & Sons, New York, 1975. 24. Karbach, V., “Validierung eines detaillierten Reaktionsmechanismus zur Oxidation von Kohlenwasserstoffen bei hohen Temperaturen,” Master thesis, Interdisziplina¨res Zentrum fu¨r Wissenschaftliches Rechnen, Universita¨t Heidelberg, Germany, 1997. 25. Selle, S., and Riedel, U., Ann. N. Y. Acad. Sci. 891:72– 80 (1999).
1185
26. Hirschfelder, J. O., Curtiss, C. F., and Bird, R. B., Molecular Theory of Gases and Liquids, John Wiley & Sons, New York, 1964. 27. Ern, A., and Giovangigli, V., J. Comp. Phys. 120:105– 116 (1995). 28. Yos, J. M., Transport Properties of Nitrogen, Hydrogen, Oxygen and Air, tech. memo. RAD TM-63-7, AVCO Corp., Wilmington, MA, 1963. 29. Deuflhard, H. E., and Zugck, J., One-Step and Extrapolation Methods for Differential/Algebraical Systems, SFB 123 report No. 318, Universita¨t Heidelberg, Germany, 1985. 30. Riedel, U., Maas, U., and Warnatz, J., IMPACT Comput. Sci. Eng. 5:20–52 (1993).
COMMENTS Dumitru Oancea, University of Bucharest, Romania. After the energy deposition, there is a continuos evolution of temperature, pressure and composition around the initial channel, and so on. How is the ignition time defined in your model? Author’s Reply. The total current input is modeled as step function or as linearly decreasing. However, the model is able to simulate other temporal profiles of energy deposition, too. The ignition time is the interval between breakdown and the end of the current input, that is, the energy deposition time.
ployed, an elliptical ignition kernel is formed. This ignition kernel overexpands in the radical coordinate normal to the direction of the laser. The overexpanding ignition kernel forms a torus in shape. However, if the laser pulse has a relatively large energy, a hot jet is observed leaving the inside of the torus. This jet evolves into a lobe. If the breakdown occurs in a flammable mixture, this lobe can increase the burning velocity significantly. Author’s Reply. The model is able to simulate laser-induced ignition as has been shown in a previous publication [1], and indeed it would be interesting to simulate this phenomenon.
● Volker Sick, University of Michigan, USA. Is there experimental evidence for the large difference in computed initial temperatures when including ionization? Author’s Reply. The differences in temperature between the two computations (with and without ions) depend on the energy input (ignition time, current input). Unfortunately, to the best of our knowledge, there is no time-resolved experimental data available for the temperature. ● Russel Lockett, City University, UK. It would be interesting to simulate a laser pulse ignition with the model you have developed. If a laser with low pulse energy is em-
REFERENCE 1. Behrendt, F., Goyal, G., Maas, U., and Warnatz, J., Proc. Combust. Inst. 24:83–91 (1992). ● A. Mokhov, University of Groningen, The Netherlands. At relatively low temperatures the negative ions such as OH, Oⳮ and HCOⳮ 3 are very important in kinetics of formation of charged components. Are you going to incorporate these ions in your model? Author’s Reply. Negative ions will be included in the model in future simulations.
NUMERICAL SIMULATION OF SPARK IGNITION INCLUDING IONIZATION ¨ RGEN WARNATZ1 and ULRICH MAAS2 MAREN THIELE,1 STEFAN SELLE,1 UWE RIEDEL,1 JU 1
Interdisziplina¨res Zentrum fu¨r Wissenschaftliches Rechnen Universita¨t Heidelberg Im Neuenheimer Feld 368, D-69120 Heidelberg, Germany 2 Institut fu¨r Technische Verbrennung Universita¨t Stuttgart Paffenwaldring 12, D-70569 Stuttgart, Germany
A detailed understanding of the processes associated with spark ignition, as a first step during combustion, is of great importance for clean operation of spark ignition engines. In the past 10 years, a growing concern for environmental protection, including low emission of pollutants, has increased the interest in the numerical simulation of ignition phenomena to guarantee successful flame kernel development even for lean mixtures. However, the process of spark ignition in a combustible mixture is not yet fully understood. The use of detailed reaction mechanisms, combined with electrodynamical modeling of the spark, is necessary to optimize spark ignition for lean mixtures. This work presents the simulation of the coupling of flow, chemical reactions, and transport with discharge processes including ionization in order to investigate the development of a stable flame kernel initiated by an electrical spark in methane/air mixtures. A transport model taking into account the interactions of charged particles has been incorporated in the flow model. This model is based on the Chapman-Enskog theory with an extension for polyatomic gases and considers resonant charge transfer and ambipolar diffusion for the computation of the transport coefficients. A two-dimensional code to simulate the early stages of flame development, shortly after the breakdown discharge, has been developed. The modeling includes an equation for the electrical field. The spark plasma channel left behind by the breakdown is incorporated into the initial conditions. Due to the fast expansion of the plasma channel, a complicated flowfield develops after the emission of a shock wave by the expanding channel. The second phase, that is, the development of a propagating flame and the flame kernel expansion, can last up to several milliseconds and is dominated by diffusive processes and chemical reactions.
Introduction One of the major tasks in combustion technology is the reduction of pollutant emission, such as NOx, CO, or soot [1]. Designing modern internal combustion engines includes the optimization of the cylinder and piston geometry as well as exploiting new technologies such as direct injection and lean combustion. By approaching the lean ignition limit, the formation of pollutants can be reduced, but some cycles do not burn completely or do not ignite at all. As the first step in a combustion process, ignition has great impact on the development of the flame kernel. The discharge in the small gap between the electrodes depends on many parameters, such as spark energy, gas composition, heat losses, the flowfield, and many more. In the first phase of the ignition process, the properties of the combustible mixture are dominated by the shock wave, which is emitted during the expansion of the hot plasma channel. It induces a characteristic flowfield and a growing hot channel, from
which ignition is initiated. In a complete discharge process, only a part of the deposited energy is available for the initiation of the ignition, whereas the other part is lost due to heat conduction processes, the expanding blast wave, and radiation [2]. Experiments (e.g., Ref. [1]), show that rapid deposition of energy during the first phase reduces the induction period. The influence of different temperature profiles was simulated by Ref. [3], but the work excluded detailed chemical kinetics. Several authors have investigated the influence of spark energy and power on the structure of the kernel [3–7]. The geometry of the electrode (width, shape, and gap size) has a strong impact on the development of the flame kernel itself, as shown experimentally in Ref. [8] or numerically in Refs. [9,10], but also on its shape [11,12]. The electrodes account also for heat losses through conduction [13,14], which might even cause quenching. Not only the electrical and geometrical parameters influence the development of the flame kernel, but also the flowfield, governed by the shock
1177
1178
COMBUSTION IN ENGINES
Fig. 1. Computational region, electrode geometry, location of the plasma channel, and boundary treatment.
wave [15–17]. Shortly after breakdown, chemical reactions must produce enough energy to overcome heat losses; the kernel has to grow beyond a critical size or quenching will occur [8,18]. Experimental studies divide the spark ignition process into three phases: breakdown, arc, and glowdischarge, which all have particular electrical properties [2]. Theoretical studies sometimes separate the overall process according to numerical aspects in order to simplify the mathematical treatment. The heating of the plasma channel is sometimes considered as a constant volume process followed by the shock-determined expansion phase [19,20]. Then, as soon as the rapid movement has slowed down, the diffusive phase starts. In this work, we perform numerical simulations of the spark ignition process, including the calculation of the electric field, detailed chemical reactions, and transport models of neutral as well as ionized species, taking into account the coupling of the gasphase processes with the discharge via the electric conductivity. The interaction of these processes is very close. All of them have to be considered in order to obtain a broader understanding of the spark ignition process. Model Description The early development of the flame kernel initiated by an electrical spark is simulated numerically by solving the set of conservation equations for reacting flows, including detailed chemistry [21]. Because the early shape of the flame kernel is assumed to have a rotational symmetry, we use a cylindrical computational domain. This reduces the geometry
to a two-dimensional coordinate system, which is shown in Fig. 1. The electrodes are modeled as cylinders with 45 tips. Setting up an appropriate set of conservation equations for the reactive flow (see below), we assume that the electrodes have a uniform temperature distribution, that the influence of the magnetic field is negligible, and that the plasma is transparent to its own radiation. To describe the reactive flow, the complete compressible Navier-Stokes equations together with the energy and species mass conservation were employed [22]. In the energy conservation equation, the coupling with the electrodynamic equations is modeled by a source term q˙e, which accounts for energy deposition by the discharge due to resistive heating and is given by q˙e ⳱ r • |grad U|, where U is the potential and r the conductivity, and assuming that the plasma is electrically neutral [23]. The gradient of the potential can be obtained by a solution of Maxwell’s equations [23]. Assuming that the medium is isotropically conducting, and that the field attains a quasi-steady state, these equations can be reduced to the simple, time-independent equation div(ⳮr grad U) ⳱ 0
(1)
After the breakdown, the plasma channel is modeled as a cylinder of hot gases [2,17]. Therefore, during this phase one can assume that electrical current flows in the z direction only, that is, the radial component of the electric field inside the kernel is much smaller than the axial component. Outside the channel, the heating by the electric current is negligible due to vanishing conductivity. In this case, the Joule heat source can be calculated directly from the total electrical current I [16,17] by q˙e ⳱ r
I2
冢冮
2
Rc
0
冣
(2)
2prr(r, z)dr
As soon as the heated gas leaves the electrode gap, the radial component of the electric field becomes more and more important, because the flame kernel develops a spherical shape. The assumption of a negligible radial field component is therefore no longer justified, and equation 1 has to be solved to recover the electrical field. Initial and Boundary Conditions The calculated properties of the plasma channel shortly after breakdown [17,19] were incorporated as an initial profile at the beginning of the simulation. It was assumed that local thermal equilibrium exists, which is true a few nanoseconds after the breakdown [2]. The first set of computations in this work uses the temperature and pressure profiles calculated in Ref. [19] as input. There, the channel has cooled down to about 6000 K, with a pressure of
SPARK IGNITION INCLUDING IONIZATION
1179
Fig. 2. Validation of the results of our simulations with experiments [2] and simulations [17, 19]. Left, CH4/air, stoichiometric (this work, [2], [18]); right, CH4/air (this work), pure air ([12,18]). In both figures, the upper set of lines shows pressure development. The lower set of lines shows flame kernel development.
around 10 bar. The dimensions of the computational region are 2.5 ⳯ 2.5 mm. The electrodes’ diameter is 1 mm, and the gap size is 1 mm. From now on, this geometry is referred to as the first configuration. The second set of computations uses the channel properties measured in Ref. [2]. The maximal temperature at 60 ns after breakdown is about 35,000 K, with a corresponding pressure maximum of about 100 bar. The computational regions’ dimensions are 4 ⳯ 4 mm. The electrodes’ diameter is 2 mm, and the gap size is 2 mm. This geometry is referred to as the second configuration throughout the paper. Initial profiles for the species were calculated from the unburned gas according to the channel’s properties. The cold mixture surrounding the plasma channel has an initial temperature of 300 K and is at atmospheric pressure. The initial mixture is quiescent. The electrodes are incorporated as walls in the computational region (Fig. 1). The boundary r ⳱ 0 is a line of symmetry for all variables, including the scalar potential U. The same is true for z ⳱ Z, except for U. Because a reference potential has to be defined for U, the equipotential surface U ⳱ 0 is located at z ⳱ Z. The other two boundaries are treated as outer boundaries, where no-slip conditions are implemented for the velocities and vanishing normal derivatives for the other variables. For the symmetry lines, vanishing normal derivatives are implemented for all variables except the velocity component normal to the symmetry line, which is set to zero. Special boundary conditions are formulated at the electrode boundaries to account for heat transfer and for current input through the electrode at z ⳱ 0.
To account for heat losses, three different types of boundary conditions were implemented: no heat loss to electrodes (Neumann condition), T/n ⳱ 0; constant wall temperature, T ⳱ Tb ⳱ 300 K (Dirichlet); ¯ ⳮ Tb) ⳮ k grad T ⳱ 0. For or mixed conditions h(T the electric current input, a time-dependent profile was assumed for the total current I, which is either modeled as linearly decreasing [19] (second configuration) or as step function (first configuration):
冦
冧
6 A, t ⱕ 0.3 ls I1st(t) ⳱ 0.2 A, 0.3 ls ⱕ t ⱕ 100 ls 0, t 100 ls
冢
I2nd(t) ⳱ 80 mA max 1 ⳮ
冣
t , 0 tign
This leads to the boundary condition U/n ⳱ ⳮ I(t)/(rpr2ign) for the scalar potential. For the one-dimensional electric source model, U is not needed, since the source term can be computed directly from equation 2. Reaction Mechanism and Transport Properties A detailed reaction mechanism consisting of 35 species and 288 elementary reactions for methane/ air mixture was used [24]. The transport model is based on the theory of dilute gases [21]; details on the evaluation of the transport coefficients can be found in Ref. [25]. A comparison between three multicomponent formulations was carried out. The frequently used (1,ⳮ1) mixing rules [26] were tested against the method proposed in Ref. [27], which uses the first
1180
COMBUSTION IN ENGINES
Fig. 3. Traveling pressure wave for the first configuration (top) at 0.2 ls and 2 ls and for the second configuration (bottom) at 0.5 ls and 5 ls for a CH4/air mixture (U ⳱ 0.8). Ignition time, 100 ls. Vectors indicate velocities. Stream traces indicate flow and vortice development.
SPARK IGNITION INCLUDING IONIZATION
1181
Fig. 4. Development of the flame kernel at 0.3 ls, 2 ls, 10 ls, and 30 ls for a CH4/air mixture (U ⳱ 0.8). First configuration.
iterate in a conjugate gradient solution of the transport system, and against solving the complete transport linear system [26]. Since the solution of the complete linear system is very time consuming compared to the other two and the improvement is negligible if the first iterate is used instead of the (1,ⳮ1) rules, all computations were done using the approximations in Ref. [27]. The electrodynamical model must be completed by an equation to derive the electrical conductivity
r. Since apart from transport of ionized air very little data can be found in the literature, the transport models are based on the properties of ionized air [28]. As a first approximation for the electrical conductivity, a two-dimensional piecewise polynomial fit in P and T of degree two based on this data was chosen. This fit is limited to temperatures lower than 7000 K. For higher temperatures, shape-preserving cubic piecewise Hermite interpolation in T together with a spline interpolation in P was used.
1182
COMBUSTION IN ENGINES
Fig. 5. Profiles of some major species in a CH4/air mixture (U ⳱ 0.8) at t ⳱ 1 ls and t ⳱ 100 ls. Ignition time, 100 ls.
Numerical Solution The differential-algebraical equation system (DAE) presented above has to be solved numerically. The extremely differing time and length scales lead to very stiff systems. Discretization in space and time was treated separately by the method of lines. Space discretization was done with central finite differences using five-point stencils for all derivatives but mixed second derivatives, as in the dissipative terms, where nine-point stencils were employed. Since the pressure wave has an important influence on the flowfield in the beginning, it has to be treated accurately by the scheme. In order to cope with the steep gradients at the shock fronts, an artificial viscosity method was used, which broadens the wave but does not influence the flowfield. Implicit time integration methods were applied due to stiffness. The extrapolation solver LIMEX [29], designed for stiff DAEs, which includes step size and order control, was used for solving the timedependent system. The discretization scheme leads to block-nonadiagonal structures of the Jacobian and the iteration matrices [30]. Results Validation During the breakdown phase, the gas temperature in a very small channel (40 lm) rises to several thousands of degrees. Because the time is too short for the pressure to equilibrate, it rises immediately to several tens of mega-Pascals. A shock wave is
formed and the plasma channel starts expanding violently, causing the channel to cool down. During the first 5 ls, the flowfield is dominated by the processes induced by the shock wave. To validate our code, we compared the experimental results from Ref. [2] for a stoichiometric methane/air mixture with our computations for the same setting during this initial phase. For the first configuration, Fig. 2 (left) shows the pressure wave and the flame kernel development from Ref. [2] together with the simulations in Ref. [17] in comparison with our results. Fig. 2 (right) shows the same for the second configuration, together with the simulations in Ref. [19]. The pressure wave is well resolved in both cases, whereas in our computations the flame kernels radius is a little too big, but within tolerable agreement, because the actual speed of the pressure wave is consistent with the experiment. Lean Methane/Air Mixture Figure 3 shows a two-dimensional contour plot of the traveling pressure wave for a lean (U ⳱ 0.8) CH4/air mixture, an ambient pressure of 1•105 Pa and an initial temperature of 300 K for both configurations. The development of the shock wave, the beginning rarefaction wave, and the transition of the expanding channel from cylindrical to spherical configuration is clearly visible in both cases. The pressure maxima decrease very fast, as typical for cylindrical expanding shocks. Due to the very high initial temperatures in the plasma channel, the pressure wave is much stronger in the first case. The rotating vortices which develop near the electrodes’ corners as well as the beginning of recirculation caused by
SPARK IGNITION INCLUDING IONIZATION
1183
Fig. 6. Influence of ions on the temperature (M1, ----; M2, •• –) and on the temperature source term (M1, – – –, M2, ——) (upper part). Mole fraction of key species (lower part) along the symmetry line between the electrodes at 10 ls (left) and 25 ls (right). Symbols: O, ⵧ; OH, 䉭; CH3, 䉮; CHO, (no symbol); O2, 䊊; and NOⳭ, -----. Scale: 10 ls: xNOⳭ • 40 (M2); 25 ls: xCHO • 5, xNOⳭ • 5 • 104 (M2).
fresh gas flowing inward along the electrode walls are clearly visible in Fig. 3 (bottom right), for the second configuration. Due to increased friction, rotating vortices do not appear for the first configuration. The shock wave as well as the hot kernel undergoes a transition from cylindrical to spherical geometry (see Fig. 4). Due to the increase in axial velocity and the decrease in radial velocity, when the shock leaves the gap, the expansion slows down. After the initial phase, a flame kernel is formed. Fig. 4 shows four stages of its development at 0.3 ls, 2 ls, 10 ls, and 30 ls, starting out as cylinder, but changing into spherical shape. While the ignition source
still supplies energy, the temperature stays around 5000 to 6000 K inside the hot kernel. After switching off the ignition source at 100 ls, the maximal kernel temperature gradually decreases. The expansion velocity of the kernel also decreases. The speed of the kernel expansion approaches the burning velocity after the spark is turned off. Figure 5 shows the mass fractions of some species for the first configuration at t ⳱ 1 ls and at t ⳱ 100 ls in radial direction near the line of symmetry between the electrodes. At 1 ls, high radical concentrations are found at the transition region between the plasma channel and the surrounding unburned mixture, but also in the plasma channel itself. At 100
1184
COMBUSTION IN ENGINES
ls, high radical concentrations mark the location of the flame front. Species such as O, H, or CO can mainly be found inside the channel, where the temperature is still higher than 4000 K. The steep profiles of O2 and CH4 indicate the location of the flame front. The profiles show how the dissociation mechanism and the combustion reactions intersect. In the spark channel (T 3000 K), species, such as H2O, which are combustion products, are dissociated. The Influence of Ions To investigate the influence of ions in the igniting mixture, the second configuration was reexamined for the lean methane/air mixture and the same initial condition as discussed in the section above. The reaction scheme was extended by seven species (NⳭ 2 , NOⳭ, NⳭ, OⳭ, H3OⳭ, HCOⳭ, COⳭ). No ions of CH4 species were included because air is the main component of the mixture, and little data are available on ionization of CH4. The reaction scheme including the ions is referred to as M2; M1 is the scheme without ions. Figure 6 shows the profiles of temperature and chemical source term (upper part) as well as the profiles of several species (lower part) along the symmetry line between the electrodes after 10 ls and 25 ls. The ignition time in this example is 10 ls. It can be seen that in the early phase of the process the formation of ions (consuming energy) leads to an about 150 K lower maximum temperature. However, this energy used for the ionization is released later when the temperature decreases. After 25 ls (Fig. 6, upper right), the flame kernel of the calculation M2 is about 100 K hotter than in the M1 case. It can also be seen from the lower part of Fig. 6 that with scheme M2, significantly more radical species, such as O, are formed, and the concentration of reactants, such as O2, is less in the center region. To summarize, neglecting ionization in the simulation, a faster cooling of the inner kernel and a smaller activated volume were observed.
Conclusions A model for the coupling of a two-dimensional cylindrical reactive flow with electrodynamics was presented. This allows a study of the early phase of the development of a stable flame kernel, including the emission of the pressure wave by the expanding hot, conducting channel. The detailed description of the electric field reveals insight into the behavior of the discharge process. The computations show that the influence of the shock wave on all processes dominates during the first few microseconds (less than 5 ls) and is then followed by a diffusion-controlled phase during which a flame front develops at the rim of the plasma
channel. The main emphasis of the work, however, has been to present the results of a model which allows to couple the electrical process with chemical reactions, flow, and transport during spark ignition. Further calculations using this model will be able to provide a better understanding of the complex processes associated with spark ignition. Acknowledgments This work was financially supported by Deutsche Forschungsgemeinschaft, the Fonds der Chemischen Industrie, and the Max-Buchner-Forschungsstiftung. REFERENCES 1. Maly, R., and Vogel, M., Proc. Combust. Inst. 17:821– 831 (1978). 2. Maly, R., Fuel Economy in Road Vehicles Powered by Spark Ignition Engines, (J. C. Hilliard and G. S. Springer, eds.), Plenum Press, New York, 1984. 3. Bradley, D., and Lung, F. K.-K., Combust. Flame 47:71–93 (1987). 4. Sher, E., and Rafael, S., Proc. Combust. Inst. 19:251– 257 (1982). 5. Lim, M. T., Anderson, R. W., and Arpaci, V. S., Combust. Flame 69:303–313 (1987). 6. Maly, R., Proc. Combust. Inst. 18:1747–1754 (1981). 7. Pischinger, S., and Heywood, J. B., A Study of Flame Development and Engine Performance with Breakdown Ignition Systems in a Visualization Engine, SAE Paper 88-0518. 8. Xu, J., “Untersuchung der Funkenzu¨ndung mittels 20Laserinduzierter Fluoreszenz von Otf-Radikalen,” Dissertation, Institut fu¨r Technische Verbrennung, Universita¨t Stuttgart, Germany, 1995. 9. Ishii, K., Tsukamoto, T., Ujiie, Y., and Kono, M., Combust. Flame 91:153–164 (1992). 10. Kono, M., Ishii, K., Niu, K., Tsukamoto, T., and Ujiie, Y., Prog. Astronaut. Aeronaut. 131:55–70 (1991). 11. Pitt, P. L., Clements, R. M., and Topham, D. R., Combust. Sci. Technol. 78:289–314 (1991). 12. Kono, M., Niu, K., Tsukamoto, T., and Ujiie, Y., Proc. Combust. Inst. 22:1643–1649 (1988). 13. Pischinger, S., and Heywood, J. B., How Heat Losses to the Spark Plug Electrodes Affect Flame Kernel Development in an SI Engine, SAE report 90-0021. 14. Ko, Y., and Anderson, R. W., Electrode Heat Transfer During Spark Ignition, SAE report 89-2083. 15. Borghese, A., Diana, M., Moccia, V., and Tamai, R., Combust. Sci. Technol. 76:219–231 (1991). 16. Akram, M., AIAA J. 34(9):1835–1842 (1996). 17. Kravchik, T., Sher, E., and Heywood, J. B., Combust. Sci. Technol. 108:1–30 (1995). 18. Akindelle, O. O., Bradley, D., Mak, P. W., and McMahon, M., Combust. Flame 47:129–155 (1982). 19. Scha¨fer, M., “Der Zu¨ndfunke,” Dissertation, Institut fu¨r Physikalische Elektronik, Universita¨t Stuttgart, Germany, 1997.
SPARK IGNITION INCLUDING IONIZATION 20. Sher, E., Ben-Ya’ish, J., and Kravchik, T., Combust. Flame 89:186–194 (1992). 21. Hirschfelder, J. O., and Curtiss, C. F., Proc. Combust. Inst. 3:121–127 (1949). 22. Maas, U., and Warnatz, J., IMPACT Comput. Sci. Eng. 1:394–420 (1989). 23. Jackson, J. D., Classical Electrodynamics, John Wiley & Sons, New York, 1975. 24. Karbach, V., “Validierung eines detaillierten Reaktionsmechanismus zur Oxidation von Kohlenwasserstoffen bei hohen Temperaturen,” Master thesis, Interdisziplina¨res Zentrum fu¨r Wissenschaftliches Rechnen, Universita¨t Heidelberg, Germany, 1997. 25. Selle, S., and Riedel, U., Ann. N. Y. Acad. Sci. 891:72– 80 (1999).
1185
26. Hirschfelder, J. O., Curtiss, C. F., and Bird, R. B., Molecular Theory of Gases and Liquids, John Wiley & Sons, New York, 1964. 27. Ern, A., and Giovangigli, V., J. Comp. Phys. 120:105– 116 (1995). 28. Yos, J. M., Transport Properties of Nitrogen, Hydrogen, Oxygen and Air, tech. memo. RAD TM-63-7, AVCO Corp., Wilmington, MA, 1963. 29. Deuflhard, H. E., and Zugck, J., One-Step and Extrapolation Methods for Differential/Algebraical Systems, SFB 123 report No. 318, Universita¨t Heidelberg, Germany, 1985. 30. Riedel, U., Maas, U., and Warnatz, J., IMPACT Comput. Sci. Eng. 5:20–52 (1993).
COMMENTS Dumitru Oancea, University of Bucharest, Romania. After the energy deposition, there is a continuos evolution of temperature, pressure and composition around the initial channel, and so on. How is the ignition time defined in your model? Author’s Reply. The total current input is modeled as step function or as linearly decreasing. However, the model is able to simulate other temporal profiles of energy deposition, too. The ignition time is the interval between breakdown and the end of the current input, that is, the energy deposition time.
ployed, an elliptical ignition kernel is formed. This ignition kernel overexpands in the radical coordinate normal to the direction of the laser. The overexpanding ignition kernel forms a torus in shape. However, if the laser pulse has a relatively large energy, a hot jet is observed leaving the inside of the torus. This jet evolves into a lobe. If the breakdown occurs in a flammable mixture, this lobe can increase the burning velocity significantly. Author’s Reply. The model is able to simulate laser-induced ignition as has been shown in a previous publication [1], and indeed it would be interesting to simulate this phenomenon.
● Volker Sick, University of Michigan, USA. Is there experimental evidence for the large difference in computed initial temperatures when including ionization? Author’s Reply. The differences in temperature between the two computations (with and without ions) depend on the energy input (ignition time, current input). Unfortunately, to the best of our knowledge, there is no time-resolved experimental data available for the temperature. ● Russel Lockett, City University, UK. It would be interesting to simulate a laser pulse ignition with the model you have developed. If a laser with low pulse energy is em-
REFERENCE 1. Behrendt, F., Goyal, G., Maas, U., and Warnatz, J., Proc. Combust. Inst. 24:83–91 (1992). ● A. Mokhov, University of Groningen, The Netherlands. At relatively low temperatures the negative ions such as OH, Oⳮ and HCOⳮ 3 are very important in kinetics of formation of charged components. Are you going to incorporate these ions in your model? Author’s Reply. Negative ions will be included in the model in future simulations.