A simulation model including intake and exhaust systems for a single cylinder four-stroke cycle spark ignition engine

Int. J. mech. 8¢L Pergamon Press. 1975. Vol. 17, pp. 97-124. Printed in Great Britain

A SIMULATION MODEL I N C L U D I N G I N T A K E A N D E X H A U S T SYSTEMS FOR A SINGLE C Y L I N D E R F O U R - S T R O K E CYCLE S P A R K IGNITION E N G I N E R . S. BENSON,* W . J . D . ANNAND~" a n d P . C. BARUAtt* *Mechanical Engineering Department, University of Manchester Institute of Science and Technology, Sackville Street, Manchester; and tPollution Research Unit, Manchester University, Manchester, England

(Received 15 March 1974, and in revised form 19 August 1974) S u m m a r y - - A comprehensive simulation model is presented for a spark ignition engine including intake and exhaust systems. The power cycle simulation requires only one empirical factor to correct for turbulent speed of the flame front in order to complete the cycle calculation including NO. The exhaust pipe gas dynamics include chemical reactions along p a t h lines. Calculations are presented which compare well with experimental results. The model predictions compare favourably with previous work. NOTATION speed of sound aq Annand constant for convective heat transfer aA speed of sound after isentropie change of state to reference pressure a r e f reference speed of sound A non-dimensional speed of sound (a/aref) Aa non-dimensional a~(a~/aref) C~ specific heat a t constant pressure C~ specific heat at constant vohune D diameter f friction factor -rw/½pu2 F area g specific Gibbs free energy function k specific heat ratio (C~/C~) K thermal conductivity of fluid K, rate constant for equation (i) K~ equilibrium constant a

m

mass

n P

rotational speed pressure reference pressure rate of heat transfer per unit time per unit mass of fluid total heat flux pressure ratio :P/Pref o r radius gas constant universal gas constant Reynolds number one-way equilibrium rate for reaction (i) time temperature for i t h zone gas temperature particle velocity or specific internal energy non-dimensional particle velocity (u/aref) or internal energy 97

Pref

q

Q T

R Rmol

Re R~ t t~ T u

U

98

R . S. BENSON, W . J . D , AI~,kbID a n d P . C. B~-RU),H Ui U(z,j)

V lF X

X

X~ Xref

Z Zeoeff

Ae dZ P 'Tw

¢

internal energy coefficient for ith term internal energy coefficient for ith term and jth species volume or velocity work distance non-dimensional length (x/Xref) molar fraction of ith species reference length non-dimensional time [aret(t/Xref) ] viscosity coefficient crank angle Riemann variable A + [(k - 1)/2] U Riemann variable A - [(k - 1)/2] U stoichiometric coefficient increment of a increment of Z density wall shear stress equivalence ratio viscosity

Su2~xes

stagnation condition burned product or backward cylinder C equilibrium or exhaust e / forward i inlet in incoming m unburnt mixture n previous time step n + l current time step outgoing out P product "tO wall 0 b

INTRODUCTION SIMULATIO~ models can be of great assistance to the engine designers if t h e y give a good r e p r e s e n t a t i o n of the engine system. T h e i r economic value is in the r e d u c t i o n in time and costs for the d e v e l o p m e n t o f new engines a n d their technical value is in t h e identification of areas which require specific a t t e n t i o n as the design s t u d y evolves. Simulation models for complete engine systems including i n t a k e a n d e x h a u s t systems have been successfully developed for diesel engines. 1, ~ On the other hand, models for spark ignition engines h a v e been confined to the power cycle, in particular with the identification o f the p a r a m e t e r s which influence the e x h a u s t emissions, a I t is well k n o w n t h a t t h e q u a n t i t y o f residual gases in the cylinder a t the c o m m e n c e m e n t o f t h e compression stroke influences the nitric oxide emissions, a n d b y artificially increasing the residuals t h r o u g h e x h a u s t gas recirculation the nitric oxide emissions can be reduced. 4 Qualitative, and to a certain e x t e n t q u a n t i t a t i v e , work on h y d r o c a r b o n emissions indicates t h a t the u n b u r n e d h y d r o c a r b o n s leave the cylinder in two peaks, one a t the beginning o f the e x h a u s t stroke during blowdown and the other at the end of the e x h a u s t stroke.a 5 A b o u t

A model for spark ignition engine with intake and exhaust systems

99

half the total mass emissions comes from each peak. At present, the model simulations for both nitric oxide and hydrocarbon emissions neglect the interaction of the gas dynamics in the intake and exhaust system with the cylinder conditions. Since it would appear that control of exhaust emissions c a n b e a t t e m p t e d t h r o u g h t h e d e s i g n o f t h e i n t a k e a n d e x h a u s t s y s t e m , i t is c l e a r t h a t a s u i t a b l e m o d e l t o s i m u l a t e t h e s e s y s t e m s is d e s i r a b l e . I n t h i s p a p e r a m o d e l is p r e s e n t e d f o r a s i n g l e c y l i n d e r f o u r - s t r o k e s p a r k i g n i t i o n e n g i n e . I n a s e c o n d p a p e r t h e m o d e l is a p p l i e d t o a c r a n k c a s e s c a v e n g e d t w o - s t r o k e e n g i n e , e T h e p o w e r c y c l e w i l l b e first d e s c r i b e d f o l l o w e d b y t h e g a s e x c h a n g e p r o c e s s . N i t r i c o x i d e e m i s s i o n s o n l y will b e e x a m i n e d . T h e m o d e l w i l l b e compared with experimental results. POWER

CYCLE

MODEL

A number of power cycle models have been suggested. That of Lavoie et a l Y is generally recognized as giving the closest approximation to the combustion process in a spark ignition engine. However, the model depends either on a known pressure-time diagram or a semi-empirical burning law which would also include the ignition delay. Lavoie et a l 3 considered a temperature gradient in the cylinder related to the burning rate. I n this paper we propose to use the chemical reaction equations suggested b y Lavoie b u t the cylinder volume is restricted to two zones, a burnt and an unburnt zone. The two zones are separated b y the flame front. Furthermore, we propose to use a simple method for predicting the ignition time for the commencement of combustion. W i t h these two modifications to Lavoie's method we avoid the empirical burning law and simplify the computational procedure. However, the new method introduces a new variable, namely, the flame propagation speed. I t is to be hoped t h a t introduction of this variable m a y enable a more rigorous prediction of the combustion process t h a n a burning law. Since the rates of the energy producing reactions in the burned gas zone are so fast t h a t the burned gases are close to thermodynamic equilibrium, we assume t h a t the pressure and temperature in the burned gas can be computed from equilibrium thermodynamics. The nitric oxide concentrations, only, are computed from the rate equations. Other species concentrations are calculated b y numerical solution of the equilibrium equations based on a modification of Vickland's method. 8 The volume of the burnt gas zone during the combustion process is determined b y methods reported b y Annand. 21 Contrary to Lavoie et al. 7 the analysis includes heat transfer during the combustion process. Since equilibrium calculations are iterative, it is necessary to make close estimates of the product temperature if rapid convergence is to be achieved. Some of the techniques will be described. These are based on numerous calculations covering normal engine conditions. Before discussing the cycle calculation, we will briefly review the model for the mechanism of flame propagation and ignition. We will then review the ecluilibritnn thermodynamics and the rate kinetics. MECHANISM

OF

FLAME

PROPAGATION

Despite the excellent qualitative work, particularly from combustion photography, there still remains uncertainty about the fundamental details of the structure and behaviour of the combustion zone in a spark ignition engine. One area in which quantitativo d a t a cannot yet be reliably predicted is the mechanism of the flame propagation. The theory of laminar flame propagation through a homogeneous air fuel mixture is well established. The rate of propagation depends p a r t l y on the physical properties (specific heat, thermal conductivity, mass diffusivity) of the substance making up the unburned mixture and products and p a r t l y on the chemical nature of the reaction and the relation between the controlling reaction rates and temperature. Spalding et al. 1° have shown t h a t the propagation speed can be calculated if the reaction can be fully described. However, under turbulent conditions the flame front is modified and the propagation

100

1%. S. BENSOlV, W. J. D. ANNAND and P. C. BARUAH

