A new ignition model for spark-ignited engine simulations

Parallel Computing 27 (2001) 179±200 www.elsevier.com/locate/parco A new ignition model for spark-ignited engine simulations O. Yasßar a,b a Department of Computational Science, State University of New York, 350 New Campus Drive, Room 294 FOB, Brockport, NY 14420, USA 1 b Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA Received 6 December 1999; received in revised form 25 May 2000 Abstract The amount of spark energy deposited into the combustion chamber is key to an optimum ignition as one can end up with mis®res when this energy is low or with other undesired e€ects on engine performance and byproducts when it is high. Experimentally, up to now, no one has been able to correlate the combustion outcome accurately to the spark parameters in a controllable way. Theoretical investigation and computer modeling is leading to a better understanding of how spark ¯ames propagate. A new computational approach to ignition dynamics is presented here for spark-ignited (SI) engine combustion simulations. Our computational model, using the MPI communication library, attempts to solve temporal and spatial equations of the electromagnetic (EM) equations in conjunction with the well-known NavierStokes equations of the standard KIVA-3 engine code. The interaction between the gas and the ¯ame (plasma) kernel in the spark region is computed through the momentum and energy exchange between these two ®elds. Preliminary results show a distinct spatial distribution of physical quantities at the ¯ame front and within the in¯ammation zone. A slight change in the spark discharge current has signi®cant impact on the combustion and emissions. Enhanced accuracy of spark ignition modeling might help us better compute the early ¯ame propagation and its in¯uence on the cyclic variability of engines, potentially leading to design of new spark plugs. Ó 2001 Elsevier Science B.V. All rights reserved. Keywords: Engine combustion simulations; Spark ignition modeling; Parallel computing; Computational ¯uid dynamics 1 Correspondence address. E-mail address: oyasar@brockport.edu (O. Yasßar). 0167-8191/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 8 1 9 1 ( 0 0 ) 0 0 0 9 4 - 6 180 O. Yasßar / Parallel Computing 27 (2001) 179±200 1. Introduction Aside from economical concerns attached to automotive engineering, the environmental concerns and government regulations also necessitate design of better engines. There is ever-increasing pressure on vehicle manufacturers to improve fuel eciency and reduce emission of harmful by-products. The ®eld of spark ignition thus far has su€ered greatly since the random nature of spark discharge has made it dicult to correlate spark energy input and the combustion process output. The breakdown of sparking events depend on the availability of the ®rst free electron in the background gas and this availability is a process that depends on the random motion of particles in the gas [1] unless a laser or other forms of radiation is used to control this mechanism [2,3]. The other option is to raise the applied voltage to a high enough value to force the breakdown. However, higher voltages mean higher corrosion of the spark plug. As a result, the spark gap increases over time, requiring even higher voltages to have a breakdown. This process continues until the electrical limits of the spark system is reached. The ability to reliably ignite the fuel mixture at low voltages with a smaller gap would extend the time before one reaches the highest gap distance and voltages allowed by the system prior to a system breakdown (end of plug lifetime). More work has to be done on spark dynamics to obtain optimum conditions. Experimentally, the combustion output has been observed to be signi®cantly dependent on the changes in plug wire electrical characteristics, spark gap distance, gas pressure, and other variables but up to now no one has been able to correlate the outcome speci®cally and accurately to one of these variables in a controllable way [4]. Despite many studies to analytically relate breakdown voltage and deposited energy to the spark gap parameters and the gas conditions (pressure, temperature), diculties remain in assessing what phase of the spark (namely breakdown, arc, and glow) is more important to the energy delivery [5]. Computer modeling is paving the way to a better understanding of how spark ¯ames propagate, yet most models use rather simple and empirical methods for spark energy deposition. While the determination of the total energy delivery is important, it is also necessary to examine in detail the time evolution of the energy delivery to the gas. This is because the spark gap does not remain at a constant current and impedance during the discharge process, but rather experiences three (or more) distinct phases, breakdown, arc, and glow [5] for which energy is delivered in di€erent amounts at di€erent time scales. Thus, the assumption of a single energy input is insucient. A direct coupling of spark energy dynamics and the gas dynamics is needed. Spark energy deposition occurs at a timescale much shorter than the gas timescale. Current simulation models march with the gas (¯uid) timescale, thus ignoring the spark evolution. They rather use an estimate or a semi-empirical quantity for the amount of spark energy deposited into the ¯uid. It has been observed that the combustion models are strongly sensitive to the values used for spark energy deposition as input data [6]. Several research groups have attempted to incorporate time scales [7,8] into the ignition process. Studies most relevant to our discussions here are Kravchick et al. O. Yasßar / Parallel Computing 27 (2001) 179±200 181 [9], Maly [5], Lu et al., [6], Mantel [10], and zur Loye and Bracco [11]. Lu and coworkers examine e€ect of engine operating conditions on combustion using the standard KIVA ignition model [28]. They vary the XIGNIT parameter (an estimate for spark energy deposition rate until the temperature reaches combustion conditions ± 1600 K), to observe e€ects on the combustion and emissions. The inadequacy of the standard KIVA ignition model is clearly stated by Lu et al. Mantel [10] examines the e€ect of spark plug geometry on the ¯ame propagation and lists important observations. A signi®cant progress to date in incorporating greater detail of the spark dynamics into a combustion model has been recently made by Kravchik et al. [9]. Although formation of a conducting channel was incorporated into their model, several assumptions were made; including a weak self-induced magnetic ®eld, a low current density, and a spatially-uniform electric ®eld. These assumptions warrant further examination. We have observed, in this paper, that the in¯uence of magnetic ®eld and the ¯uid velocity on the current density is signi®cant. Temporal and spatial variations of electromagnetic (EM) ®elds and properties are certainly a signi®cant source of in¯uence on the spark energy deposition, peak plasma temperature and pollutants. More complete studies are needed to ®nd conditions with minimum spark power to operate with lean mixtures so that both the emission and spark lifetime concerns are answered. Lack of a complete breakdown/ignition/ combustion model, according to Kravchik et al., restricts the fundamental understanding of the process and the ability to improve the means toward ecient energy conversion. Validation of spark ignition models is a dicult and perhaps an impossible task due to the inherent uncertainty in the sparking process and thus the lack of con®dence in the accuracy of the experimental data against which a set of computations can be compared. As we stated before, spark input energy in conventional ignition systems occurs over a time scale from tens of nanoseconds (breakdown phase), to microseconds (arc phase), to milliseconds (glow phase), which means that the spark dynamics begin at time scales much shorter than the ¯uid time scale and extends to times comparable to ¯uid ¯ow. The same equipment and chamber conditions would not produce the same spark energy every time [4,12±16]. This uncertainty is inherent (unless controlled) and primarily depends on how and when the ®rst electron starts the breakdown process. One can comfortably assess that the understanding of the processes involved in ignition and self¯ame propagation in internal combustion engines still has a long way to go. On the one hand, e€orts involving engine experiments with measurement instrumentation have concentrated on quick success and the complexity of the mutually dependent in¯uences of the large number of di€erent operational parameters of the real engine prevented any signi®cant breakthrough in the past, thus leading to crude estimates for spark energy and ignition processes. On the other hand, experimental and theoretical work outside the engine has been su€ering from inadequate measurement instrumentation and crude simulation of the engine conditions. Thus progress in the understanding of spark ignition has been slow and much of the e€ort has gone into more easily understood generation of the required high-voltage pulses [5]. 182 O. Yasßar / Parallel Computing 27 (2001) 179±200 2. Computational model A realistic approach to the spark dynamics would be to model its evolution through the governing EM equations and take into account the interaction of the electrical and magnetic ®elds with the ¯uid ®eld. With a sizable magnetic ®eld, the impact of the magnetic pressure on the spark kernel front and its propagation into the chamber might be important. Our experience [19,26] with plasma channels have shown a signi®cant e€ect on the plasma hydrodynamics by heat transfer and magnetic pressure. The heat transfer from the spark to the gas then will be more accurately computed by the amount of work done on the gas by the electrical and magnetic ®elds. One would then march the simulation on the ignition time scale that is much shorter than the ¯uid time scale. Marching with smaller time steps will certainly cost more, but this will ensure a much more accurate energy deposition into the gas. The amount of energy deposited into the gas has a great e€ect on the combustion that follows. The dynamics of self-sustaining and propagating ¯ame front (a thin reaction sheet where exothermic combustion chemical reactions occur) has a profound e€ect on how uniformly the combustion takes place in the chamber. The amount of unburned hydrocarbons then will be a direct result of ¯ame dynamics that is tightly related to an accurate account of the spark formation and its interaction (such as heat transfer and EM forces) with the gas. Di€erent assumptions of the form of the electric ®eld and the form of the mobility of ions has lead to di€erent results [20]. Because of limited computing power in the past, it has not been possible to completely address these issues to permit a comprehensive description of the discharge event during the breakdown, arc and glow phases. We believe that this is now in the realm of possibility with the advent of high performance computing and the development of the our massively parallel KIVA-3 code [21]. The sparking event in KIVA-3 is currently modelled through a single input energy parameter. This parameter, (XIGNIT), is an empirical approximation of spark-energy deposition rate into the chamber, resulting in a model that lacks the details of spark ignition dynamics and its interaction with the ¯uid ®eld and with electrode surface. During the ignition time (speci®ced by the crank angle), the speci®c internal energy of every spatial cell in the ignition region is increased by a factor of XIGNIT  DT for every time interval (DT) until end-of-ignition is reached. If the temperature in the ignition cell(s) reaches a certain temperature (1600 K) before the end of the ignition window, then the energy deposition is terminated. The choice of XIGNIT and ignition temperature is based on crude estimates and cannot account for what really happens. Such estimates and the lack of details mostly come from the randomness in spark breakdown that precludes one from understanding the sparking process both qualitatively and quantitatively. This randomness, coupled with other variables between engine cycles, has tended to confound previous attempts to draw relations between ignition system inputs and engine outputs. Use of experimental techniques which provide a high degree of control and repeatability will signi®cantly aid in the validation of the model development. A sparkcontrol mechanism developed by Sauers et al. [2] has provided a key motivation for O. Yasßar / Parallel Computing 27 (2001) 179±200 183 our e€ort of adding a new spark ignition model into the KIVA-3 engine code. The spark-control mechanism developed by Sauers et al., can produce the same spark energy with an uncertainty of less than 1%, and indeed, to our knowledge, is the only experiment with such a control mechanism. We view this spark ignition experiment as a valuable tool that will enable us to adjust our theory and calibrate our model, as well as permitting us to ®nd ways to incorporate the spark-energy stabalization technology into a real engine. The fundamental equations to model a turbulent, radiative, and magnetized gas phase reactive ¯ow are the continuum time-dependent equations for conservation of the total mass density q, the individual chemical species densities qm , the momentum density qu, and the internal energy density ef . These equations may be written as [3,21,22,28] oq ‡ r
qu† ˆ q_ s ; ot 1†   oqm q ‡ r
qm u† ˆ r
qDr m ‡ q_ c ‡ q_ s dm;s ; ot q 2†     o 1 R qu ‡ 2 q ‡ r
quu ‡ ~ PR ot c ˆ rp A0 r 2=3qj† ‡ r
r ‡ Fs ‡ qg ‡ o f e ‡ eR † ‡ r
ef u ot ˆ pr
u ‡ 1 A0 †r : ru JB ; c r
qf ‡ qR † ‡ A0 q ‡ Q_ c ‡ Q_ s ‡ J
E; 3† 4† where q_ s and q_ c are source terms due to sprays and chemistry , and dm;s ˆ 1 when species m is the same as species of the spray droplet. qR ; ~ PR , and eR are radiation related quantities: heat ¯ux, pressure tensor, and energy density, respectively. J
E is the Joule heating term and is equal to E0
J0 , the rate of Joulean dissipation in ¯uid frame, plus u
J  B=c†, the rate at which the force J  B=c does work on the ¯uid. D is the di€usion coecient, A0 is the ¯ow parameter (0-laminar,1-turbulent),  and j are the turbulent kinetic energy and dissipation rate, and r is the Newtonian viscous stress tensor consisting of the ®rst and second coecientsP of viscosity, l and k. Finally, qf is the ¯uid heat ¯ux vector as qf ˆ KrT qD m hm r qm =q†, with T as the temperature and hm as the speci®c enthalpy of species m. The solution of ¯uid density, momentum and energy ®elds involve interaction terms between the ¯uid ®eld and the others such as spray particles, turbulence, radiative transfer, chemical reactions, and EM ®elds. Consequently, one would need to solve governing equations for those ®elds as well, depending on their in¯uence on the ¯uid motion. When turbulence is considered A ˆ 1†, two turbulence transport equations are solved for j and . The standard j± equations may be written as 184 O. Yasßar / Parallel Computing 27 (2001) 179±200 oqj ‡ r
quj† ˆ ot 2=3qjr
u ‡ r : ru ‡ r
oq ‡ r
qu† ˆ ot    l rj Prj q ‡ W_ s 5† and ‡ 2 c 3 1    l r c3 qr
u ‡ r
Pr  c r : ru j 1 c2 q ‡ cs W_ s ; 6† where W_ s is external source due to interaction with the spray, and c1 ; c2 ; c3 ; Prj , and Pr are constants determined through experiments and theoretical assumptions. The assumption of j± for turbulence modeling in combustion systems is being questioned frequently. The central assumptions in the j± model are the isotropic eddy viscosity concept and the gradient di€usion model. Also a local equilibrium assumption is made, which is valid with limitations at very high Reynolds numbers. Turbulence persists as one of the major theoretical and computational problems of ¯uid dynamics. P P The chemical reactions are symbolized by m amr vm () m bmr vm , where vm represents one mole of species m and amr and bmr P are integral stoichiometric coecients for reaction r. These coecients must satisfy m amr bmr †Wm ˆ 0 to conserve mass in chemical reactions. The database for chemical reactions (kinetic and equilibrium) has to accompany the model and we currently follow the data available in KIVA-3. The interaction of gas with particles is a complicated problem and the interaction terms appear in most governing equations. In essence, one has to compute the particle density ®rst and then derive its moments to obtain the source terms that appear in the ¯uid equations given above. We follow the work in [28] for spray dynamics. The time evolution of f, the particle distribution function, is given by of o o  _ o o _ ‡ ‡ rx
f v† ‡ rv
f F† ‡ fR† ‡ f y† f y † f Td ‡ ot or oTd oy oy_ _ ‡ fbu _ ; ˆ fcoll 7† where the quantities F, R, T_d , and y are the time rates of change of spray droplet's _ and fbu _ are _ The terms fcoll velocity, radius, temperature, and oscillation velocity y. sources due to droplet collisions and breakups. As noted here, f has ten independent variables in addition to time. These variables are the three droplet position components x, three velocity components v, equilibrium radius r (the radius of the droplet if it were a complete sphere), temperature Td , distortion from sphericity y, _ Droplets break up if y is larger than unity and coalesce and the time rate of change y. if the collision impact parameter b is less than a critical value brc . The space-time evolution of the EM ®elds for a medium moving with velocity u is derived from the Maxwell's equations, Ampere's Law, Faraday's Law and Ohm's Law [29,30]; O. Yasßar / Parallel Computing 27 (2001) 179±200 oB gc2 4p ˆ r  u  B† r  J; r  B; r  B ˆ ot c 4p   1 oB 1 rEˆ ; Jˆr E‡ uB ; c ot c 185 8† where r and g are the electrical conductivity and the resistivity (r ˆ 1=g) of the plasma and c is the speed of light. The ®rst equation describes how the magnetic ®eld lines are convected and di€used in the non-relativistic and low-frequency plasma ¯uid. Others basically constitute a relation between the magnetic ®eld B, the current density J, and the electric ®eld (E) measured in the laboratory frame. The e€ects on the ¯uid, as seen in these governing equations, are both EM and mechanical. The magnetic ®eld creates a force on the ¯uid that is equivalent to a magnetic hydrostatic pressure B2 =8p† via the J  B term in the momentum equation. There may be a need for other auxiliary equations or boundary conditions to solve the EM ®eld, depending on the problem under consideration. In case of spark-ignition in an engine or plasma channels in inertial con®nement fusion reactors, either the current pro®le or the breakdown voltage and spark characteristics and channel parameters (conductivity, geometry) need to be known for computing the current density J. The governing equations presented above all have a common form as oQ ˆ ot r
Qu† ‡ r
DrQ ‡ S; 9† where the terms on the RHS are generally convection, di€usion and source terms. These equations need to be discretized both in time and space. There are two paths to discretize the given equations: the di€erential approach and the control volume approach. In the di€erential approach, the equations are considered to be PDEs and it follows a pure mathematical track, and the physical signi®cance of the variables can be lost. In the control volume approach, the conservative nature of the equations is preserved. Therefore, spatial di€erences are done via the control volume here, integrating the terms over the control volume. Most of the governing equations involve gradient and divergence terms as generalized in Eq. (9). Volume integrals are converted into surface area integrals using the divergence theorem and computation of such terms constitute the heart of the di€usion solvers in the code. Block-structured (tensor-product) grid system used in KIVA-3 o€ers the ability of representing the problem domain with multiple blocks as long as the connectivity between blocks is maintained properly. The connectivity is de®ned through indirect addressing and this makes it no longer necessary to store the grid points in a particular order. The connectivity arrays that describe who the neighbor points are in all directions (left, right, back, top, bottom and front). The mesh in each block is made up of arbitrary hexahedrons (cells), the corner of which are vertices. Each block is independent and surrounded by ghost cells in all directions to handle the in¯ow/ out¯ow boundary conditions. Also, cell-face boundary conditions for all six faces of a cell permit greater ¯exibility and simpli®cation in the application of the boundary conditions. The use of ghost cells and cell-face boundary conditions make it possible to apply the same physics, numerics and boundary conditions to the smallest units of O. Yasßar / Parallel Computing 27 (2001) 179±200 186 the domain (cells or blocks) as well as to the whole domain. This generality forms a convenient foundation for a block-wise distribution of our computational model (described in detail at an earlier publication [18,21]) on systems of independent processors. In a domain decomposition among parallel processors, two adjacent processor communicate the values on the shared boundary. The communication primitives (now in MPI rather than the Intel NX as reported earlier by us at [21]) know where to access elements (on the shared boundary) in the memory. Each send is met with a receive operation. Since the access pattern for vertices and cells on the processor boundaries need to be known in advance for communication requirements, one needs to sort storage arrays, separating the elements on the left, right, and so on. The boundary shared between neighbor processors has the same grid points on both sides to assure proper communication. The suggested sorting is done in the pre-processor. The physical grid is surrounded by ghost cells in all six directions. These ghost cells correspond to the real cells of the adjacent blocks residing on other processors, thus creating a bu€er area for storing boundary information that reside on the adjacent processors. The solution of EM equations (in conjunction with other governing equations) goes to the heart of the new computational ignition model reported here. Therefore, we will describe in some detail how Eq. (8) is solved on the same spatial (but slightly di€erent temporal) grid used in our overall combustion scheme. If we integrate the magnetic di€usion (Eq. (8)) over the control volume V,   Z Z oB gc2 dV ˆ r  B dV: r  u  B† 4p V V ot Volume Rintegrals are converted into surface area integrals using the divergence H theorem V ‰r  AŠdV ˆ S dS  A. This results in the following:   Z I oB gc2 dV ˆ rB : dS  u  B† 4p V ot S Although most of our equations are solved in the cartesian (x; y; z) coordinates, we have opted for polar coordinates (r; h; z) for the magnetic di€usion equation. An SI engine piston can be considered as a cylinder with a spark plug placed at the center of the cylinder head. With such a geometry, the spark discharge current ¯ows in the ^z direction (along the axis of the cylinder), creating a magnetic ®eld in the h^ direction. Therefore, we can assume B ˆ Bh ; oBh =oz  0: Since the computational grid is a cylindrical one, there is no diculty in converting between our polar coordinates and the grid related quantities such as the surface area (dS) of the control volume and the ¯uid speed. The terms under the surface integral become rBˆ 1 o rBh †^z and r or u  B ˆ uz Bh ^r ur Bh ^z: Now, multiplying these terms with the area vectors dS ˆ dSr ^r ‡ dSh h^ ‡ dSz ^z will yield O. Yasßar / Parallel Computing 27 (2001) 179±200 oBh gc2 o rBh † dV ˆ dSr ot 4pr or dSr ur Bh dSz uz Bh ; which can be discretized as  gi‡1=2 c2 1 ri‡1 Bn‡1 Bin‡1 Bni i‡1 V ˆ dS i i‡1=2 tn‡1 tn ri‡1 4p ri‡1=2  gi 1=2 c2 1 ri Bn‡1 i dSi 1=2 4p ri 1=2 ri ‡ dSi n‡1 1=2 ui 1=2 Bi 1=2 187 ri Bn‡1 i ri  ri 1 Bn‡1 i 1 ri 1 n‡1 dSi‡1=2 ui‡1=2 Bi‡1=2 ;  10† where the indexes n and i represent the temporal and spatial steps. With the exception of the magnetic di€usion equation, the temporal di€erencing is done in three di€erent phases. Phase A is where most of the source terms are calculated, Phase B is where the di€usion terms are calculated and Phase C is where the advection is computed. All these three phases apply to the same time interval, but they account for di€erent contributions. The solution of the magnetic di€usion equation, on the other hand, is found using an implicit ®nite-di€erencing scheme. The timescale of magnetic di€usion is very short and the use of implicit di€erencing is needed to avoid extremely small timesteps. Finally, the magnetic di€usion equation is organized as n‡1 n di Bn‡1 ‡ ai Bn‡1 i 1 ‡ ci Bi i‡1 ˆ si ; where di ; ci ; ai , and si are cell-related quantities (i.e., conductivity) for control volume i. The solution is found using the Thomas algorithm [3] Bin‡1 ˆ ei‡1 Bn‡1 i‡1 ‡ fi‡1 ; where i ˆ 1; 2; . . . ; M 1 and ei‡1 ˆ ai =di ei ‡ ci and fi‡1 ˆ si di fi †= di ei ‡ ci †. Here, M is the number of spatial cells. The spatial mesh is swept from M 1 down to n‡1 1, therefore requiring BM to be known in advance. Considering that most of the spark discharge current is ¯owing within a region internal to the spatial mesh M, the BM can be simply equal to magnetic ®eld created by a coaxial cable a distance r from n‡1 the center of the cable, that is BM ˆ 2
I=c
rM . Here, I is the total current ¯owing within the spark channel and c is the speed of light. The value of this current depends on spark plug parameters, the applied voltage, and conditions in the chamber. In this paper, we use a typical current pro®le for standard spark plugs, however, we are currently working on a time-dependent computation of the current as a function of circuit parameters (impedance, resistance, applied voltage) and gap parameters (area, resistance) for our future simulations. Regardless of what method we choose to predict the spark current pro®le, one needs plasma conductivity (r) in our magnetic di€usion equation. The conductivity of the spark gap is a function of the background plasma [3] and will have to be calculated at each time step as the gas conditions change. Once J and E are known, then the electrical heat generation can be found using J
E. This heat term appears as a source term in the ¯ow energy equation [21]. Another exchange O. Yasßar / Parallel Computing 27 (2001) 179±200 188 term between the EM and ¯ow ®elds is 1c J  B force term that appears in the ¯ow momentum equation [21]. The electrical resistivity (g) presently used in our model is from classical transport theory [27,30]. It includes both the e€ects of Coulomb collisions and electron-atom collisions. The expression used is [31] g ˆ gW ‡ gF ; 11† where gW is (Weak) Coulomb collision, and gF is (Full) electron±atom collision term as gF  5:799  10 15 †
Z T †
ln K
T 3=2 and gW  me1=2 rea T † 1=2 T : e2 a T † Here, rea is the electron±atom collision cross section, me is the electron mass, e is the electron charge, Z is the average ionic charge state, T is the plasma temperature, and a is the degree of ionization of the plasma. The degree of ionization is a strong dependent of the plasma temperature. When the degree of ionization falls below 0:001, the behavior of g is controlled by gW . For a greater than 0:1, the dominant term is gF . Most combustion codes have a limited number of reaction rates based on the temperature range they operate in. In real spark-ignited (SI) engines, the electrical discharge creates a very high-temperature (tens of thousands of Kelvins) plasma kernel [5,23] which develops into a self-sustaining and propagating ¯ame front (a thin reaction sheet where the exothermic combustion reactions occur). We have extended the equation of state (EOS) data in KIVA-3 to temperatures beyond 5000 K by extrapolation of enthalpy data. However, this was a quick ®x and it needs to be further examined. We are planning to use the EOS data based on ionization e€ects at such high temperatures. We have already done such EOS computations for fusion plasmas and expect to have a similar approach [24]. Heat losses to the electrodes is an important factor to consider. However, spark plug mesh generation e€ort is not a trivial task as one has to guarantee smoothness [10]. Our results in this paper do not include a separate grid for electrodes, but we are currently planning to integrate a mesh for electrodes into the overall chamber. 3. Simulations The spark-ignition model we described here has been applied to a test case that is distributed with the standard KIVA-3 code [17]. We examine the combustion output (i.e., temperature, NOx density) against spark current (energy) and fuel density variations. The qualitative nature of an engine response against varying degrees of spark energy input and air/fuel ratio has been predicted, but a direct quantitative corre- O. Yasßar / Parallel Computing 27 (2001) 179±200 189 lation between engine output and spark discharge parameters has never been done. Early ¯ame formation and propagation has a vital impact on the engine performance and emissions. The energy input scenerio to the system determines the combustion eciency. While lean combustion is very desirable from an eciency and emissions standpoint, it is constrained by the potential for partial burns and mis®res that result in poor performance and higher (unburned hydrocarbon) emissions. By operating with lean mixtures, the eciency and speci®c fuel consumption of the engine can be increased, while reducing NOx and other toxic emissions. However, the decrease in the fuel/air ratio is limited by the increasing cyclic variability which results in some cycles either not burning completely or not igniting at all. Below are the results of our ignition-enhanced KIVA-3 for a test engine. The basic geometry of this reciprocating internal combustion engine can be seen in Fig. 1. The measurement of a spark discharge current represents a signi®cant task. For our simulations, we use Fig. 2, measured by our colleagues Sauers et al. [2]. This does not necessarily represent the distinct phases illustrated in [25] for a much expected behaviour of a spark discharge. Its slow rate of current rise and fall can be categorized as an TCI (energy stored in the inductance of a coil, that is essentially a glow discharge). Although there are many variables to test with our SI model, including the e€ect of di€erent phases of sparking event (breakdown, arc, glow phases and the duration of them), here we only examine how changes in a typical glow discharge current peak value a€ect the engine dynamics. Our future work will involve a close examination of further re®ned current measurements and the e€ect of each sparking phase. Fig. 1. Schematic of an engine. 190 O. Yasßar / Parallel Computing 27 (2001) 179±200 Fig. 2. A representation of spark discharge current pro®le. Besides varying the current peak value and the fuel amount, we have examined the in¯uence of ¯uid velocity ®eld on the spark energy deposition. It has been experimentally shown by zur Loye and Bracco [11] that ¯ow velocity ®eld a€ects ¯ame kernels signi®cantly. Our e€ort of examining the velocity e€ects also relates to previous modeling studies, such as by Kravchik et al. [9], where several assumptions have been made on the spatial and temporal gradients of the EM ®elds. It was assumed in Kravchik [9] that the current density was J ˆ rE. This is true for a system at rest. For a medium moving with velocity u relative to the laboratory, Ohm's law assumes the form J ˆ r E ‡ 1cu  B†. We have run cases with and without the second term that involves the velocity ®eld u and B to examine the di€erence. A signi®cant di€erence has been noted between two cases, which is consistent with observations made by zur Loye and Bracco [11]. Again, these results have been obtained on a baseline engine with a spark plug centrally located at the cylinder head. Temporal and spatial data was obtained at every 5° crank angle for a computational mesh of approximately 5500 elements. The engine speci®cations can be found in [17]. The enormity of the data makes visualization a necessity for scienti®c interpretation. The computational model produces 1.6±5.0 Gigabytes of data per piston cycle, depending on the computed parameters and based on granularity of spatial and temporal output. We had to convert the standard KIVA-3 output to other formats for post-processing purposes. One of two post-processing packages we use creates 3-D graphs and the other package (a virtual environment) is used to render these customized 3-D visualizations at immersive speeds, and navigate interactively. This virtual environment has proven valuable in understanding the relationships of computed parameters and the time-transient effects of changes in initial conditions. Viewing the 3-D simulation results in a dynamic way produces a smooth motion of the piston and a better understanding of the O. Yasßar / Parallel Computing 27 (2001) 179±200 191 compression and expansion strokes of the engine as a whole. The availability of temporal and spatial distribution of all physical quantities (density, pressure, velocity ®eld) via a virtual engine and user's access to represented data within seconds timeframe shortens the process of understanding the observed phenomena. Due to our limited capability to present these color graphs here, we include only a summary of our observations via Figs. 3±9. In these ®gures, we present peak values of temperature (T), NOx density, total density, fuel density at the ignition region, and Fig. 3. Fuel density at the ignition center. Fig. 4. Temperature at the ignition center. 192 O. Yasßar / Parallel Computing 27 (2001) 179±200 Fig. 5. Pressure at the ignition center and the cylinder wall. Fig. 6. NOx density at the ignition center. the expansion rates of certain temperature (1600 K) and NOx 1  10 5 g=cm3 † zones. Although the rate of expansion of pre-de®ned T and NOx zones does not exactly tell us the total amount of heat or NOx generated, it does tell about the dynamics of the mechanical and chemical processes. Figs. 3±9 are given for four di€erent cases. Each case is represented on the legend with several numbers indicating the discharge current (peak value) in Amps and the total amount of fuel injected in milligrams (mg). Two cases have the same discharge O. Yasßar / Parallel Computing 27 (2001) 179±200 193 Fig. 7. Total density at ignition center. Fig. 8. Radius of a radially expanding temperature zone (1600 Kelvins). current and fuel amount but slightly di€erent physics as J ˆ rE for case-A and J ˆ r E ‡ 1cu  B† for case-B. The fuel injection and ignition processes are started before the end of the compression stroke in an e€ort to maximize the air-fuel mixing and to enhance turbulence e€ects on ¯ame propagation. Fuel injection takes place from 52° to 40° Crank Angle (CA). Fig. 3 shows fuel density at the ignition center. In all cases, the fuel density is about the same up to the time of the ignition, which lasts from 27° to 194 O. Yasßar / Parallel Computing 27 (2001) 179±200 Fig. 9. Radius of a radially expanding NOx density zone (1  10 5 g=cm3 ). approximately 17° CA. The fuel depletion rate and amounts are di€erent for all cases. The energy deposition and subsequent chemical energy release and heat diffusion starts as soon as the spark discharge is applied, sending a heat wave away from the ignition center. As the piston continues to move up, the pressure di€erence between the ignition center and the cylinder wall goes up slightly (Fig. 5), causing a ¯ow towards the center of the chamber. This di€erence is due to the existence of a bowl at the center of the piston head. It is this pressure di€erence that slows down, stops and reverses radially outward heat and NOx waves. However, as we reach the top dead center (TDC) and the subsequent expansion stroke, the energy and species concentrations are released to the rest of the chamber. The energy deposition rate and amount are di€erent between four cases presented here. The case with the highest current peak value (1.1 A) deposits more energy than others, causing a higher peak temperature at the ignition center. EM and mechanical properties play an important role during the early ignition phase whereas chemical processes are not as signi®cant in the early stages as they are later. The di€erence in temperature between cases with di€erent fuel densities is only around several hundred degrees (in Kelvin) despite a notable di€erence in their fuel amounts. Towards the end the compression stroke, the ¯uid velocity is mostly in the radial direction (^r) in the ignition region. Since the magnetic ®eld lines ¯ow in the azimuthal direction for a discharge current in the z direction, the induced e€ects (u  B) lead to an increase in the current density (J). This increases the amount of EM heat (J
E) for case-B in comparison to case-A where the induced e€ects are assumed negligible. The di€erence in the peak temperature between case-A and case-B due to the interaction of the ¯uid velocity ®eld and the magnetic ®eld seems quite large during the initial stages of ignition where the discharge current is relatively high. In fact, case-A has the least amount of energy deposition. O. Yasßar / Parallel Computing 27 (2001) 179±200 195 After spark ignition has ended, the dominant processes are chemistry, di€usion and mechanical expansion of the early ¯ame. One particular observation is that the temperature pro®le (Fig. 4) for case-A crosses others and stays at a slightly higher value for the rest of the engine cycle. We believe this is because case-A has a slightly higher gas and fuel density (which continues to burn) at the ignition center due to a slower expansion of the early spark ¯ame propagation. Another density e€ect on the temperature is for case-D where it drops more signi®cantly than others around 10° CA. Total density for case-D is lower than others (naturally due to lower amount of fuel spray) and as soon as the piston moves into the expansion stroke, this di€erence becomes more notable and leads to a higher drop in the temperature at the center. The level of NOx is strongly dependent on the peak temperature. Due to less amount of electrical heat generation in case-A, the total amount of NOx generation (Fig. 6) is lower in comparison to case-B. The heat di€usion and convection in case-B is also stronger due to higher heat generation and the increased radius of both the heat and the NOx waves. The di€erence in the heat and mechanical expansion rates (case-B expands more due to higher energy deposition to start with) a€ect the spatial variations in the NOx density. The peak value of NOx density, as plotted in Fig. 6, occurs at the ignition region where energy deposition and fuel injection take place. The density goes up as the compression stroke approaches the TDC. Notice that NOx density at the cylinder center is slightly higher for case-A than case-B. We attribute this to less expanded NOx zone as seen in Fig. 9 as well as to a slightly higher density at ignition center in comparison to case-B. The NOx density pro®le at the ignition center re¯ects e€ects of compression stroke as well as those of electrical, chemical ones and even of early ¯ame mechanical processes. Although these pro®les show peak values at a certain point (ignition center), they are somehow an indicator of the total NOx amount in the activated volume, except for case-A. Changes in fuel density produce noteworthy changes in the system. Less fuel produces a lower peak temperature, thus NOx production decreases with less fuel usage. Another observation is the distribution of the fuel in the chamber. Injection dynamics and the mechanical expansion of the early ¯ame propagation and the compression stroke change the ®eld density in the chamber as intended. Fuel density at the ignition center is shown in Fig. 3. There is an increase between 52° and 40° CA due to injection. The injected fuel moves from the injection region (center) to the rest of the chamber and there is a decrease at the center due to that between 40° and 30° CA, just before the ignition. A more signi®cant drop occurs during the ignition due to burning. A typical trend here is that fuel depletion is proportional to the amount of spark energy input. The e€ect of chemical processes from 30° CA and on, up to the TDC, is well seen here. The fuel density pro®le would have followed a similar trend as the total density (an increase as one approaches TDC), but instead it ¯attens out and then goes to zero. The compression stroke replaces fuel that is burned, maintaining a ¯at slope in fuel density at the ignition center. Eventually, fuel depletion is complete, occurring at di€erent CAs for each case, depending on the spark input energy. An important observation is made here about the case with lower fuel input (case-D). Although this case has a lower amount of fuel to start with, the injected fuel does not move to the rest of the chamber as fast as the other 196 O. Yasßar / Parallel Computing 27 (2001) 179±200 cases because of a slower expansion of activated volume whose acceleration depend strongly on the amount of early energy deposition. The higher the energy deposition at the early ¯ame propagation times, the further the leftover fuel is pushed away. The concentrated fuel region moves down to the bottom of the piston bowl as it gets burned. Less fuel in the system produces less resistance to the mechanical expansion of the ¯ame. The mechanical conditions of a chamber are known to greatly a€ect the ¯ame propagation. The ¯ame propagation distance is directly related to air/fuel ratio as well as the pressure and physical length scales of the system. Yet, as compression dynamics becomes dominant and pushes activated volume back to the center, the fuel density in case-D at the ignition center falls below its counterpart (case-B). Due to advance ignition timing (before the TDC) and the geometry of the combustion chamber, an interesting phenomena occurs with the mechanical expansion of the ¯ames. It is a common practice to start injection and ignition before TDC to increase fuel mixing and cylinder pressure before the expansion stroke. We have selected a certain contour for temperature (1600 K) and NOx (1  10 5 g=cm3 ) to detect di€erences between cases presented here. Since the T, NOx , and fuel density pro®les at the ignition center are not indicative of the total amount of heat and NOx formation, perhaps the expansion zones can provide additional information. Differences in ¯ame dynamics (activated volume) before and after TDC are notable. Despite di€erent expansion rates during the early phase of ignition, the ¯ame kernel submits to the same con®nement rate during the late compression and early expansion phases when the pressure has rised in all parts of the cylinder due to piston motion. This phenomena might have occurred quite di€erently if the piston was completely ¯at. The existence of a bowl in the chamber causes variations in the pressure throughout the cylinder. Around 15° after TDC, however, the trapped heat and species concentration release at di€erent rates due to di€erent amount of energy stored inside the wave. The cylinder pressure and total density is a major factor determining this behaviour. The amount of heat trapped determines how much further the wave expands and how early it loosens up. Despite a comparable spark energy input, chemical energy release is lower with case-D and thus the ¯ame propagation is slower. Other cases have the same amount of fuel (11.6 mg), but di€erent amount of spark energy input, which leads to di€erent expansion rates for activated volume in each case. Although we have varied the peak current value just slightly (10%), this produces signi®cantly di€erent computational results. Higher current produces more energy in the combustion chamber, resulting in higher temperatures which, in turn, produce more NOx . As more energy is deposited into the background gas, the rate of ¯ame expansion increases. We have attached contours of temperature and NOx pro®les for both case-B (1.0 A/11.6 mg) and case-C (1 A/11.6 mg) cases here in Figs. 10±13. Due to cost and space limitations only a sample of the images is presented. While an attempt was made to select typical images, these images are not representative of what takes place throughout the engine cycle. Additional data can be found by contacting the principal author. These ®gures show the 1600 K temperature and 3 10 5 g=cm NOx density contours. The spatial variations within the activated volume 3 are not well represented due to the maximum limits (1600 K, 10 5 g=cm ) we have O. Yasßar / Parallel Computing 27 (2001) 179±200 197 Fig. 10. A centrally expanding temperature zone (T P 1600 K): Spark current peak is 1.1 A. Fig. 11. A centrally expanding temperature zone (T P 1600 K): Spark current peak is 1 A. Fig. 12. A centrally expanding NOx density zone (P1  10 5 g=cm3 ): Spark current peak is 1.1 A. imposed on them. Anything above these values is represented in one color. The maximum temperature at this crank angle actually is well above 6500 K and the maximum NOx value is on the order of 10 4 g=cm3 . It is quite clear from these graphs and our visualizations that both heat and NOx zones in the higher current O. Yasßar / Parallel Computing 27 (2001) 179±200 198 Fig. 13. A centrally expanding NOx density zone (P1  10 5 g=cm3 ): Spark current peak is 1 A. case expands more in all directions. This is desirable for having a more complete and uniform combustion, but it will produce more NOx . The peak value of NOx production, the amount and expansion rate of the NOx zone, and fuel consumption are higher for the higher current case as well. 4. Conclusion A detailed model with time- and spatially-dependent characteristics has been developed for spark ignition in conjunction with engine combustion and ¯uid ¯ow. This model has been integrated into the KIVA-3 engine model. Contrary to several assumptions in previous modeling e€orts, temporal and spatial variations of the spark related EM ®elds are important source of in¯uence on the ¯ame dynamics and subsequent combustion. A virtual engine environment has been developed for timely scienti®c and technical interpretations of results. A spark-control mechanism has been developed by our colleagues [2] to help study spark ignition in more details, validate model, and develop a breakdown control mechanism for real engine conditions. The spark-control mechanism would not only reduce the voltage variations but also would help bring down the breakdown voltage value to more acceptable levels that could minimize electrode wear and the level of emissions. All these new developments are expected to help design better engines with higher eciency, lower emissions and more consistent combustion. A direct correlation between combustion output and the spark discharge current has been established through numerical simulations. Conditions favoring a low discharge current that might provide enough energy for a complete combustion need to be found for a set of design parameters. Our focus has been on the accuracy of spark energy deposition so far. Our techniques will be applied to a real engine and spark plug design with more speci®cations on the plug under usage. The heat loss through electrodes and the e€ect of the spark geometry on the ¯ame propagation will be investigated. E€ect of the ¯uid velocity will be investigated in more details. E€ect of 1c J  B force on the shape of spark kernel will be investigated. Among future modeling tasks are an improved O. Yasßar / Parallel Computing 27 (2001) 179±200 199 plasma resistivity model, an expanded EOS and chemical database for high temperatures as well as mesh generation and ®ne-grain problems on both symmetric shared-memory and distributed-memory multiprocessor systems. The massively parallel version of KIVA-3 will be used in solving large-scale problems on scalable systems. We plan to explore the use of in-house distributed-computing tools to facilitate dynamic data passing between our combustion model and post-processor. These tools also increase the possibilities for a collaborative environment where researchers located at remote locations can view and steer the simulations through a secure network. Integration of our spark ignition model into our previously developed massively parallel KIVA-3 as well as development of a portable version of parallel KIVA-3 have been successful. Comparisons with real experimental data is necessary to validate the model. Integration of our spark-control mechanism into real engine conditions is also being investigated. Acknowledgements This work was sponsored by the Energy Eciency and Renewable Energy Program and the Energy Research Oce of the US Department of Energy, under contract DE-AC05-84OR21400 with Lockheed Martin Research Corporation. Computing resources and experimental equipment were made available, respectively, through the Center for Computational Sciences, Advanced Propulsion Technology Center and Life Sciences Division of the Oak Ridge National Laboratory. References [1] N. Wiegart et al., Inhomogeneous ®eld breakdown in GIS: the prediction of breakdown probabilities and voltages, IEEE Trans. Power Delivery 3 (3) (1988) 931±938. [2] Sauers, D.D. Paul, J.W. Halliwell, C.W. Sohns, Energy delivery test measurement system using an automotive ignition coil for spark plug wire evaluation, Patent Disclosure ESID No. 1774XC, S-83,831. [3] O. Yasar, A computational model for Z-pinch plasma channels, Ph.D. Thesis, University of Wisconsin-Madison, December 1989. [4] H. Albrecht et al., New Aspects of Spark Ignition, SAE 770853. [5] R. Maly, Spark ignition: its physics and e€ect on the internal combustion engine, in: J.H. Hilliard, G.S. Springer (Eds.), Fuel Economy, Plenum Press, New York, 1984. [6] J. Lu, A.K. Gupta, E.L. Keating, E€ect of IC engine operating conditions on combustion and emission characteristics, Trans. ASME 115 (1993) 694. [7] D. Bradley, F.K.K. Lung, Spark ignition and the early stages of turbulent ¯ame propagation, Combust. Flame 69 (1987) 71. [8] R. Herweg, R.R. Maly, A fundamental model for ¯ame kernel formation in SI engines, SAE Paper No. 922243, 1992. [9] T. Kravchik, E. Sher, J.B. Heywood, From spark ignition to ¯ame initiation, Combust. Sci. Technol. 108 (1995) 1±30. [10] T. Mantel, Three dimensional study of ¯ame kernel formation around a spark plug, SAE Paper No. 920587, 1992. 200 O. Yasßar / Parallel Computing 27 (2001) 179±200 [11] A.O. zur Loye, F.V. Bracco, Two-dimensional visualization of ignition kernels in an IC engine, Combust. Flame 69 (1987) 59±69. [12] M. Nakai et al., Stabilized Combustion in a spark ignited engine through a long spark duration, SAE 850075. [13] R.M. Moiden, M.D. Checkel, J.D. Dale, The e€ect of enhanced ignition systems on early ¯ame development in quiescent and turbulent conditions, SAE 910564. [14] M.S. Hancock, D.J. Buckingham, M.R. Belmont, The in¯uence of arc parameters on combustion in a spark-ignition engine, SAE 860321. [15] D.J. Fitzgerald, Pulsed plasma ignitor for internal combustion engines, SAE 760764. [16] Hall et al., Initial studies of a new ignitor, SAE 912319. [17] A.A. Amsden, KIVA-3: A KIVA Program with block-structured mesh for complex geometries, Los Alamos Technical Report LA-12503-MS, 1993. [18] O. Yasar, T. Zacharia, Distributed implementation of KIVA-3 on the intel paragon, in: S. Taylor (Ed.), Parallel Comput. Fluid Dyn., North-Holland, Amsterdam, 1995. [19] O. Yasar, G.A. Moses, R.R. Peterson, Plasma channels for the LIBRA reactor designs, Fus. Technol. 19 (1991) 669. [20] S.C. Brown, Introduction to Electrical Discharges in Gases, Wiley, New York, 1966. [21] O. Yasar (Ed.), A scalable algorithm for chemically reactive ¯ows (special issue on Advanced Computing on Intel Architectures), J. Comput. Math. Appl. 35 (7) 1998. [22] W.F. Hughes, F.J. Young, The Electromagnetodynamics of Fluids, Wiley, New York, 1966. [23] J.B. Heywood, Internal Combustion Engine Fundamentals, McGraw-Hill, New York, 1988. [24] J.J. MacFarlane, P. Wang, O. Yasar, Non-LTE radiation from micro-®reballs in ICF target explosions, Bull. Amer. Phys. Soc. 34 (1989) 2151. [25] R. Stone, Introduction to internal combustion engines, SAE Publications (1993) 150±153. [26] O. Yasar, G.A. Moses, Z-pinch plasma channels for the LIBRA reactor design, Nucl. Fus. 31 (1991) 273. [27] O. Yasar, G.A. Moses, Explicit adaptive grid radiation magnetohydrodynamics, J. Comput. Phys. 100 (1) (1992) 38. [28] A.A. Amsden, KIVA-II: A computer program for chemically reactive ¯ows with sprays, Los Alamos Technical Report LA-11560-MS, 1989. [29] J.D. Jackson, Classical Electrodynamics, Wiley, New York, 1975. [30] J.J. Watrous, G.A. Moses, R.R. Peterson, Z-PINCH ± A multifrequency radiative transfer magnetohydrodynamics computer code, Technical Report UWFDM 584, University of WisconsinMadison, Fusion Technology Institute, 1985. [31] S.V. Dresvin (Ed.), Physics and Technology of Low Temperature Plasmas, Iowa State University Press, Ames, IA, 1977.