Computer simulation of combustion in a stratified charge engine

C O M P U T E R S I M U L A T I O N OF C O M B U S T I O N IN A S T R A T I F I E D CHARGE E N G I N E A. A. BONI, M. CHAPMAN, J. L. COOK, AND G. P. SCHNEYER Science Applications, Inc., P.O. Box 2351, La Jolla, California 92037 A two-dimensional (axisymmetric), time-dependent simulation of the compression and power strokes of a divided-chamber stratified charge engine based on an axisymmetric approximation to the Honda CVCC configuration has been performed. Global methane/air kinetics and constant turbulent diffusivities have been employed. Results indicate that significant variations in the CH4/air equivalence ratio occur near the spark plug as a result of the velocity field induced by the moving piston. These variations could cause cycle-by-cycle variations in engine performance. Vortex-like motion induced by the piston is presented along with the spatial distribution of ca4, products, NO~, temperature and velocity at various times in the engine cycle. Unburned CH 4 is found in the prechamber at the end of the calculation.

1. I n t r o d u c t i o n Better methods of a n a l y z i n g c o m b u s t i o n d y n a m i c s in internal c o m b u s t i o n (IC) engines are needed due to increasing d e m a n d s for both greater thermal efficiency and lower pollutant output. Consequently, various improved diagnostic devices and analytical techniques are b e i n g developed. 1 The stratified charge engine (SCE), w h i c h has the potential of a c h i e v i n g both objectives, has recently received r e n e w e d attention, 2 It is apparent that large spatial and temporal gradients of temperature, pressure, velocity, and species concentrations exist in such a device. F l a m e p r o p a g a t i o n t h r o u g h a stratified charge depends to a great extent on the hydrod y n a m i c mixing and transport processes, on the chemical reaction m e c h a n i s m s and rates, a n d on the pressure, d e n s i t y and temperature, all of which vary spatially and temporally w i t h i n the engine cylinder. Consequently, a mathematical m o d e l c a p a b l e of d e s c r i b i n g flame propagation t h r o u g h a turbulent, stratified fluid is n e e d e d to accurately determine b u r n i n g (energy release) as well as pollutant p r o d u c t i o n rates. Most models of IC engine combustion are t h e r m o d y n a m i c models, i.e., the combustion c h a m b e r is considered a control volume, the fluid d y n a m i c m o t i o n is neglected and the energy release rate specified. 3-s These models

have been used successfully for many applications. However, fluid mechanical processes are m u c h more important in engines such as the SCE, so d e v e l o p m e n t of spatial models w h i c h include such effects are necessary. Sirignano, 6 Rosentweig-Bellan and Sirignano, v Bracco 8'9 and Boni, et al lOhave presented one-dimensional, t i m e - d e p e n d e n t simulation models a p p l i e d to conventional Otto cycle, Wankel and stratified charge engines. Most of this work e m p l o y e d algebraic turbulence models, assumed gas phase combustion, and utilized global chemical kinetics formulations. Bracco included two-phase droplet combustion n and recently presented some preliminary two-dimensional (2D) results as well. ~2 Boni, et al ~3 first reported on the 2D aspects of the d i v i d e d chamber, S C E flow field. Their calculations w h i c h ignored m a i n chamber comb u s t i o n indicated the existence of a separated, vortex flow in the m a i n chamber. T h e purpose of the present paper is to extend our earlier work 10.13b y s t u d y i n g the compression and power strokes with a time-dependent, 2D (axisymmetric) computer simulation model. In the following paragraphs, the numerical simulation m o d e l for a d i v i d e d chamber, SCE is described. The compressible conservation equations with constant turbulent diffusivities c o u p l e d to a simple, global chemical reaction model have been formulated for solution using an arbitrary Lagrangian-

1527

1528

MATHEMATICAL MODELING

Eulerian (ALE) numerical technique. A simplified NO c h e m i c a l kinetics model is used to calculate NO~ concentrations. The simulation yields a p r e d i c t i o n of the flame propagation dynamics l e a d i n g to temperatures, pressures and species concentrations as functions of spatial location, time, and various engine parameters. 2. A p p r o a c h We obtain a solution of the conservation equations b y use of the A L E technique. T h e method uses a finite-difference mesh with grid points w h i c h m a y move with the fluid (Lagrangian), be h e l d fixed in space (Eulerian), or be moved in any other prescribed m a n n e r ) 4 An implicit f o r m u l a t i o n permitting accurate solutions to be o b t a i n e d for flows at all speeds is also incorporated. The governing equations include conservation of mass, m o m e n t u m , energy and c h e m i c a l species, as well as thermal and caloric equations of state. Initial a n d b o u n d a r y conditions must also be specified. In the A L E approach, it is convenient to express the conservation equations in integral form for an infinitesimal (finite-difference grid element) volmne w h i c h m a y be in motion with arbitrary, p r e s c r i b e d v e l o c i t y ) ~ Denoting the volume e l e m e n t b y V(t), the surface of V(t) b y S(t), and the outward normal vector b y fi, the equations expressing conservation of mass, m o m e n t u m , and energy, respectively, are

pdV-

Here, p, ~i and e T are the fluid density, velocity and specific total energy, respectively; 0 is the grid velocity; P the pressure; K t is the turbulent thermal conductivity; T t h e temperature; e i(T) the specific internal energy of species i i n c l u d i n g the energy of formation; D, the multi-species diffusion constant of species i; and F~ the mass fraction of species i. ~ is the viscous stress tensor defined in terms of the laminar a n d t u r b u l e n t viscosity coefficients, ix I and ixt, respectively, as ~- (Ix l + Ixt) def ti

Conservation of the ith chemical species is given b y the e q u a t i o n

~

dV-fs.>PF,(O-~')r*dS

r

] . (oD~VF~)dS (5)

where 5

E

Fi= 1

and