speed increases with increasing turbulence. Until recently it was held t h a t the increase of flame speed was due to the wrinkling of the flame front (which increases its area) n in combination with enhancement of heat and mass transfer rates, and the turbulent propagation speed ut would be expressed as a simple function of the laminar propagation speed us and the r.m.s, turbulent velocity u'. This has been found to be true only for low intensity of turbulence. When u" becomes comparable with us it appears that the flame front is broken up and replaced by a much thicker zone containing small volumes of unburned, burning and burned material stirred together by turbulence. Sokolik et al. 12 and Chomiak is have presented some experimental evidence for this. Howe and Shipman 14 have proposed an elaborate statistical model characterizing the burning zone as spherical parcels of unburned gas dispersed through burned gas with burning at periphery of the parcels. At the moment, there is no precise relation between turbulent and laminar flame speeds for turbulence of high intensity. In their computer simulation work, Phillips and Orman 15 reviewed a number of theories for laminar flame propagation and their final choice was for the thermal theory of Mallard and Le Chatelier. 15 Bailey 17 a tt e m p t ed to use the Mallard and Le Chatelier expression for his work on combustion limitations of gaseous fuels in reciprocating engines. However, some difficulty was encountered relating to the adjustment of temperatures. This was attributed to the exclusion of a pressure term in the expression. KuehP s used measurements made on a flat burner to develop an empirical expression for propane-air mixture which correlated experimental flame speed data over a wide range of pressures and temperatures. This was successfully used by Bailey. 1~ Kuehl's expression for laminar flame speed us in a propane air mixture is ---- [ us

1"087 x 105

] p-O-OgS76, (1)

[((104/Tb)÷(900/Tu))4.gasj

where T,, T 5 are the temperatures of unburned and burned mixtures, K ; us is the laminar flame speed, cm/see ; and p is the pressure in. Hg. I n the present model Kuehl's expression for laminar flame speed is used as the basis of the calculation. To allow for turbulence a multiplying factor is used. This factor will be some function of the turbulence level in the cylinder and the engine speed. The turbulent flame speed is then given by u, = ff x u.

(2)

where f f is a flame factor determined by one of two methods. The best approximation of the value of f f can be obtained by carrying out a number of calculations with different values of f f and comparing the computed pressure-time diagrams, and the experimental pressure-time diagram. When the peak pressures in both the diagrams are matched the value of f f is obtained. The other approximate method is followed in the case when no experimental diagram is at hand. In this f f is adjusted until the computed pressure diagram gives a symmetrical burning time about T.D.C. Having obtained the flame speed it is assumed t h a t the flame propagates spherically from the spark plug. IGNITION

AND INITIATION OF TWO ZONES COMBUSTION SPACE

IN

In the model we shall use, ignition of the fuel is said to take place when a finite volume of the burnt mixture exceeds 0.001 times the total cylinder volume ~ (swept volume plus clearance volume). At the nominal ignition time the mixture (at state p,~, tin) is assumed to burn adiabatically at constant volume to the product temperature, tp. Fo r this purpose an internal energy balance is made for unit mass of the mixture. The f i r s t e s t i m a t e of the product temperature t~ is based on the following expressions: t~=tm+2500x~xf~

for~
and tr

= t m+ 2500 x ¢ ×f, -- 700(¢ -- 1.0)f,

for ¢ > 1.0,

where f~ is the fraction fresh mixture and ¢ is the equivalence ratio.

A model for spark ignition engine with intake and exhaust systems

101

An energy balance is made and t~ adjusted until the specific internal energy of the products = t h e specific internal energy of the reactants, and t, determined. The flame speed (us) is calculated from equations (1) and (2) where T~ = t v and T~ = t~. The volume of the burnt mixture V~ is given by

v, = ~ s , where r = ut(Aa/360n) is the radius of the flame front from the spark plug, Aa is the angle increment and n is the engine speed, rev/sec. F o r initiation of combustion V, is set equal to 0.001Ve and the angle (Aa)delar calculated from the above expression. I f the time from the nominal spark angle is less t h a n (Aa)delar the mixture is considered not to have been burnt and the compression process carries on for the next time increment da. The process continues for as m a n y time intervals as necessary until the total angle from the nominal spark timing is greater than (Aa)delar when combustion is said to have commenced. After the combustion of the small nucleus of fuel-air mixture the combustion chamber is subdivided into two zones; a burnt zone, suffix p, and an unburnt zone, snffix m. The process is initiated in three steps. Referring to Fig. 1, we show the three steps, A, B and C. Hea? loss qm

rap, fp", p','vp"

rnrn I = m

mm

tm

tin' p' v2

I p I v I

I

/

root 2

mm 2

--

tin"

(b) Initial s t a t e

mp, fp2,vp2

Final volume

vm v2 Transform at fixed volume

,

(c)

tm 2 vm 2 v2, p2 Equilibration of pressure

FIG. 1. Basis of the first combustion step calculation.

Step The temperature t~ and pressure p" are given by

p,

p~ V~ r

= v--~

•

(4)

Step B. Appearance of flame nucleus During this process t~ and p ' are kept unchanged. Using ~ as the product temperature for adiabatic constant volume combustion from t~ we can write P" = (~R----~-i~] R ' ¢ t P'.

(5)

To obtain the mass of the products m~ at this stage, we assume a slight increase in density of the products compared with the density of the mixture:

m~ =

2rrr 8 m 10 -1° -3~ + 2 x V1,

(6)

where m is the total mass and r, the flame front radius. The second term in equation (6) is the slight increase in density which has been found from a large number of calculations to give a good starting point for the calculation of this stage of the process. Now 8 m,,, = m-m~ (7)

102

R . S . BENSON, W. J . D. ANNAND and P. C. B ~ V A H

and

so t h a t t h e t o t a l internal energy is t h e n

U = mu'~ = m~u'~+m~u'~.

(9)

Step C. Equilibration of pressure This step is assumed to be an adiabatic c o n s t a n t v o l u m e process. The total internal energy is U = m~, urn, + m~ u~s. (10) Balancing t h e t o t a l internal energy before a n d after step C m~, u~ + m~ u~ = m~,, urn, + m~ uv,

( 11 )

which, after r e a r r a n g e m e n t , becomes

mm, Cvmt

-1

= m~Cv t~ 1 -

.

(12)

As before

tm

[p2~¢km-1)/~m

t'-~"~= [P-;]

'

(13)

t~ = V ~ !

'

(14)

and

P~ p~

P2 P ' = P_~,R,, t~ p ' p~ p ' R~ t~ "

(15)

Combining equations (12), (14) and (15) and letting

t~n, ---- 8,

k,,,

t~

kv- 1

k,~----i

k--7- = ~

and

C vp t"~, m~ m~ C~, t"

_

A,

{R., t'] ~p-1)/kp \R~!

=

B

wc h a v e - 1 = A(1 - B S - ) ,

(16)

which can be solved for 8 b y an iterativo technique. THERMODYNAMICS

Equilibrium thermodynamics The following species are considered to be present in the p r o d u c t in t h e cylinder and in t h e exhaust gases. T h e y are referred to b y the n u m b e r in brackets a p p e a r i n g against t h e i r names.

(1) H20,

(2) I-I2,

(7) N,

(8) CO~, (9) CO,

(3) OH,

(4) H ,

(I0) 0 2,

(5) N2,

(6) NO,

(ll) O,

(12) AR.

A r g o n is s e p a r a t e d off f r o m atmospheric n i t r o g e n in order to d e t e r m i n e t h e 1~0 concentration. The composition is calculated in t e r m s of m o l a r fractions of these species d e n o t e d by X~. The fuel is w r i t t e n as Cca Hhy Oox, a n d one mole of t h e t o t a l p r o d u c t s arises from a moles of fuel plus (1/¢) times t h e stoichiometric q u a n t i t y of air (~ = equivalence ratio).

A model for spark ignition engine with intake and exhaust systems

103

The equilibrium distribution of these species can be fully described b y the following seven reactions : (1) ½H~ > H, (17) (2)

½0~

> O,

(18)

(3)

½N,

> N,

(19)

(4)

2H20

> 2 H 2 + 02,

(20)

(5)

H20

> OH + ½H.,,

(21)

(6)

CO~+H~

> H2O +CO,

(22)

(7)

H,0+½N,

> H,+NO.

(23)

These reactions are the same as those used by Vickland et al. 8 The equilibrium constant K~ for the stoichiometric reaction between the substances A, B, C, D, vaA+vbB . " vcC+v~D, can be expressed asz9 /X~o X~d~ K~ = l~¢v° " - - " ~va! Z" ~v0+vd--Vo--Vo, (24) where v is the stoiehiomotric coefficient, X is the molar fraction and p is the total pressure. Using this expression of equilibrium constants for reactions given in equations (17)-( 23 ), can be written respectively as (1)

K~, = 4(P) Xd/dX,,

(25)

(2)

K~, = 4(P) Xn/dXzo,

(26)

(3)

Km ----4(P) X , N X s ,

(27)

(4)

K~, = pXzo/b 2,

(28)

(5)

K~, = ~(p) Xa/(b ~X2),

(29)

(6)

K~, = bXg/X 8,

(30)

(7)

K , , = ~/(p) X J ( b dXs),

(3])

where

b = X d X ~. The equilibrium constants (K~) for these reactions are given by z9

[

In K , = t Z

-- Z

--

(32)

The specific Gibbs function is given b y the expression

g(T)

Rmo I T

= ag(1 - I n T ) - b a T - ~

ca

do ee T2--~ Ta--~ Td-kg,

(33)

where ag, bg, cg, dg, % and ]c~ are constants, the values of which can be obtained from ref. (19). We have 13 unknowns, the 12 species fractious and the total number of tools, so we need 6 more equations for the solution. One of these is (8)

ZX,

= 1.

(34)

The remainder are the atomic mass balances for argon, carbon, hydrogen, oxygen and nitrogen. The method of solution of equations (17)-(34) is outlined in Appendix I. BASIC RATE KINETICS The following basic reaction is considered:

A+B

Kf Kb

• C+D,

(35)

104

R . S . BENS01% W. J. D. A~NAND and P. C. BARUAH

where Kf and Kb represent forward and backward rate constants. Using (A)/(A)~ = a,

(B)/(B)~ = fl,

(C)/(C), = ~,

(D)/(D)e = 8,

where the suffix e represents equilibrium, for the considered reaction, then for the forward rate = K I ( A ) (B) = aflKt(A)e (B)~ and the reverse rate = K b ( C ) (D) = ySKb(C)e (D)~. A t equilibritun the forward rate is equal to the reverse rate, K~(A), (B), = Kb(C)~ (D), = R~. Therefore, the net rate = K I ( A ) ( B ) - K b ( C )

(D) = ( o ~ f l - ~ ) R~.

NO KINETICS I t is now recognized t h a t the formation of NO in an engine combustion chamber is a non.equilibrium process. I n an earlier work lgewhall s° assumed t h a t the initial conditions for the expansion calculation are given by chemical equilibrium distribution of species corresponding to a combination of temperature and pressure representative of the engine subsequent to flame propagation. The assumption is based on two facts, firstly, the pressure and temperature at the initial stage of expansion are at a m a x i m a and the reaction rates are at a maxima, secondly, owing to the kinematics of the engine mechanism the piston velocity and hence the expansion rate is very low at the beginning of the expansion stroke. I n this quantitative theoretical analysis rate expressions for 16 elementary chemical reactions were combined to give 13 coupled nonlinear differential equations representing 13 chemical species. These equations were numerically integrated along with a coupled energy equation starting at the initial point of expansion where the chemical species were determined with equilibrium calculations. The time variations of the rate-controlled concentrations of the species from their equilibrium concentrations were shown. I n the present paper no such assumption is used and the reaction kinetics of nitric oxide are applied from the start of combustion. The rate kinetic model for NO is based on the theory developed by Lavoie et al. 7 with some modifications. The governing equations for the mechanism of NO formation are (1)

N+NO

(2)

N + 02

(3)

N+OH

~

~

'. . . .

.

N ~ +O ,

K l f = 3.1xl01O×eC-lso/r~;

(36)

* N O + O,

K~t = 6 . 4 x 106 x T x e ~-31~/r~ ;

(37)

~ NO+H,

K3f = 4.2x 101°;

(38)

~ N2+OH,

K4! = 3.0x 1010xeC-S350/T~;

(39)

N2+O2,

Ksf = 3"2× 1012xeC-ls,9°°/~;

(40)

(4)

I-I+N30 ~

(5)

O+N$O

(6)

O+N30 ~

• NO+NO,

Kel = K51;

(41)

(7)

N30+M

• N2+O+M,

K~I = 1013×e~-3o,3oo/~.

(42)

~

.

'. . . . .

The above rate constants are given in m3/kg mole per sec. At the temperatures and densities in an internal combustion engine it is a reasonable first approximation to assume t h a t the volume of the reaction zone is negligible and t h a t the gas within a cylinder consists of a burned fraction at thermal equilibrium plus an unburned fraction frozen at original condition. We will consider the process occurring in the burned gas behind the reaction zone only. The rates of energy producing reactions in a flame are sufficiently fast so t h a t the burned gases are close to thermodynamic equilibrium. I t can be assumed, therefore, t h a t H and O H are at equilibrium concentrations, and t h a t O is in equilibrium with 03. The rate equations for NO, N and N20 can then be developed. These are developed in Appendix I I ; it is there shown t h a t the

A model for spark ignition engine with intake and exhaust systems

105

rate equation for NO is the most important and is given b y 1 d

(

Rx

Re

~ h~ ((NO) V) = 2 ( l - - d e~) \ l + % [ R x / ( R s + R a ) ] - t 1 + [ R d ( R -4-+ R 5 + 27) ]]

(43)

where ~, denotes (NO)/(NO),, R is the "one-way" equilibrium rate for the i t h reaction, V is the volume of b u r n t gas and the suffix e denotes the equilibrium condition. The kinetic rates for N and NsO formation are much faster, and m a y be assumed to equilibrate instantaneously, giving values for (N) and (NsO). Along with (NO) determined by integration of the rate equation, a total nitrogen balance can be formed to yield (Ns) a t each instant. We now turn to the cycle calculation.

CYCLE

The cycle of events between air valve closing and exhaust valve opening comprises: (1) Compression stroke. (2) Ignition and propagation of flame front. (3) Expansion: (a) two zone; (b) full products. ( 1 ) Compression stroke The compression process starts at the t r a p p e d condition. The state of the gas a t this point is derived b y using a perfect mixing model for fresh charge and residuals from the previous cycle. During compression any reaction is neglected. The first law of thermodynamics is dQ

_ dT

dV

(44)

The equation of state is pV = mRT.

(45)

A combination of equations (44) and (45) after rearrangement gives

~= da

[_ (, +R)

dV t RidQll v

P -S-Ja+

~O-d-~ Jl

(46)

and dT

_/1

d~ =

~' I~ ~-g+~

dV

l~__pa)

•

(47)

F o r the unburned mixture, dt,~ (V d V + l dp~ d-'~ ---- tm " ~ p 'd"o~ "

(48)

The heat transfer rate from gas to wall is given b y Annand's equation for convective heat transfer 9

Q

N =

(Re) b ( t , . - t,o).

(49)

The work is dW

dV

d~ = p ~-g"

(50)

As the compression process continues the variables are incremented b y using the following general expression dx xn+l = x . +~--~ Aa, (51) where x is any variable. The numerical method used for this purpose is the R u n g o - K u t t a method.

106

]~.

S. BENSON, W. J. D. ANNAlUD and P. C. BAR~=AH

(2) Ignition and propagation of flame front W e h a v e already considered the process for ascertaining the ignition of the f u e l - a i r m i x t u r e . Once ignition has t a k e n place we consider the flame to be p r o p a g a t e d spherically f r o m t h e spark plug. The volumes of t h e b u r n t and u n b u r n t regions are c o m p u t e d by the m e t h o d s described elsewhere. 21 F o r this purpose, a n o t e m u s t be kept n o t only of the enflamed v o l u m e against the flame radius b u t also the flame front area and the flame c o n t a c t area w i t h t h e cylinder wall and piston. W h e n the flame is p r o p a g a t e d and t h e zones set up by the m e t h o d s described earlier, the calculation proceeds in a step-by-step m a n n e r using R u n g e - K u t t a integration. (3a) Expansion with two zones The m a j o r assumptions are : (a) the original charge is homogeneous; (b) t h e pressure at a n y t i m e is uniform t h r o u g h o u t the cylinder; (c) t h e v o l u m e occupied b y the flame reaction zone is negligible; (d) the burned gas is at full t h e r m o d y n a m i c equilibrium e x c e p t for t h e nitrogen species ; (e) t h e u n b u r n e d gas is frozen at its original composition; (f) b o t h burned and u n b u r n e d gases h a v e u n i f o r m local specific heats; (g) there is no h e a t transfer between b u r n e d and u n b u r n e d zones. I t can be seen t h a t this m o d e l is e x t r e m e l y simplified b u t experience shows t h a t the assumptions are v e r y well justified. Spark plug

Total volume = V

Pressure = P

FIG. 2. Combustion zone shape. Considering Fig. 2, the t o t a l internal energy for the s y s t e m is V

=

(52)

m m U m q - m ~ u ~.

The first law of t h e r m o d y n a m i c s is dQ dU dW d---~ = -d-~a+--~-a, dU

dum

(53) dram

du~,

dm~

(54)

and dm m dc~

dm~, do~ "

Therefore, e q u a t i o n (52) becomes dm~ dt,n _ dt~ d V dQ ----= O. (u~ - um) - - ~ + mm C~,, - ~ + m~ (/t,~,"~-~+ P -dc~ da

(55)

A model for spark ignition engine with intake and exhaust systems

107

dV dV~ dV~ d'-~ = d~ +-d-~'

(56)

~Now

and differentiating the equation of state (pV = toRT),

dV

(V~

V~dm,

m,~Rmdt,,

a-~= [-~,-~':~ ~

p

m~R~dt,

d~ ~ p

a~

Vdp

(57)

p d~"

Applying the first law to a unit mass of m~ du~ I dQ~ dV~ d a = m,~ da 'P--~-a '

(58)

where Q~ is heat transfer to the unburned mixture assumed uniformly distributed throughout the mixture and V~ is the specific volume of the mixture. We can write equation (58) as C

dt~ ~=m.~

1 dQ,n p V ~ , [ 1 dt~ d~ \t~ da

1 ~_~p) p

(59)

and after using the usual gas state relationships, we have dt~

m,~ C, -~--~= V,~ ~ + ~ - ~

.

(60)

F r o m equations (55), (57) and (60) we can obtain separate equations for dp/da, dt,~/d~ and dt~/da: F r o m equation (60), dt~ da

~

V~ dp 1 dQ +- -m~ C , . d a m~ C r . d a

(61)

°

F r o m equation (57).

d'~ = m.R--'~

\ :p

"p "] ~

pC.,,, da p C . . da

p

"

F r o m equation (55),

dm~ ×~

[C~,

C~ R,~ \ dQ,~ dQ - -

+ \ C ~ - Rr Cr,.) d~

d~ "

(63)

The heat transfer terms are determined from the following relations. F o r burnt mixture Q

a~ K~(Re) °'~ (t~ - tw) =

D

'

(64)

where (Re)~

-

p, DVptston and /~,

K, = C,,~, 0.7 "

F o r unburnt mixture

_Q = a~ K~(Re)~7 (t~-tw) 2'

D

(65)

where (Re),~ --- pm DVp~ton and /~m

K ~ = ~G*"/zm 0.7

The piston speed Vplston and cylinder bore D are used as the parameters in the Reynolds number expressions.

108

R . S . BE~so~, W. J. D. A ~ A ~ D and P. C. B ~ u A ~

The pressure changes dp/da can be obtained directly from equation (63) since all the terms are known. The temperature changes dt~/da and dt~/da can then be directly evaluated from equations (61) and (62), respectively. This completes the calculation for the two zone expansion for the time a °CA. We then calculate the new variables x at the next time step by the R u n g e - K u t t a method using the general expression dx xn+l = xn + ~ ha where x is a n y variable, the suffix n + 1 is the new time step and Aa is the time step increment in crank angle.

Time step I n the R u n g e - K u t t a method the crank angle increment Aa is the parameter which decides the accuracy of solution and the computing time. I f Aa is iarge the results will be inaccurate and if it is too small, the computing time will be higher. To make a compromise, a test was developed to m a i n t a i n the right level of accuracy while adjusting the value of As. I n this method all the variables (except the variables heat transfer and work output) during the expansion stroke are tested to see if the increment of a n y of them during a time step (As) has decreased b y more t h a n 40 per cent of the variable itself. I f this is found to occur the time step is halved and the calculation repeated. I n the actual computing, with a n average value of Aa -- 0.1-0.5 °CA during the expansion stroke, reasonable results are obtained. I t is to be noted that during the two zone period the time step is divided into a n u m b e r of sub steps.

Termination of combustion As the combustion process continues the product volume increases at a rate determined by the flame speed and correspondingly mixture volume decreases b y the same amount. The completion of combustion is detected by the current mixture volume. Theoretically, combustion should terminate when V~ ~<0. I n the numerical method it is assumed to terminate at the beginning of the time step in which the current value of V~ is just negative. (3b) Expansion with .full product Once the combustion is complete the variables are organized to calculate for single zone only. The R u n g e - K u t t a method is used. Throughout the expansion calculation a check is kept to see if NO is frozen when the NO rate kinetic calculations are by-passed. This completes the power cycle. We now t u r n to the gas exchange process.

GAS EXCHANGE

PROCESS

We will first consider the gas dynamics in the pipe then examine the cylinder boundary conditions between pipe and cylinder. Finally, we examine the cylinder thermodynamics.

Pipe Benson et al. ~adeveloped methods for studying non-steady flows in gases with variable composition and specific heats. These methods gave comparable results with the normal method for gases of constant specific heats2 2 We propose, therefore, to use the normal method ss with modifications to allow for the calculation of the gas composition along the path lines, and for the variable specific heat due to temperature and composition. We will start b y reviewing the basic equations then proceed to the method of solution with variable gas composition. For one-dimensional non-steady flow with friction a n d heat transfer the basic equations z2 are Continuity:

~p egu ~p pu dF O. ~T+P~+u~+-~ -~ =

(66)

A model for spark ignition engine with intake and ~xhaust systems

109

Momentum: ~u ~u 1 ~p ~+u~+~+~,

~,

= o,

(67)

where F ' is wall friction defined b y

F"

4"fu~ u = D 9 lut

and f =

~'w

~"

The first law of thermodynamics (energy equation) is

qpFd~ =

T + ~ + y ) j d~,

(6S)

where q represents rate of heat transfer per unit time and per u n i t mass of fluid. Using continuity and m o m e n t u m equations and simplifying, the energy equation is obtained in the form ( k - l)p(q-i-uF') = -~ . . t - u - ~ - a

(-~+U~x).

(69)

The heat transfer rate q m a y either be through the wails or be the result of longitudinal heat conduction. I n the event of a chemical reaction with change of state these will be represented b y a heat addition equal in magnitude to the constant volume heat of reaction of the actual process. ~a The equations of continuity, energy and m o m e n t u m are solved b y the method of characteristics. For the continuity, energy and m o m e n t u m equations this gives the following characteristics : dx d--/ = u _+a (characteristics) (70) and dx d-~ = u

(path lines).

Combining the continuity, energy and m o m e n t u m equations with the characteristics we obtain dp du ( 4 ] u s u \ a2pud-F , 4 j ' p a u ~ u O, (79`) d--t+-Pa'~ - ( k - 1 ) p _ q ' I D 2 lu]) + "~ ~ x ~ -~ l u l which is the compatibility condition along dx/dt = u +_a. Combining the continuity, energy and m o m e n t u m equations with path lines, we obtain _ ~ _ a 2 . _ dp ~_(k_l) ( q-~. u--~ 4 . f u--~~ Fu-~ u ] ---- O, (73) which is the compabflity condition along dx/dt --. u. Defining the R i e m a n n variables 2~ as

= A+~-~U, k-1

~ = A---ya~d A ---~

,

u -- ~=--~

we have the characteristic solutions in the form: Direction condition : dX = U+A. dZ

(74)

110

R . S . BENSOI~, W. J. D. ANNAI~D and P. C. BARUAH Compatibility condition: dA

k-IAU dFdz AdAa 2 F d-X + Aa

k - 1 2fXrefu~ U 2" D [UI ( k - 1) 2 qXrof i

Direction condition : dX -~ = U-A.

(76)

Compatibility condition : kdfl =

2

1 A U d F d Z + A d A a + k - 1 2fXre I U2 U F (iX ~ 2 D [U]

×

-Y-

ar--~ef $ dZ.

(77)

Path line characteristic: dX

= U.

(78)

k - 1 A a [ q X r e f + 2fXre t Ua ) = - - ~ --~ A ~ aaref D [ I dZ.

(79)

Compatibility condition: dA a

A numerical method for solving these equations has been described by Benson et al. 2~, 2~ I t has been tested extensively on various engines with good success. 1, 2 The heat transfer term in the )t, fl characteristics is (d~)heat transfer -----

( k - 1)' qL~2f 1 dZ 2 aref A

and in the path line characteristic is (dAa)heat transfer

k- 1 q~f.~2d Z = ~ aref A "

A computational scheme using these terms for wall heat transfer in both subsonic and supersonic flows using Reynolds' analogy is given by Benson. 25 I n the case of chemical reaction the heat of reaction is included in q. q = C/heattransfer + qreactlon.

A detailed method for including heat of reaction in the method of characteristics for solving a problem of a flame tube has been given by Benson. 24 P a t h lines

The identity of a fluid element is preserved along a path line. I n the present calculation the path line composition is defined by the molar concentration of the 12 species listed earlier. Chemical reactions are considered aJong the p a t h lines. A t each time step the composition is calculated from the thermal equilibrium equations (17)-(34), with the NO, N20 and N kinetically controlled through the expressions (36)-(43). F r o m the chemical reaction data the heat of reaction can be evaluated. I n the calculations described in this paper no heat of reaction was included, but it will become important in thermal reactors. F r o m a knowledge of the chemical composition and the temperature the internal energy of the mixture and the gas constant R can be obtained whence the ratio of the specific heats can be determined for each path line. The following procedure is used. Initial values of the species are entered at time t = 0 from the release conditions of the first power cycle.

A model for spark ignition engine with intake and exhaust systems

111

Considering a path line, Fig. 3, the pressure and temperature are calculated from the expressions ~'A+18~~ki'k-1 , P =

t"£~-~J

P'°f'

(80)

2+ From p and T the equilibrium composition is calculated from equations (17) to (34) using the same method as in the cylinder calculations, and from the equilibrium values of (NO)e, (N),, (N~O), the kinetically controlled amounts of NO, N and NsO are evaluated as in the cylinder calculations from equations (36) to (43). I n practice, the equilibrium calculations and the rate calculations are organized as subroutines which can be entered either from the cylinder or the pipe calculations. We now have all the species, the pressure and temperature; we can calculate, therefore, the ratio of the specific heats, ]c, in the following manner.

Enter species from previous time step Z+AZ / .

.

AZ /Path

.

.

.

A~k

Interpolate species ot/I and 6'

P,'T' -[

.

-

/

lines

Boundary __

Z

pj " ~

A \ -r-

B

/ / I

I

I

I

Boundary

[NO] Along pofhhne is rote controlled

Fro. 3. Calculation of path line species. The specific internal energy is u = Ul+U ~ T + u a T 2 + u l T a + u a T l

(82)

and Cv --- ( ~ )

= u2 + 2ua T + 3u, T~ + 4us T '.

(83)

V

For a mixture of species C. = ( Z x n u ( ~ , . i ) + 2 ( Z x . u ( a , n ) ) T + 3 ( Z x . u ( a , . ) )

T ~ - b 4 ( Z x . u ( a , . ) ) T a,

where n is the number of species and u is the internal energy of the mixture. Molecular weight of mixture m = ~m~x~,

(84)

(85)

where m is the molecular weight, x is the species fraction, i is the number of species, C~ = C~ + R

and

k

=

c,,IC,,.

(86) (87)

The effects of variable gas composition and variable specific heat are identified by k in the path line and characteristic compatibility equations (75, (77) and (79). I n the boundary calculations we include the effect of variable composition in the following manner. The botmdary species are entered on the boundary path lines depending on the direction of flow, The general approach at this stage is to enter the species coming from the boundary on to the boundary path line if there is flow into the pipe. Owing to the nature of calculation, the values are entered from the previous time step. I f flow is leaving the pipe the species are interpolated between the inside and outside path lines (Fig. 3). We will consider three cases: (a) cylinder, (b) junction and (e) nozzle and open end.

112

R. S. BE~so~, W. J. D. ~NNAND and P. C. BARUAH

(a) Cylinder boundary I n this case when there is outflow from the cylinder into the pipe the pipe end boundary species are filled in with the cylinder species. I f there is inflow from the pipe into the cylinder the boundary species are obtained by the interpolation of the inside and outside path lines at the pipe end. (b) Junction boundary The junction boundary uses a constant pressure perfect mixing model. The principle applied here is that only the boundary species flowing into the junction will participate in mixing (Fig. 4). The resultant mixture species will be entered for the pipe ends where flow is moving away from the junction. The pipe ends, where flow is into the junction, will retain their own incoming species.

Retain~ incominspecies g ~

=.

spectes

~

- v e flow adjust outgoing

FIG. 4. Junction boundary with species. The flow direction and the mass flow parameters axe obtained using the conventional expressions. Using suffi~ nn to denote a pipe end at the junction, we can write. Velocity parameter:

u.. = ~.-Aou~..

(ss)

Mass flow parameter:

w..

r.. F...

(s9)

Considering positive flows only, represented b y W+n (Fig. 4), we write the average concentration for species i

(X,)~, =

X (X,).. W+. X W+~,

(90)

For pipe ends having negative flow (x,)..

= (x,)~,.

(91)

A similar procedure will be applied for NO, R and k. The equilibrium a n d rate equations are not included in the junction boundary.

(e) Nozzle and open end boundary I n this boundary when flow is going out of the pipe into the ambient condition the species at the pipe end are obtained b y the interpolation of the outside and inside p a t h lines existing at the pipe end. When there is inflow into the pipe the values from the previous time step are used, to fix the pipe end species at the current time step. This is a n approximation, in actual practice a slight amount of ambient species m a y enter into the pipe. However, it is expected that mass blown back will be small and this approximation will have negligible effect on the prediction of NO. For engines with excessive back flow the calculation will have to be modified to allow for the mixing of air and gas. To obtain the integrated mass flow for a n y gas X we use the expression X - Z- m(X*) ~ , where X~ is the mass/mass concentration and m is the total mass flow rate.

(92)

A model for spark ignition engine with intake and exhaust systems

113

C&der In this model we assume that the composition of the products of combustion are frozen at release (e.v.0) and remain so during the gas exchange process in the cylinder. When the gases are in the pipes the residence time will be longer than in the cylinder and the equilibrium and rate equations are used. The value of k in the cylinder will depend, therefore, on the temperature in the cylinder until the air valve opens. At this point we assume for a four-stroke engine that the incoming fuel-air mixture mixes with the products of combustion but there is no chemical reaction. We will now require to evaluate continuously the composition of the cylinder contents to determine k. For flow through the inlet and exhaust valves, we assume a quasi-steady theory. The Riemann variable hi, and the entropy level A, at a time step are known from the previous time step. They have, however, to be modified for inflow to the pipe for the current time step. An outline procedure for the boundary conditions for the valves follows. We will then describe the thermodynamic relations for the gas exchange process in the cylinder. Flow through a valve Outflow. The boundary conditions at a given time Z will be dependent not only on the flow through the valve but also on the conditions immediately downstream of the valve. The Riemann variable Xi, is modified by the entropy A,, of the gas entering the pipe. The following expression is used :

Ano = hn.+w4

~&,-4m)~

WV

where hi, Ois the corrected Riemann variable, hi, n is the uncorrected Riemann variable, A,, is the entropy of the gas leaving the cylinder, A, is the entropy of the gas in the pipe and (A/A,) is the value of (A/A,) either at the previous time step of the current time step (usually taken as (hi, C+bUt, ,)/2Aa0 from the previous iteration). The flow through the valve is based on the constant pressure model.rs* a7 The general forms of the boundary equations are, using the parameters, A,=

a, u=;9 a,,f’

e=; 0

and *(k - 1) U” c = (l-*(k-l)r_P)“~

(94)

Subsonic flow in valve @

= [& (lc,1/($“+ PC) +)]*‘(?

(95)

Sonic boundary (A

0) >$I~,

sonic flow.

(96)

Choked jlow in valve (%-J = ~(~)‘~+l)‘z’l-l)(l_~Uz)t.

(97)

Entropy A

) (&$1)‘2,

BC

(P$.-ll,:

(93)

Axlt (99) The computational procedure for a given value of k is to first estimate the entropy Let A,, = A,, and calculate hh c from (93), then evaluate hUt O from (99). level A,,. Whence U can be determined from u

~~utc-huc (k-l)

*

ww

114

R . S . BENsoz% W. J . D. AN)CAND and P. C. BARUAH

T h e n w h e t h e r the flow is choked or n o t the pressure ratio (p~/p~) is d e t e r m i n e d from either (97) or (95). The e n t r o p y Aac is now e v a l u a t e d from (98) and c o m p a r e d w i t h the e s t i m a t e d value of A ~ . I f the two agree w i t h i n the accuracy required, t h e n the calculation is c o m p l e t e ; if not, t h e n the calculation is r e p e a t e d w i t h a n o t h e r e s t i m a t e of A , until t h e required accuracy is obtained. I n f l o w . F o r inflow we assume t h a t the flow is isentropic to the t h r o a t a n d t h a t the t h r o a t pressure is equal to the cylinder pressure. The b o u n d a r y conditions are now represented by the expression

(A,4/(k-1) _ ~ 2 ) ( , ~ , _ A , ) ~ _ ~ _

~b2(A,.O_ 1) = 0

(101)

for subsonic flow where A

=

2

'

Al*n= ~

\-~-~ ]

and

Ao*ut= ~

\-~-!

.

I t is n o t usual to h a v e choked inflow t h r o u g h a valve. F r o m the values of pC, Am and A~, we can d e t e r m i n e A~n. E q u a t i o n (101) is t h e n solved for A* and t h e n A~ut whence ~Out is fully determined. Cylinder thermodynamics

The c o m b u s t i o n products, suffix p, are considered to h a v e a fixed composition a t release. e.v.o. This remains c o n s t a n t u n t i l the air v a l v e opens, a.v.o., w h e n t h e a i r - f u e l m i x t u r e enters t h e cylinder. The a i r - f u e l m i x t u r e mixes w i t h t h e p r o d u c t s w/Shout chemical reaction to give the final products, suffix g. The molecular weighs of the fuel air m i x t u r e mw~, a n d t h e p r o d u c t s mwr are o b t a i n e d f r o m t h e expression (85). The specific heats at c o n s t a n t v o l u m e for t h e air fuel m i x t u r e Cv~ and the products C~r are o b t a i n e d from expression {84). N o w R =- Rm°l mw

Hence Rmo] R~ ~- - - , m wp

(102)

Rmo! R m = --,

(103)

m~v m

whence Cvm = Cv,,, + R m,

(104)

C~p = C~,~+ R v,

(105)

X m ~

(106)

and t h e n u m b e r of mols ~ncm mw m

X~

m e -- m~,,,

(107)

~r~ wp

where m c is t h e t o t a l mass in the cylinder and me~ is the mass of a i r - f u e l m i x t u r e . The t o t a l mass in the cylinder m c = m~g is given b y m~g = Xm m wm+ X~ ra w~ = me.

(108)

The specific heats at c o n s t a n t pressure and volttme will be C~, = X m mwm C~,, + X~ mw~ C~p,

(109)

m~) a

Cvg = X m mwm Cv'' + X ~ mwp C~,,

( 11 O)

A model for spark ignition engine with intake and exhaust systems

115

and viscosity ( X m Zcoeff m 4mw.

Zooeffe =

+ X r Zcoeffr 4mw~)

(X,,, 4mw,,, + X~, ~]m,~p)