p,=pF i

We consider five species in our formulation, viz, fuel (CH4) , oxidizer (02) , diluent (N2) , products and p o l l u t a n t (NO). We use an ideal gas equation of state 5

(t)

Pi= o i B T / M ,

P= E

(1)

of Ot

o(,dv- f

o~,(O-~,)

of

v(t)

Pi

(7)

,=1 with R being the gas constant, and

ridS

s(t)

v(t)

+f Ot

(6)

i=l

p ( ~ r - ~) 9 ridS = 0

(t)

(4)

PeTdV- f

+f

S(t)

aPdS=f PeT(U

s(t)

(fi . u) PdS = f s(t)

ei(T) = tl, dS

(2)

S(t)

fi . f i , d S + S(t)

(3)

+f

ri 9 ( K t V T ) d S S(t)

'

+ e,(29S.15~

~) ridS

S(t)

C~,(T)dT 298.15OK

(8)

Values for Cv.i(T) a n d ei(298.15~ the constant volume specific heat and energy of formation, respectively, for each species, are taken from the c o m p i l a t i o n of G o r d o n and McBride 15 or the J A N A F T a b l e s ) 6 Here we are p r i m a r i l y interested in illustrating the fluid mechanical features and their inter-relationship to the combustion processes; hence, we m a i n t a i n a simple chemical model. As shown b y Bowman, 17 it is possible, b y a judicious choice of rate constant, to o b t a i n a

COMBUSTION IN A CHARGE E N G I N E t e m p e r a t u r e p r o f i l e c l o s e l y a p p r o x i m a t i n g that o b t a i n e d w i t h a d e t a i l e d , multi-step, k i n e t i c c o m b u s t i o n m e c h a n i s m . T h u s , we l i m i t our s t u d y to the m e t h a n e / a i r system for w h i c h t h e p a r t i a l - o x i d i z a t i o n r e a c t i o n occurs via C H 4 + 1 / 2 O 2 ~ C O + 2 H 2 as t h e r a t e - c o n t r o l l i n g step in the c o m p l e t e c o m b u s t i o n of C H 4. F o r this reaction m e c h a n i s m , w e use the data s u m m a r i z e d b y M e l l o r ~8 to o b t a i n for the rate of C H 4 c o m b u s t i o n (in cgs-~ units) 6 x 10 s M C'/2 T H4 (PCH4)clae m =

-

-

12,200

p Cl/2 ,, exp H 4 I-s0 2

(9)

T

F o r the N O kinetics, w e utilize the simplified Zeldovich mechanism. k1

O+N

TABLE I Pertinent calculation parameters Engine Speed = 6000 RPM Bore = 7.400 cm

Stroke = 8.636 cm Rod Length = 12.818 cm K o = 14.06 exp

M o 2 ( 1 0 - 6 p)o.815

2~NO+N

(10)

k4

( - l15,600/ RT)

kl = 7.14 x 1013exp ( - 75,500 / R T) kz = 1.20 x 10 l~ T e x p (-5,500/RT)

N + O 2 ~ NO + O

B.D.C. Main Chamber Equivalence Ratio (qb) = 0.62 B.D.C. Prechamber Equivalence Ratio (~b) = 1.51 B.D.C. Temperature = 300~ B.D.C. Pressure = I Bar Ignition Crank Angle = 327 ~ (33 ~ BTDC) Turbulent Momentum Diffusivity = 50 cm 2 / see Turbulent Prandtl / Schmidt Numbers = 0.70

k a = 3.20 x 1 0 9 T e x p

(- 39,100/ HT) k4 =

k~

1529

1.50

><

10 la

(11)

k3

It can be s h o w n , that if d [ N ] / d t [ 0 2 ] - [O] are in e q u i l i b r i u m ,

= 0 and

d [NO] -

dt

- 2K~/2 k , [ 0 2 ] ,/2 I N 2 ]

[NO] 2

k ak4

k 4 [NO]

(12)

I + k2 [ 0 2 ]

Therefore, d [NO] ( b N O ) c , heil/ = M

N

O

-

dt

(13)

w h e r e MNO is the m o l e c u l a t w e i g h t of N O a n d K o is the e q u i l i b r i u m c o n s t a n t for the o x y g e n d i s s o c i a t i o n reaction. E q u a t i o n s (9) a n d (13) are the c h e m i c a l rate expressions s o l v e d s i m u l t a n e o u s l y w i t h Eq. (5) for [ C H 4 ] a n d [ N O ] d i s t r i b u t i o n s t h r o u g h o u t the e n g i n e cycle. Values for local [ 0 2 ] a n d [N2] are o b t a i n e d f r o m s t o i c h i o m e t r y a n d a law of mass action for the m e t h a n e o x i d a t i o n reaction. T h e r e a c t i o n rate c o n s t a n t s e m p l o y e d are reprod u c e d in T a b l e I. I n the p r e s e n t paper, the species, e n e r g y a n d

viscous (momentum) turbulent diffusivities are c o n s t a n t in s p a c e a n d t i m e e.f. T a b l e I. W e c h o o s e a " r e p r e s e n t a t i v e " v a l u e for the v i s c o u s d i f f u s i v i t y b a s e d o n the e x p e r i m e n t a l results of Witze, for a m o t o r e d s i n g l e - c y l i n d e r e n g i n e 19-2' o p e r a t i n g at 2 0 0 0 RPM. Similar m e a s u r e m e n t s , also on m o t o r e d e n g i n e s , h a v e b e e n reported. 22-2v We are u n a w a r e of meas u r e m e n t s s u b s e q u e n t to i g n i t i o n or at h i g h e n g i n e R P M . T h e use of m o t o r e d e n g i n e data, h o w e v e r , m a y be a p p r o p r i a t e in l i g h t of the r e c e n t h y p o t h e s i s of A n d r e w s , et al. 2s that t h e ratio of t u r b u l e n t to l a m i n a r b u r n i n g v e l o c i t y correlates w e l l w i t h the R e y n o l d s r m m b e r (based on the t u r b u l e n t microscale) of the c o l d gas into w h i c h a f l a m e propagates. W i t z e z~ observes, for all crank angles, n e a r l y c o n s t a n t m a g n i t u d e s for the relative t u r b u l e n c e i n t e n s i t y (u'), m i c r o s c a l e l e n g t h (k) a n d m a c r o s c a l e l e n g t h (A). T h e m e a n v e l o c i t y (U) increases s l i g h t l y d u r i n g intake a n d near top d e a d center a n d b o t t o m d e a d center. It is not k n o w n w h e t h e r this b e h a v i o r is typical of all spatial locations or just near t h e c y l i n d e r d o m e as m e a s u r e d b y Witze. N e v e r t h e l e s s , f r o m the e x p e r i m e n t a l d a t a r e p o r t e d in Ref. 20, and in s u b s e q u e n t m e a s u r e m e n t s , 21 we can o b t a i n s o m e o r d e r of m a g n i t u d e estimates for t h e diffusivities. W e e m p l o y a c o n v e n t i o n a l d e f i n i t i o n of e d d y d i f f u s i v i t y , vt, in terms of l e n g t h scale (l) a n d the t u r b u l e n t k i n e t i c e n e r g y (k = 0.5 u '2 ), viz

1530

MATHEMATICAL MODELING vt

~L t -

40% higher. A predictive turbulence m o d el is required to establish the validity of this approximation. More extensive experimental data on t u r b u l e n c e levels in motored and b u r n i n g IC engines is also needed. Th e equations can now be solved for the selected engine geometry and initial conditions, c.f., Section 3, provided that the b o u n d ary conditions are specified. We e m p l o y a no-slip condition on velocity and no-flux condition on mass and energy at the b o u n daries of the cylinder. Furthermore, the piston face moves at a prescribed t i m e - d e p e n d e n t velocity c o r r e s p o n d i n g to the desired speed.

(14)

C~. lk 1/2

-

P

C is a constant of order 0.1. 29 From Reference 20, u ' / ( Y ~- 0.5, 0 = 1000 e m / s e c , h = 0.05

cm, and A = 0.5 cm. Subsequent measurements by Witze 2~ indicate that these values scale linearly with e n g i n e RPM. We thus find that for 6000 RPM, 1 c m 2 / s e c < v t < 100 c m 2 / s e c is a reasonable order of magnitude bound. Since the t u r b u l e n t Prandtl and S c h m i d t numbers ~ 0.7, 29 we take v t = 50 c m 2 / s e c and the e n e r g y / s p e c i e s diffusivities about

J 14

i0

12

6

8

4

2

0

(a)

/ R(CM)

J

/k

J

IIII |llI

J J J

Z(CM)

J f

f _

f

14

I0

12 (b)

1 J

14

12

iO

Q

6

4

2

0

(c)

F1c. 1. Computational grids; (a) main grid at t = 0 (0 = 180~ (b) main grid at t = 4.95 msec (0 = 358~ (c) fine grid at t = 0 (0 = 180~

COMBUSTION IN A CHARGE ENGINE

1531

3. Results Equations (1)-(13) are written in finite-difference form for a 2D axisyrnmetrie coordinate system. Two separate computational grids were employed, both of w h i c h are illustrated in Fig. 1. Figures l a a n d l b show the m a i n grid at 108 ~ (BDC) and 358 ~ ('2~ BTDC), the difference being a conformal axial compression of the main c h a m b e r grid. The opposite occurs d u r i n g expansion. Zones within the nozzle and p r e c h a m b e r are fixed (pure Eulerian). Thus, a constant n u m b e r of zones are in the computational d o m a i n throughout the calculation. Figure i c shows a more finely zoned grid at 180 ~ (BDC). The latter was used only for a compression stroke calculation to obtain better resolution of the p h e n o m e n a near the cylinder wall. T a b l e I lists the pertinent geometrical, transport, thermodynamic, and chemical constants, as well as the initial conditions e m p l o y e d in the calculations. At the b e g i n n i n g of the compression stroke (BDC) the piston accelerates into the m a i n chamber reaching, for example, velocities of 93 c m / s e c , 726 c m / s e c , and 1468 c m / s e c at times of 8.57 x 10 5 see (183~ 5.58 x 10 4 see (200 ~ and 1.22 x 10 3 see (224~ respectively. Pressure waves are generated at the piston face which compress the gas in the cylinder. A plot of the i n d u c e d velocity field for the fine-grid geometry of Fig. l c is illustrated in Fig. 2 for a time of 2.77 msec (279.71~ (Note that the headless vectors run a w a y from the crosses d e n o t i n g their origin; vector lengths are proportional to the c o m m o n logarithm of the speed.) We observe the formation of a complex, v o r t e x - l i k e flow field in the piston face-cylinder wall corner, but for a piston velocity-acceleration range well b e y o n d those e m p l o y e d in the works of Tabaczynski, et al., 3~ and Daneshyar, et al. 31 Indeed, the latter have noted that whereas at low engine speeds of ~ 6 0 0 RPM the vortex size may be significant with respect to the cylinder bore, the vortex size for engine speeds greater than 2400 RPM is c o n s i d e r a b l y smaller. Our calculations deal with the outer limit of the latter regime. At such high velocity and acceleration levels, the analysis of Ref. 3 i is invalidated due to neglect of the inertial terms in the m o m e n t u m equation and transition to turbulent flow. Our calculations show no clear, large-scale, recirculating vortex formation since lateral velocities only of order 10% of the longitudinal velocity are found. This behavior is illustrated in Fig. 2. The same behavior was found for both computational meshes.

VELOCITVECTORS Y 12 t~r lit~

~ lilt&

II

I

I I

10 ] l 11 11 1 1 11 11 111 I

I I I l~l.

I I 11111, 1 I I lla, I I lllll,

I I 111i~

s 41 1 1 1 1 z l

111 1 1 11 I 1 1 11 11 I 6 1111 I 1111 I

0

lll/a 1 1111~ I 1 Ilia 11111L I lllia I IIIIL

2

Ftc. 2. Velocity vector plot of fine grid case in middle of compression stroke (t = 2.77 msec, 0 = 280 ~) showing oscillating radial velocity component along main chamber wall.

The fine-mesh calculation e m p l o y e d 3 zones adjacent to the wall (of thickness 0.05 cm, 0.06 era, and 0.10 cm) to resolve the wall b o u n d a r y layer while the m a i n grid used a single 0.37 cm thick zone. Close examination of Fig. 2 also shows an oscillatory behavior of the radial velocity near the cylinder wall. Again, this behavior is obtained for both c o m p u t a t i o n a l meshes, c.f., Figs. 2 and 3a, a n d appears only d u r i n g the compression stroke. It is not clear whether this behavior is numerical in origin or is the result of some form of h y d r o d y n a m i c instability p o s s i b l y associated with the interaction of the vortex and the wall in the presence of compression waves e m a n a t i n g from the piston. A detailed study of this p h e n o m e n a is b e y o n d the scope of the present paper; however, we should point out that if the p h e n o m e n a is indeed physical, its effect on the q u e n c h layer behavior and the resultant u n b u r n e d hydrocarbon emissions c o u l d be significant. To our knowledge, such oscillatory behavior has not been observed or reported previously. Figure 3 shows a velocity vector plot (3a) and a C H 4 mass fraction contour plot (3b)

1532

MATHEMATICAL MODELING Z( cm ) 15

T > ~

,~ 7.6T~E-03 7.595E-02 3.59~E-03

DQ =

R( CFI )

13 ~

13

1 //.

1!

/.I-. /.

11-

1.

'lI

9.

)

2 (a)

2 (b)

FIr 3. Representative plots at ignition point (t = 4.09 insec, 0 = 327~ (a) velocity vector plot; (b) CH 4 mass fraction contour plot.

at 327.23 ~ (4.09 msec), the ignition point.* Flow velocities at this time range up to 3000 c m / s e c . The large velocity gradients coupled to the effect of species diffusion have resulted in a significant modification in the equivalence ratio distribution from that initially charged into the engine. Figure 4 shows the radial distribution of C H 4 mole fraction at the top of the prechamber for two crank angle positions around the time of ignition. Note the significant radial gradient and its change with time due to diffusion and the velocity field. These variations may be important in terms of spark plug location and ignition characteristics of SCE's a n d have been cited previously as major causes of cycle-by-cycle variations in automotive c o m b u s t i o n behavior, ae,aa I n our previous work, 1~ we effected ignition by i n s t a n t a n e o u s l y reacting all of the 0 2 * In all figures showing contours of variables, the highest (H) and/or lowest (L) values are indicated on the plot and as Q(H), Q(L). Uniform contour increments, DQ, can be employed to calculate intermediate values.

in the last prechamber zone. This method, hereafter termed the "fast b u r n , " was applied in the present calculation to the three innermost radial (0.30 cm radius) and two uppermost axial (0.76 cm thickness) zones of the computational grid. Figures 5a and 5b show the velocity vectors and c a 4 m a s s fraction contours at 4.55 msec (343.9 ~ just after u n c h o k i n g of the nozzle. It can be seen that the flow is streaming i n a jet-like m a n n e r into the m a i n chamber with ignition of the m a i n chamber charge b e g i n n i n g at about the 11 cm longitudinal position a n d just off centerline. To examine the importance of the ignition simulation t e c h n i q u e on the b u r n i n g characteristics, we also applied an alternate technique, termed "slow b u r n , " to the same zones at the same crank angle. In the slow b u r n technique, the chemical reaction rate in these few cells is set to be that at the real pressure and molar concentrations, but at a fictitious temperature approximately equal to the C H 4air upper ignition temperature (1025~ This artificiality continues until the real tempera-

C O M B U S T I O N IN A C H A R G E E N G I N E

0

[',,. 0,1 ~

9

1533

03 ~1 Cq

O

I

I

I O

I I

\

I

\ \

\

\

\

\

\ \

\

\ \

\

\

\

o ~z ~

\ \

\ \

\

\

\

\

\ \

r O

\

\

\

\

\

\

\ \ I

c~

\

-4

\

\

~2

\

\

\ \

\

\

\

\

\ \

\ I 0

4

I

i

r-4

\ O

! 9 ~

9

O

tD

,r~

r-.4 r-4

c~

~4 r--I

1534

MATHEMATICAL MODELING

,•

Z( CM )

Q(L)

=

3.219E-03 6.438E-02 3.219E-03

Q(H) -:

R( Cr0 )

12

12

\ E,".>.,\ I0

10

i

2 (a)

3

0

I

1

2 (b)

i

3

Fro. 5. Representative plots at point at which flame brush has entered main chamber (t = 4.55 msec, 0 = 344~ (a) velocity vector plot; (b) CH 4 mass fraction contour plot. ture equals or exceeds the fictitious value. All other ceils must actually exceed the CH4-air lower ignition temperature (923~ before chemical reaction commences. Figure 6 shows the CH 4 mass fraction contour plot for the slow b u r n at 352.57 ~ (4.79 reset). Comparison with Fig. 5b shows that the slow b u r n acts like a spark delay of slightly greater than 8 - 2 / 3 ~ (0.24 reset), with the two ignition procedures yielding nearly identical results. Thus, the ignition simulation technique is not of basic importance. Figure 7 shows the velocity vector a n d temperature contour plots at 4.69 msec (348.7~ Note that the piston is still traveling upwards. A distinct 2D flame front surrounds the prechamber jet and b u r n i n g continues into the main chamber. If b u r n i n g had not occurred in the main chamber, a strong vortex w o u l d have formed as previously reported, la However, the velocity field in the m a i n chamber is driven m a i n l y by the c o m b u s t i o n process and is not strongly affected by the prechamber. It should be noted that since the nozzle is now unchoked, both hot c o m b u s t i o n

products and lean charge flow through the nozzle and into the prechamber. A late-time (6.19 msec, 402.7 ~ snapshot of the temperature a n d NO xdistribution is shown in Figs. 8a and 8b. The flame front has reached the piston face, b u t has not yet progressed into the corner. The NO~ m a x i m u m resides in the prechamber. Earlier, the NO ~m a x i m u m contour was p u s h e d into the main chamber due to the jet flow field; however, it has n o w receded into the prechamber as a result of the mixing processes b e t w e e n the two chambers. We also observe that at times s u b s e q u e n t to 6.19 msec, the b u l k of the u n b u r n e d m e t h a n e resides in the corner of the prechamber. We find u n b u r n e d [ C H 4 ] approaching 200 p p m at that location. Due to the simplicity of the hydrocarbon oxidation reaction employed in our calculations, we believe that this result is qualitatively correct but its quantitative validity is questionable. Finally, in Fig. 9, we show the temporal evolution of total NO~ (in ppm) residing in the engine enclosure. We find that the NO x peaks out at about 8.8 ppm and is about a

COMBUSTION IN A CHARGE ENGINE

\M.f

0(L~ = 0(H)

1535

2.6q8E-03 5.297E-02 2.6qSE-03

= =

12

Z( CM ) I0

i

i 1

i

2

3

Fic. 6. CH 4 mass fraction contour plot for "slow burn" ease at t = 4.79 reset, 0 = 352.6~

factor of 4 higher than reported by us previously in a 1-D calculation. The 1-D calculation neglected m u l t i e o m p o n e n t species diffusion which may have a slight effect, b u t the m a i n difference in the two calculations is the mixing between the m a i n chamber charge and the preehamber charge after u n e h o k i n g of the nozzle. The presence of O 2 results in c o n t i n u e d c o m b u s t i o n of the rich preehamber charge and a higher temperature (by about 300~ in the preehamber. This result highlights the importance of chamber geometry and nozzle configuration. 4. D i s c u s s i o n and Conclusions We have presented a m u l t i d i m e n s i o n a l computer simulation of the compression and power strokes of a divided-chamber stratified charge engine. Our numerical procedure represents a considerable improvement over other techniques in that the arbitrary LagrangianE u l e r i a n coordinate system allows us to inelude arbitrary piston motion and variable, nonorthogonal, computational zones resulting in a more realistic simulation of the complex e n g i n e geometry.

We have observed that the vortex motion i n d u c e d in the piston face-cylinder wall corner for the 6000 RPM operating point leads to a complex flowfield whose dynamics are determined by interaction of the free stream flow, the wall b o u n d a r y layer and the pressure waves originating at the piston face. A high-resolution computational grid and an accurate turbulence model are necessary to resolve the boumdary layer flowfield with sufficient accuracy to determine the wall vortieity. Thus, while our computed compression-stroke flowfields await these improvements before quantitative predictions ean be made, we believe that the qualitative features are correct. In recent work, we reported significant vortex motion associated with the entrance of the preehamber jet into the high-expansion ratio m a i n chamber and its ultimate interaction with the piston. ]3 This p h e n o m e n o n is modified as a result of the propagation of a b u r n i n g jet into the main chamber and ignition of the lean mixture. The flowfield i n d u c e d by the propagating flame, coupled with the withdrawal of the piston, precludes the formation of a strong vortex. However, mixing between the two chambers still occurs with the NO x peak shift-

1536

MATHEMATICAL MODELING

lq

Q(L) = 5.950E+02 Q(H) = 2.199E+03 DQ = 8.~IE+OI

14 A( CM )

12

12

-.... -..,....I0

i/ . 11 0

I

2 (a)

,o

i

i 1

2

i 3

(b)

FIG. 7. Representative plots at t = 4.69 msec, 0 = 348-3/4~ (a) velocity vector plot; (b) temperature contour plot.

ing from p r e c h a m b e r to m a i n chamber a n d then back again. We believe that the shifting of the N O x peak b a c k into the p r e c h a m b e r is primarily a result of the adveetion of oxygen from the main c h a m b e r into the prechamber. We also find that, in spite of mixing b e t w e e n the chambers both before and after ignition, a large pocket of u n b u r n e d fuel remains in the p r e e h a m b e r corner at the end of the calculation. It appears that the geometry of the two chambers and their interconnecting nozzle p l a y a large role in engine performance. This observation agrees with the experimental observations reported p r e v i o u s l y b y Sakai, et al., 34 Date, et al., a5 a n d Yagi, et al. a6 Also, as a result of large radial gradients observed in the f u e l / a i r mixture ratio in the prechamber, spark p l u g location a n d ignition t i m i n g m a y be critical. The m i x i n g between lean m a i n chamber charge and rich p r e c h a m b e r charge a n d its effect on ignition timing has b e e n reported in the experimental study b y Purins. a7 These p h e n o m e n a c o u l d result in large cycleby-cycle variations in both performance a n d emissions.

Extensive parametric sensitivity calculations for variations in engine geometry, transport parameters a n d kinetic rate constants have not been p e r f o r m e d to date. However, we believe that each of the conclusions reached above is primarily a result of macroscopic f l u i d motion and will be most strongly affected b y engine geometry a n d o n l y marginally sensitive to transport parameters w h e n they are varied over realistic ranges. At the present time, we retain simplified combustion chemistry and turbulence formulations. An improved, multistep methaneair kinetics model a n d a three-equation turbulence model are currently under d e v e l o p m e n t and will represent the next improvements in the level of sophistication. These developments, in conjunction w i t h experimental information, should p r o v i d e for quantitatively accurate simulations with which parametric sensitivity studies can be performed. In the interim, however, some c o m m e n t should b e made regarding the interpretation of our results with respect to m o d e l assumptions, viz. spatially and temporally constant diffusivities a n d a single global reaction rate expression with

COMBUSTION IN A CHARGE ENGINE

15

1537

,5

\DO =

Z( CR )

l

13

11

=

3.95qE+0( 7.051E+0: 3.503E+00

13 R(CM)

11

9

2

0 fa)

2 (b)

FIG. 8. Representative plots near completion of main burning (t = 6.19 msec, 0 = 402-3/4~ (a) temperature contour plots; (b) NO concentration (in PPM) contour plot.

simplified Zeldovich NO~ kinetics. Also of some concern is the effect of numerical diffusion on the computed results. It has been recognized that in the numerical solution of the Navier-Stokes equations in E u l e r i a n coordinate systems, finite difference approximations to the advective terms all have the effect of introducing an artificial diffusion term into the solution. 1 For differencing with first order spatial accuracy, there exists a second order derivative with diffusion coefficient of order u A x / 2 . Thus, a cell Reynolds number, Re c - u A x / v t can be defined which is a measure of the ratio of "artificial diffusion" to the real, t u r b u l e n t diffusivity. Under ideal circumstances and when a t u r b u l e n c e model which yields accurate t u r b u l e n t transport coefficients is employed, it is highly desirable to m a i n t a i n Re c of order unity. I n our calculations, with the zoning and diffusivities employed, we estimate Rec to be between order one and ten. This did not lead to oscillations in the velocity or t h e r m o d y n a m i c fields. However, it is clear that numerical diffusion effects have led to deviations of the total diffusivity from the n o m i n a l value of 50

c m Z / s e c assumed for v t. This effect will be most p r o n o u n c e d in regions of high velocity and large computational zones, i.e., in the m a i n chamber around 90 ~ into the compression stroke a n d in the latter parts of the power stroke, as well as in the nozzle (high velocity), It is interesting to note, however, that the numerical diffusion tends to act as a "turbulence model," generating augmented diffusivities in regions of high shear. There are no currently available differencing schemes in E u l e r i a n coordinate systems which have been shown to eliminate numerical diffusion altogether, but reducing the zone size will tend to m i n i m i z e the effect. One could employ a purely Lagrangian formulation (where the advective terms do not appear in the equations), but this is not realistic for two-dimensional calculations with high shear rates. At the present time, since we have not employed a realistic t u r b u l e n c e model and are uncertain of the diffusivities within an order of magnitude, the effect of the enhanced diffusivity on the calculated results is difficult to assess rigorously. However, two points should be mentioned with respect to our findings. Our computed

1536

MATHEMATICAL MODELING

2

%

2D

/

/

/

/

/

/ /

/ !

/

!

3.0

!

!

/

I I I I

2.0

/

/

I I I I I

1.0 .9

.7

/ /

I

/

! !

.4

/

/

4.0

4.5

l 5.0

/

/

I

/

!

/

/

/

/ / /

/

!

/ .2

/

/

I !

.3

/

/

/

.5

f

/

/

.•

/

/ / / / /

i

I

5.5 6.0 t (~SEC)

I ~.5

7.0

7.5

Fro. 9. Total Engine NO vs. Time

b u r n i n g times are about 35 ~ crank angle, w h i c h is quite reasonable. Secondly, the real Reynolds n u m b e r of the flow is large due to the presence of the nozzle. Therefore, inertial effects tend to d o m i n a t e over viscous effects anyhow, thus r e d u c i n g the sensitivity of the computed results to variations in v c This conclusion may not hold for engines without prechambers. With respect to the assumption of a single constant value for vt, we would expect that this assumption leads to an underestimate of t u r b u l e n t transport near walls, at the throat, and in the vicinity of the flame (due to flame i n d u c e d turbulence). O n the other hand, i n regions of low velocity, such as in the early part of the compression stroke, we clearly

overestimate the t u r b u l e n t transport. One other "model a s s u m p t i o n " and its affect on the computed results requires discussion, i.e., the chemical model. Clearly, a large variation in the preexponential factor or in the activation energy can lead to significant differences in the computed flowfield and NO~ levels. We demonstrated previously that a (non-physical) increase of the preexponential factor of Eq. (9) resulted in such rapid energy release rates that a detonation wave rather than a deflagration wave resulted. 1~ This observation clearly illustrates the strong c o u p l i n g between the chemistry and the hydrodynamics, and the parametric sensitivity to energy release rates w h i c h result from a (clearly non-physical) large increase in the rate con-

COMBUSTION IN A CHARGE ENGINE stant. Bather than varying the rate constant employed in Eq. (9), we believe that it is more realistic to employ a two-stage combustion model comprised of a m e t h a n e b u r n stage followed by oxidation of CO, c.f., Dryer and Glassman. a8 Utilization of the single ratel i m i t i n g global step should result in a slight overestimate of the C U 4 b u r n i n g rate, thus leading to an overestimate in NO~ levels. However, neglect of s u p e r - e q u i l i b r i u m Oatom concentration in the NO~ model should yield an underestimate. Quantitatively accurate predictions of the pollutant level thus await use of an improved chemical formulation. The results of the numerical simulation indicate that many essential features of the reciprocating, spark ignition engine can be computed and are in qualitative agreement with available experimental observations. As such, multi-dimensional numerical simulation appears to be a new, useful and relatively inexpensive (~$600 of C D C 7600 time for a calculation) tool to be used in conjunction with experimental and design studies associated with IC engine research a n d development. With the development of improved turbulence, chemistry and finite difference formulations, m e a n i n g f u l parametric sensitivity studies can be performed. Finally, we must m e n t i o n the critical need for experimentation on IC engine combustion dynamics. Currently no adequate data on incylinder combustion variables exists with which to compare our results. Recent research performed in the automotive industry research laboratories aa-ar,ag,4~ can provide some information on both divided-chamber and openchamber SCE performance, b u t measurements on turbulence, species concentrations, etc. are lacking. A detailed discussion of research needs in the automotive c o m b u s t i o n area was presented by B o n i )

Acknowledgment This work was sponsored by NSF/RANN under Grant A'ER75-08441 under the cognizance of Dr. D. Senich.

REFERENCES 1. BONI, A. A.: Workshop on the Numerical Simu-

lation of Combustion for Application to Spark and Compression Ignition Engines, NSF/RANN Report (October 1975). 2. BaACCO, F. V.: Preface to a Special Issue on Stratified Charge Engines in Combustion

1539

Science and Technology 8, 1-3 (1973). 3. LAvoIE, G. A., HEYWOOD,J. B. AND KECK, J. C.: Combustion Science and Technology 1, 313 (1970). 4. BLUMBERG, P. N.: Comhustion Science and Technology 8, 5-24 (1973). 5. TABACZYNSKI,R. J.: "A Review of Pollution Control in Internal Combustion Engines." Paper presented at the Workshop on the Nmnerical Simulation of Combustion for Application to Spark and Compression Ig'nition Engines, NSF/BANN Report, 3-79-3-104 (October 1975). 6. SImGNANO,W. A.: "One-Dimensional Analysis of Combustion in a Spark-Ignition Engine." SAE Paper 719010 presented at the 1971 Intersociety Energy' Conversion Conference, Boston, Mass. (1971). 7. ROSENTWEIG-BELL~N, J. AND SlHIGNANO, W. A.: Combustion Science and Technology 8, 51-68

(1973). 8. Beacco, F. V.: Coinbnstion Science and Technology 8, 69-84 (1973). 9. BrtAcco, F. V.: "Introducing a New Generation of More Informative Combustion Models.'" SAE Paper No. 741174 presented at the International Stratified Charge Engine Conference, Troy, Michigan (1974). 10. BONI, A. A., CH.ad~MAN,M. AND SGIINEYER, G. P.: "A One-Dimensional Variable Area Computer Simulation of Combustion in a DividedChamber Stratified Charge Engine." ASME Paper No. 75-WA/DGP-1 (1975), 11. GUPTA, H. C., BRACCO,F. V., AND WESTBR(X)K, C. K.: "Matheinatical Modeling of Two-Phase, Unsteady, Two-Dimensional Flows." Presented at the 5th International Colloquium on Gasdynamics of Explosions and Reactive Systems, Bonrges, France. To appear in Acta Astronautica. In press (1976). 12. BKacco, F. V., GUPTA,H. C., KRISIINAMUR'I'HY,L., SANTAVICCA,D. A., STEINBEltGER, R. L., .aNDW*~ SHAW,V.: "Two-Phase, Two-Dimensional, Unsteady Combustion in Internal Combustion Engines; Preliminary Theoretical-Experimental Results." SAE Paper No. 760114 (Feb. 1976). 13. BONI, A. A., CHAPMAN',M. AND SCHNEYER, G. P.: "Computer Simulation of Combustion Processes in a Stratified Charge Engine." Acta Astronantica 3, 293-307 (1976), 14. HINT, C. W., AMSD~N,A. A. AND COOK, J. L.: J. Comp. Phys. 14, 227-253 (1974). 15. GORDON, S. AND McBUIDE, B. J.: "Computer Program for the Calculation of Complex Chemical Equilibrium Compositions, Rocket Performance, Incident and Reflected Shocks and Chapman-Jouguet Detonations." NASA Report SP273 (1971). 16. STU~, D. R. ANDPROPr~ET,H.: "JANAF Thermochemical Tables," 2nd Ed., National Bureau of

1540

17.

18.

19.

20.

21. 22. 23. 24.

25. 26.

27.

MATHEMATICAL M O D E L I N G

Standards Publication NSRDS-NBS37, U.S. Govt. Printing Office, Washington, D.C. (1971). BOWMAN, C. T.: Emissions from Continuous Combustion S~tstems, W. Cornelius and W. G. Agnew (Eds.), p. 99, Plenum Press (1972). MELLOR, A. M.: Emissions from Continuous Combustion Systems, W. Cornelius and W. G. Agnew (Eds.), p. 23, Plenum Press (1972). WITZE, P. O.: "Hot-Wire Turbulence Measurements in a Motored Internal Combustion Engine." SAND75-8641, Sandia Laboratories Livermore (April 1975). WITZE, P. O.: "Hot-Wire Measurements of the Turbulence Structure in a Motored Spark-Ignition Engine." AIAA Paper No. 76-37, presented at the 14th AIAA Aerospace Sciences Conference, Washington, D.C. (January 1976). WITZE, P. O., Private Communication (1976). SEMENOV,E. S.: Instruments and Experimental Techniques 1, 102 (1958). SEMENOV, E. S. AND SOKOLIK, A. S.: Izv. Akad. Nauk SSSR, Otd. Tekhn. Nauk. 8, 130 (1958). SEMENOV,E. S.: "Studies of Turbulent Gas Flow in Piston Engines." NASA Tech. Translation F-97 (1963). IVANOV, V. N.: Izv. Vysshikh Ucheb. Zaved Mashin 3, 91 (1964). MOLCHANOV, K. K. AND MOSKOVSKOGO,T.: Auto. Inst. Autotransizdat, Moscow No. 17, p. 85 (English Translation, Shell Trans. No. 1019). HORVAT1N,M. ANDHUSSMAN,A. W.: Measurement

of Air Movements in Internal Combustion Engine Cylinders. DISA Information No. 8, 13-22 (July 1969). 28. ANDREWS, G. E., BRADLEY, D. AND LWAKABAMBA,

S. B.: Combustion and Flame 24, 285-304

(1975). 29. LAUNDER, B. E. AND SPALDING, D. B.; Lectures in Mathematical Models of Turbulence, Aeadeinic Press, New York (1975). 30. TABACZYNSKI,R. J., HOULT, D. P., AND KECK, J. C.: J. Fluid Mech. 42, 249-255 (1970). 31. DANESHYAR,H., FULLER, D. E., AND DECKKER, B. E. L.: Int. J. Mech. Sci. 15, 381-390 (1973). 32. PATTERSON,D. J.: "Cylinder Pressure Variations. A Fundamental Combustion Problem." SAE Transactions, 75, Section 1, 621-632 (1967). 33. WINSOR, R. E. AND PATTERSON, D. J.: "Mixture Turbulence--A Key to Cyclic Combustion Variation." SAE Transactions, 82, Section 1, 368-383 (1973). 34. SAKAI,Y., KUNII, K., TSUTSUMI, S. AND NAKAGAWA, Y., "Combustion Characteristics of Torch Ignited Engines," SAE Paper 741167 (1974). 35. DATE, T., YAGI, S., ISHIZUYA, A., AND FUJII, I., Research and Development of the Honda CVCC Engine," SAE Paper 740605 (1974). 36. YAGI, S,, FuJII, I., WATANABE, M. AND NARASAHA, S., "On the Emission-Combustion Temperature Relationship in the CVCC Engine," SAE Paper 760109 (1976). 37. PURINS, E. A., "Pre-chamber Stratified Charge Engine Studies," SAE Paper 741159 (1974). 38. DRYER, F. L., AND GLASSMAN, I.: Fourteenth Symposium (International) on Combustion, The Combustion Institute, 987-1003 (1973). 39. BRANDSETTER,W. R., AND DECKER, G.: Combustion and Flame, 25, 15-23 (1975). 40. LAVO1E,G. A., ANDBLUMRERG,P. N.: Combustion Science and Technology 8, 25-38 (1973).

COMMENTS Alan Mitcheson, University of Leeds, England. Firstly, let me congratulate the authors on a well presented and comprehensive piece of work. However, the model only relates to one quarter of a 4 stroke cycle. As the remaining 3 strokes will form the initial conditions for the stroke being modeled, would the authors expect such difference in their solutions if such effects were accounted for?

Authors" Reply. We have simulated both the compression and the power strokes of the cycle, but have not calculated the intake and exhaust strokes explicitly. To do this would require extensive modeling of the intake and exhaust system components, which must ultimately be done. While the calculations performed to date have been "initialized" with

an assumed level of turbulence and a fuel/air mixture ratio distribution in a quiescent state, perturbations from these conditions, including a velocity distribution, could be assumed on an ad hoc basis to investigate the effect of such variations on the engine performance. We w o u l d expect that a high level of turbulence caused by the intake stroke would affect the ignition and initial flame propagation results. However, experimental results obtained by Witze (cited in the paper) indicate that the turbulence generated by the inlet port has pretty much been dissipated by the beginning of the compression cycle. More extensive experimental data over a wide range of parameters would be useful here. Similarly, volumetric efficiency effects will change the state somewhat, but probably not significantly. We expect that the most significant effect will result from the

C O M B U S T I O N IN A C H A R G E E N G I N E residual gases (burned and unburned) that remain in the chamber rather than being swept out as in an "ideal engine." However, before one attempts to model the exhaust cycle, we believe that an improved cylinder wall boundary layer/hydrocarbon-quench layer formulation should be incorporated into the model. Clearly, the full cycle must ultimately be simulated before a quantitative predictive capability can be claimed.

William.[. McLean, Comell University, USA. The very important contribution from the inlet port flow to the initial turbulence properties of the cylinder gases has not been included in this model. Also, the prechamber flow usually induces swirl in the main chamber gases. Has this effect been included?

Authors" Reply. This question is related to tlaat posed by Mitcheson and has been partially answered there. As we indicated, neither effect has been included, but for differing reasons. Witze's experiments, cited in the paper, have indicated that the turbulence generated by the inlet port has almost dissipated by the beginning of the compression stroke. Thus, it was not included explicitly in our calculations, although it could be on an ad hoc basis. Swirl is induced in the H o n d a CVCC engine by the non-axisymmetric geometric offset of the prechamber. Our axisymmetric model cannot accommodate that effect since it is 3-D. While extension of the model to 3-D is possible (and not difficult in principal) we have not done so because we believe

1541

it is premature. Furthermore, computer storage restrictions would limit resolution a n d sub-grid scale model detail. Work is in progress to include a swirl velocity (re) in a 2-D context, i.e., without any azimuthal (0) variation in parameters, thus retaining the axisymmetrie approximation. With this capability, an ad hoc introduction of an initial swirl could be included in the simulation. We hope to report on such calculations during the next year.

D. B. Spalding, Imperial College, England. Comparisons with experimental data are certainly desirable; but these need not be, indeed at first .should not be, those having the full complexity of the engine. The grid employed in your work is suitable for predicting steady flows in venturi nozzles, for which experimental data are already available. Have you made comparisons with these or similarly simple situations? Authors" Repl~j. We have not used the computer model to simulate steady venturi or nozzle flows and agree that this is a good idea. Simple spherically symmetric flows (burning and non-burning), shock tube flows and channel flows have been simulated with success. On the basis of our experience with this and other similar computer codes, we believe that the fluid flow predictions are valid. In our opinion, however, the more important question has to do with the accuracy of simulating reactive flows rather than non-reactive ones. Unfortunately, such data are sparse and the need for such data, particular for simple flow configurations, is great.