---~ Zcoeff e ×

(111)

'

(112)

T0"S4~.

Specific heat ratio:

k¢, = C=,t,/C~,,.

(113)

The mass flows into and out of the cylinder (dm,/d~ and dme/da) are controlled by the boundary conditions at the valves. The cylinder pressure (Pc) is calculated by the first law of thermodynamics in the form

Vc 'dpc k c - 1 da =

kc

dV~

a~o din,

kc--~--lPc-~"-~ k i - 1

da

a~, dme

dQ

(114)

k e - 1 d~ ~ d a '

when there is outflow through the inlet valve, the second term on the right-hand side becomes a~ dm, {din, ) k~--i da ' \ d a is negative . The mass balance in the cylinder is dmc = dmi da d~

dme da '

(mc)~+l = (me).+ (dmdda) An, (Pc)~+l = (Pc).+ (dpdda) Aa and (T~).+ 1 = (Pc VdRmc).+~. The heat transfer rate is calculated using Annand's convective equation (65)* with

t m = tg and C ~ = C~a. The speed of sound ae. is determined from the kee and Reo and the cylinder temperature. Thus, the basic wave action and gas exchange computer programs developed for compression ignition engines can be modified by the adjustments indicated. I n essence this includes the addition to the p a t h line compatibility condition A~, the concentration (X,) of each of the species. These in turn lead to the updating of k for each path line as a function of p, T and X~. The computer logic is more complex and the chemical reaction calculations along path lines are time consuming. I n practice, below certain fixed temperatures the constituents are almost constant. In these cases the computer program by-passes the equilibrium thermodynamic expression and the rate equations. TESTS

ON SIMULATION

MODEL

The model was developed in two parts and then brought together and tested as an integrated whole. The power cycle was developed and tested against work carried out oll a small gas engine on which engine research had previously been carried out on combustion limitations of gaseous fuels. 1~ The gas exchange process had been fully developed for compression ignition engines. This program was modified to include the chemical reactions along p a t h lines as well as the gas exchange process in the cylinder to allow for 12 species. The small gas engine was set up as shown in Fig. 5. No carburetter was fitted, the fuel and air being mixed in the chamber shown. The inlet pipe was 0.600 m long. The exhaust pipe was 0.346 m long with a sampling station at 0.219 m from the engine. The receiver at the outlet end of the exhaust pipe was connected to the laboratory exhaust system, there was therefore a back pressure on the system. The engine bore was 0"105 m, the stroke 0-152 m and the compression ratio 7"7. Tests were carried out at constant speed of 920 rpm and variable air-fuel ratios, as part of an investigation of sampling techniques. 2s We only report here the NO results obtained with the N.D. Infra-red Analyser. The purpose of the calculation tests at this stage were (i) To establish whether the simplified model for the closed cycle gave reasonable prediction of NO in exhaust pipe. (ii) To establish whether the simplified model for the open cycle gave reasonable prediction of the effects of the gas exchange process on NO in the exhaust gases. In the engine under test, the valve overlap was 18 ° CA (EVC 372 °, AVO 354°}.

116

1~. S. B E N S O N , W . J . D . A N N A N D a n d P . C. B ~ u a . K

I f t h e m o d e l is to give a good simulation of t h e engine under t e s t t h e n u m b e r of e s t i m a t e d variables in t h e i n p u t d a t a should be a m i n i m u m . The flame factor f f is one u n k n o w n , which m u s t be based on previous experience. H o w e v e r , to reduce t h e c o m p u t a tional t i m e a good e s t i m a t e of the residuals will help. I t should be pointed out t h a t t h e calculation will a u t o m a t i c a l l y adjust the residuals. I n order to t e s t h o w sensitive the

~kJ~

Fue ~ )

Cyhnder

Samphng stahon

.u...o./®

°

/

Air/

~ Mixing chamber

Receiver

P=pe No.

Length

D=clmeter

(m)

I ]3

TIT

(m)

0 219

0 044

0 127 0'600

0 044 0 044

F1o. 5. F o u r - s t r o k e cycle engine configuration.

60

50

at f..~

4O

// O x E

/

\\

\ \ \

&3o

\ \

! I

•r- 2 o

\

Fraction residual 0 04 ~

10

o 07 . . . .

Flame facfor = 4 0

,\ \

\ \

\ \

i

14

15

16 17 12 19 20 Atr fuel ratio Fro. 6. Effect of fraction residual on nitric oxide prediction. prediction of N O is to t h e residuals a n u m b e r of p o w e r cycles were calculated w i t h t h e flame factor set a t 4.0 based on t h e previous w o r k of Bailey 1~ on t h e s a m e engine. T h e results are shown in Fig. 6. I t is quite clear t h a t t h e p e a k N O is strongly influenced b y t h e residuals. This is n o t u n e x p e c t e d since N O f o r m a t i o n r a t e is d e p e n d e n t on t h e m a x i m u m t e m p e r a t u r e a n d this is influenced b y t h e internal energy of t h e inert gases. The greater t h e qnAutity of inert gases t h e lower t h e t e m p e r a t u r e rise during combustion. The flame factor f f o f 4.0 was based on B a i l e y ' s w o r k 17 which did n o t include w a v e action calculations. W e therefore carried out a n u m b e r of calculations t o establish t h e

A m o d e l for spark ignition engine w i t h i n t a k e and e x h a u s t systems

117

influence of this factor on t h e N O prediction w h e n w a v e action is included a n d t h e b a c k pressure in t h e receiver set to t h e e x p e r i m e n t a l value. T h e initial v a l u e of t h e residuals was selected a t 0.055 a n d f f set a t 3.0, 3"15, 3'5 and 4.0. The calculation was carried o u t for some 8 consecutive cycles (16 revolutions} for each v a l u e o f f f i A t y p i c a l result is shown in Fig. 7. I t will be seen t h a t t h e residuals soon settle to a c o n s t a n t value and at t h e

o 25,o

&

E ¢o --

2490

~ 13_

-

4200

~E 4000

--104

®E m

--102

-0060 ~ --0056

I

I

I

I

I

I

2

3

4

5

6

7

u. =

8

NO. of c o m p u t e d c y c l e s Eqmvalence ratio=09

Flame factor=315

Fro. 7. N u m b e r of cycles to s t e a d y state. eighth cycle t h e N O emission at t h e sampling point is constant. I n Fig. 8, t h e N O results a t t h e eighth cycle a t t h e sampling p o i n t are shown. H e r e it is clearly seen t h a t t h e factor f f has a strong influence on the results. B y c o m p a r i n g the results at an equivalence ratio of 0-9 w i t h t h e experiments, a flame factor of 3.15 was selected as t h e base for the speed of 920 r p m . Calculations were carried out o v e r an equivalence range of 0.8-1.1. The results for t h e e i g h t h cycle at t h e sampling p o i n t are shown in Fig. 9. On t h e same g r a p h t h e e x p e r i m e n t a l results are also shown. I t will be seen t h a t there is good a g r e e m e n t b e t w e e n e x p e r i m e n t s and calculations. I t m a y be deduced, therefore, t h a t t h e i n t e g r a t e d m o d e l gives good predictions of t h e N O emissions. The mass flow predictions were also v e r y good. H o w e v e r , no m e a s u r e m e n t s were m a d e of cylinder pressures or power so we are not able to c o m p a r e these results. A further p a p e r will e x a m i n e these parameters, e I n Fig. 10 t h e calculated influence of receiver b a c k pressure on t h e emissions is shown. I t will be seen t h a t emissions are reduced as t h e b a c k pressure is increased. This is due to the increase in residuals in the cylinder which reduces t h e p e a k t e m p e r a t u r e . All t h e results show t h a t t h e calculated emissions at the sampling point for N O were v e r y n e a r l y t h e same as those calculated in the cylinder. This was due to t h e short distance b e t w e e n t h e sampling p o i n t and t h e cylinder a n d t h e low t e m p e r a t u r e s . Details of t h e calculation (Fig. 11) show trends which h a v e been r e p o r t e d b y other workers. This gives confidence in the modelling. The t r e n d of flame speed w i t h air fuel ratio (Fig. 9) agrees w i t h a n u m b e r of published results. ~9 The increase in flame radius w i t h t i m e (Fig. l l ( a ) ) agrees w i t h results given b y Bailey z? and I i n u m a . 3° The v a r i a t i o n of flame speed w i t h angle (Fig. 11 (b)) agrees w i t h test results of I i n u m a . 8° The prediction of m a x i m u m N O a n d its relationship to t h e t i m e for freezing (Fig. ll(c)) agrees well w i t h t h e w o r k of L a v o i e et al. ?

R . S. B ~ . ~ s o N , W . J. D . ANNAND a n d P. C. B A ~ U ~

118

The method is now being extended to four-cylinder four-stroke engines and to twostroke crankcase compression engines. I t is hoped to report later on the four-cylinder work, the results from the first stage of the crankcase compression engines are given in another paper. 6

Equivalence r a t i o - 0 9 6O

Com

'

=

/ ( ,

~o x

u 40

~0 2700

.~ o 2600 .._6...-- A ' ' - - - 2400

,~ 150 E .

1:30 I10

o

--

90

30 55 Flame factor

40

F I o . 8. I n f l u e n c e o f flame factor.

2700

l

!

Flamefactor =3.15

28001-i I ~,...~I

2500

E

2300

-~

P~ak temperal"ure 0.

24O0

b M X E

.91 P

.~1 Temperature at which *9°1 / n i t r i c oxide freezes

2ooc

~3o

(hi

~

i~

f - , F J-

18(X?

160C -

I ~ , ~ l a m e

20

15 E

speed

Z

14

I

15

I II

16

I

17

18

19

I

10

I I

I 7 ~ U..

[NO] Cm'nputed •

I

20

I

Air fuel ratio 14

I

15

II

Flame factor :3 15 I

I

~6 17 t8 A i r fuel r a t i o

FIG. 9. Influence of alr-fuel ratio.

•

Experiment

I

19

I

20

A model for spark ignition engine with intake and exhaust systems

119

2500

t

2480 420C

~

4000

Z •3 8 0 0 0 064 0062

~

0 060

~0058 0056

I

~0

H Back pressure,

T

1.2 bar

FIG. 10. Influence of back pressure.

45

~D

~ouE --

~ol-- Abp . . . . t delay: 3~/z:CA | ( F s t i m a t e d by o x t r a p o i a t l o n ~ 81"- of established part ~ | curve on to tlmo /

" 61-"i" 4L ~ /

~

"1~ ;

L.

2LI

oL~l 0 I]

~

/II

]

r.pm 920 I Eq. . . . Ionc, rotio091 Flame factor 3 15 I Ignillon 15e BTDC I

/ I0

li If 3o1~ t l

o

I

I

I

30

20

40

I

I

50

60

Travel 1"line from ignH'ion~

"CA

-12~oo I

',,

~

-

i.J

I

7~

) ) ~ " ~

',\ ', \

2700

I

/

/

(a) I

(c)

~

Computed flame speed

14

1500 -

12 I::

1300 -

lo 3~ 6o 85 Iio 135 Crank angle, ATDC (d)

6 th cycle

,,

~,] ,,'~ i.o 900 •o o

4

4o

~- i- 0 5 /// ' t / 0 705

720

I 15

Crank angle,

~ 46 ~_ Two zone termi~ Ot'On'~'~ 2 30

45

~ 1--

5oo-

~

300

~L

--

2

I 138 EVO

ATDC

FxG. 11. Performance prediction.

I 360 C r a n k angle r

ATDC

,572 AVC

120

R.S. BENSON, W. J. D. ANNAND

and P. C. BARUAH

CONCLUSIONS A comprehensive simulation model of a four-stroke cycle spark ignition engine is presented. The model combines a full power cycle simulation with the prediction of NO emissions with a comprehensive gas d y n a m i c model for the cylinder a n d ducts allowing for chemical reactions in the e x h a u s t pipe. Calculations are presented comparing the predicted NO with test results from a single cylinder engine. I t is shown t h a t good a g r e e m e n t is o b t a i n e d between the predicted a n d measured NO over an equivalence range of 0.8-1-1 when the flame speed is corrected at the equivalence (0.9) corresponding to the p e a k NO. The NO emission is sensitive to the residuals a n d the e x h a u s t b a c k pressure b o t h of which are a u t o m a t i c a l l y included in the overall simulation model. F u r t h e r w o r k is required to o b t a i n a b e t t e r q u a n t i t a t i v e u n d e r s t a n d i n g of the factors which influence the flame speed. I n the m e a n t i m e the selection of an empirical f a c t o r f f which when multiplied b y the l a m i n a r flame speed gives the t u r b u l e n t flame speed is suggested. The p r o g r a m written from this simulation model can be used to assist in the design of intake a n d e x h a u s t manifolds as well as to s t u d y m a n y problems such as location of t h e r m a l reactors, catalytic devices, e x h a u s t gas recirculating valves; effects of misfire a n d maldistribution of fuel-air mixture. L a t e r papers on these applications will be presented as the research progresses. Acknowledgements--The authors wish to acknowledge with thanks the assistance of

Mr. I. Milne who made available the results of his experiments on the National Gas Engine, Mr. J. Nicholson for his assistance in the computing work and Dr. B. J. Tyler for his advice on NO measurements. The work has been carried out in both of the Mechanical Engineering Departments of the University of Manchester and U.M.I.S.T. with their funds. It is now being sponsored by the Science Research Council. All the calculations were carried out on the U.M.R.C.C.C.D.C 7600 computers. We acknowledge with thanks the assistance of the Director and his staff.

REFERENCES 1. R. S. BENSON,S.A.E. Auto. Engng Cong., Detroit, Michigan, Paper No. 710173 (1973). 2. R. S. BENSO~ and P. C. BA~UAH, S.A.E. Combined Vehicle & Powerplant Meetings, Chicago, Paper No. 730667 (1973). 3. J. B. HEYWOOD, S. M. MATHEWSand B. OWEN, S. A. E. Paper No. 710011 (1971). 4. J. B. HEYWOODand J. C. KECK, Environ. Sci. Tech. 7, 216 (1973). 5. R. J. TABACZYNSKI,J. B. HEYWOODand J. C. KECK, S.A.E. Paper No. 720112 (1972). 6. R. S. BENSON, P. C. BARUAHand B. WHELAN, ~. M~ch. Engng. 7. G. A. LAVOIE,J. B. HEYWOODand J. C. KECK, Combustion Sci. Tech. 1, 313 (1970). 8. C. W. VICKLAND,F. M. S~aANGE, R. A. BELL and E. S. STARKMAN,S.A.E. Trans. 70, 785 (1962). 9. W. J. D. ANNAND,Proc. I. Mech. E. 177, 973 (1963). 10. D.B. SPALDING,P. L. STEP~NSON and R. G. TAYLOR,Combustion Flame 17, 55 (1971 ). 11. K. AKITA,Int. Chem. Engng II, 739 (1971). 12. A. S. SOKOLIK,V. P. KARPOVand E. S. SE~NOV, Combustion, Explosion Shock Waves 3, 36 (1967). 13. J. CHOMIAK,Combustion Flame 15, 319 (1970). 14. N. M. HOWE and C. W. SHIPMAN, Tenth Symp. (Int.) on Combustion, p. 1139 (1964). 15. R. A. PHILLIPS and P. L. ORM~, Advances in Automobile Engineering, Vol. 4, Pergamon Press, Oxford (1966). 16. E. MALLARDand H. LE CHATELIER,Annales de Mines 4, 274 (1883).

A model for spark ignition engine with intake and exhaust systems

121

17. A. C. BAILEY, Ph.D. thesis, University of Manchester (1971). 18. D. K. KUEHL, E i g h t h Syrup. (Int.) on Combustion (1962). 19. R. S. BENSO~, A d v a n c e d Engineering T h e r m o d y n a m i c s . Commonwealth and International Library. 20. H. K. NEWHAI~, Twelfth Syrup. (Int.) on Combustion (1968). 21. W. J. D. ANNA~D, J . Mech. E n g n g Sci. 12, 146 (1970). 22. R. S. BENSON, 1~. D. GARG and D. WOOLLATT, I n t . J . mech. Sci. 6, 117 (1964). 23. R. S. BE~SON, A. WZLD and D. WOOLLATT, Proe. I . Mech. E . 184, 276 (1969-70). 24. R. S. BE~SON, Report No. 1 for C.E.G.B., U.M.I.S.T. (1970). 25. R. S. BENSON, I n t . J . mech. Sci. 14, 635 (1972). 26. E. JENNY, Doctoral thesis, E.T.H., Zurich (1949). 27. R. S. BENSON and K. GALLOWAY, Proc. I . Mech. E . 183, 253 (1968-69). 28. I. L. MILNE, M.Sc. diss., U.M.I.S.T. (1971). 29. W. J. D. ANNA:~D, U.N.I.C.E.G. Discussion Draft (1971). 30. K. IINUMA, J a p . A u t o . RCS. I n s t . Tech. Memo, No. 2 (1971).

APPENDIX NUMERICAL

PROCEDURE EQUILIBRIUM

I

FOR THE SOLUTION EQUATIONS

OF

The equilibrium equations (25)-(31) can be rewritten as : - (1)

K . = K~,,/4p,

(AI)

(2)

K b -- K~,,/~lp,

(A2)

(3)

K c = K,,./~]p,

(A3)

(4)

Ka = K , , / p ,

(A4)

(5)

K. = K,,/~]p,

(A5)

(6)

K t = K,,,

(A6)

(7)

Kg = K,,/.,]p,

(A7)

(8)

Xx,

= 1.o.

(A8)

The atomic balances are (9)

argon

(10)

carbon

(11)

hydrogen

(12)

oxygen

(13)

nitrogen

Xl~ = av,

(A9)

X 8 + X 9 = aw,

(A10)

2X1+ 2X~ + X 3 + X 4 = ax,

(All)

Xx + X3 + X6 + 2X8 + X9 + 2Xlo + X l l = ay,

(A12)

2 X ~ + X s + X 7 = az,

(A13)

where v, w, x, y and z are the number of atoms of A, C, H, O and N corresponding to 1 mole of fuel plus the air, 1 mole of product arises from a mole of fuel plus 1/¢ stoichiometric quantity of air. The procedure of solution of these 13 equations is based on Vickland s modified to reduce computational time. The calculation starts with the preparation of the equilibrium constants in terms of Ka, Kb, Kc, Kd, Ke. K 1 and Kg. For one particular state (p, T) these are fixed. Values of a and b are selected and all the values of Xj are evaluated and the equations (AS) and (A12) are checked with these calculated values for balance. I f they do not balance within a stipulated accuracy a N e w t o n - R a p h s o n adjustment is made to a and b and calculations repeated. F r o m a large number of calculations the following are suitable values for a and b for the first estimate

122

R. S. B~.~so~, W. J. D. AN~A~D and P. C. BAttUAH

Value of a When ¢~> 1.0,

a=

1.3 [ca + 0.5hy + 1.863(2ca + 0.Shy - ox)/¢] exp (0.13T/1000)

(A14)

and when ~ < 1.0, 1.3 a = [0.25by + 2.363(2ca + 0 . S h y - ox)/~ + 0.5ox] exp (0.13T/1000)

(AI5)

Value of b

When temperature T (in °K) > 3000 b = exp ( 1 0 . 3 - ( 3 . 1 - 0 . 1 7 logP) T/IO00)

(A16)

and when T ~<3000 b = exp ( - 9 . 0 + 0 . 5 1 o g P + 3 0 0 0 0 / T ) .

(A17)

Adjust b for the following case B X -- 2 . 0 - 9 . 0 1 o g ~ .

I f B X > 3.5, B = exp (exp (3.5)+0.25 logP)

(A18)

B -- exp [exp (2.0-9.0 log ~) T 0.25 log P]

(A19)

b -- B,

(A20)

I f B X <~3.5, I f b > B , then use where ca is the n u m b e r of carbon atoms in the fuel, hy is the n u m b e r of hydrogen atoms in the fuel, ox is the n u m b e r of oxygen atoms in the fuel, ~b is the equivalence ratio of fuel (actual fuel-air ratio]theoretical fuel-air ratio), P is the pressure in atmospheres and T is the temperature in °K. Having fixed the values of a and b, X 2 can be calculated by using the equilibrium constants of equations (25) and (29) in equation ( A l l ) which will be 2 b X 2 + 2X~ + K~ b ~/Xz + K~ ~X~ = ax.

(A21)

Then follow the calculation of X 1, X8 and X 4. From (27) and (31) inserted in (AI3), we have 2 X 5 + Kg b ~ X 5 + Kc ~]X6 = az, (A22) whence X 5 followed b y Xs and X~. A combination of equations (30) and (A10) gives X s + K I Xs/b = aw,

(A23)

whence X s, X , and X10 are evaluated from equation (28) and X l l from equation (26). X12 is given directly b y equation (Ag). These values are now tested in equations (A8) and (A12), neither of which has so far been used. I f they do not balance within a stipulated accuracy the calculations are repeated with new values of a and b adjusted b y the N e w t o n Raphson method, in the following procedure. Equation (A8) is written as es = ~ X s - 1. (A24) Equation (A12) is written as er =

X1 + Xa + 2Xs + Xs + 2X1° + X l l + X° a - . Y

(A25)

The right-hand sides of the above equations are functions of a and b only, therefore, we write er ----F(a, b), (A26) es = G(a, b).

(A27)

A model for spark ignition engine with intake and exhaust systems

123

The required solution is obtained when er = 0 ,

(A28)

es = 0 .

(A29)

The solution starts with approximate values of er and es where

F(an, bn),

(A30)

es = O(an, bn).

(A31)

er =

The new values of the variables are adjusted for the next step b y an+l --~ an-- Aa,

(A32)

bn+l = bn-Ab,

(A33)

where Aa = er (O]~b)(es)--es (~/~b) (or)

J(an, bn) Ab = er -- (O/tOa)(or) Aa

(~/~b) (or) and

J(a., b.) -- ~a (er) ~-~( e s ) - ~ (er) ~a (es). The above procedure will be repeated for new values of an+l and b~+1 until the solution within stipulated precision is obtained.

APPENDIX RATE FORMATION

II

FOR NITRIC

OXIDE

The rate formation for nitric oxide given b y equation (43) is derived as follows. Let K~I denote the forward rate constant for the i t h reaction, K~b the backward rate constant for the ith reaction and R~ the "one-way" equilibrium rate for the i t h reaction. Then

(NO)

(N)

( N O ) . = ~''

~

(N,O) = ~°'

(N,o)----: = y ' ;

where the suffix e represents equilibrium and V is the volume of burned gas.

Calculation of (I/V) (d/d0 ((NO) V) From equation (36) the net rate is -

KxI(N) (NO) + Klb(N~) (O) = - ~ 8. Klt(N)o (NO)e + Klb(N2), (O)e.

But KlI(N)~ (NO)~ = Klb(N~). (O)o = R1, so t h a t the net rate becomes - a, f~, R 1+ R 1. Using the similar terms for equations (37), (38) and (41), involving NO, we finally have sa

((NO) V) =

-e,(~,RI+R2+Ra+2c~,R.)+RI+[3.(R~+Rs)+2y, R..

(A34)

Calculation of (l/V) (d/dr) ((N) V) Considering equations (36)-(38) and using similar procedure we have 1 d

d-~((N) V) --- -ft.(a. R 1+ R2 + R,) + RI + a.(R, + Ra).

(A35)

124

R . S . Bv.NSON, W. J. D. AN~AND and P. C. BA-RUAH

C a l c u l a t i o n o f (l/V) (d/dr) ((N20) V)

Considering equations (39)-(42), we have ((N20) V) = - y , ( R 4 + R s + R 6 + R ~ ) + R 4 + R ~ + a ~ R ~ + R

(A36)

7.

I t has been foun& t h a t the relaxation times for equation (A35) and equation (A36) are severM orders of magnitude shorter t h a n t h a t of equation (A34). Therefore, steady state m~y be assumed for (N) and (N20) which means the right-hand sides of equations (A35) and (A36) can be sot equal to zero. Then from equation (A35) Rt + ae(R2 + R3)

(A37)

[3e = a ~ R I + R 2 + R a

and from equation (A36)

R4 + Rs +or~.Rs + R~ Y~'=

R~+R6+Re+R~

(A38)

"

Using these values of fl~ and Y0 in equation (A34) we have Rx

R,

)

1 d ((NO) V) - 2(1 -- a~) (1 + ot.[Rz/(R, + Rs) ] + 1 + [ R e / ( R 4 + R 5 + R~)],, ' "V which is the final rate equation for NO.

(A39)