A non-equilibrium analysis of the Otto cycle engine from the viewpoint of combustion and emissions.


Material Information

A non-equilibrium analysis of the Otto cycle engine from the viewpoint of combustion and emissions.
Physical Description:
xv, 335 leaves. : ill. ; 28 cm.
Srivatsa, Shesh Krishna, 1948-
Publication Date:


Subjects / Keywords:
Internal combustion engines -- Combustion   ( lcsh )
Automobiles -- Motors -- Exhaust gas   ( lcsh )
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )


Thesis--University of Florida.
Includes bibliographies.
Statement of Responsibility:
Shesh Krishna Srivatsa.
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 000582660
notis - ADB1038
oclc - 14158355
System ID:

This item is only available as the following downloads:

Full Text








I would like to express my thanks to Dr. V. P. Roan, the

Chairman of my supervisory committee, for suggesting this study,

and for his guidance and constructive criticism during the course

of this work. Thanks are also due to Professors R. B. Gaither,

R. K. Irey, D. Kim, and U. H. Kurzweg for serving as members of

the supervisory committee. Finally, the help of Mrs. Jeanne Ojeda

in typing the manuscript is gratefully acknowledged.



ACKNOWLEDGEMENTS ................................................. ii

LIST OF TABLES..................................................... vii

LIST OF FIGURES...................................................viii

KEY TO SYMBOLS ..................................................... x



I INTRODUCTION........................................... 1

EMISSIONS............................................... 3

2.1 Introduction... ..................................... 3

2.2 Crankcase Blowby............................. ...... 3

2.3 Evaporation Losses.......................... ....... 4

2.4 Formation of Pollutants............................ 5

2.4.1 Hydrocarbons........................... ... ... 5

2.4.2 Carbon monoxide ............................. 7

2.4.3 Oxides of nitrogen.......................... 7

2.5 Effects of Various Parameters on Engine Emissions.. 8

2.5.1 Hydrocarbons ................................. 8

2.5.2 Carbon monoxide............................. 9

2.5.3 Oxides of nitrogen.......................... 10

2.6 Emission Standards................................. 10



2.7 Control Methods................................. 11

2.8 Description of the Combustion Process in an
Otto Cycle Engine............................... 13

2.9 Laminar Flames: Introductory Remarks........... 14

2.10 Closing Remarks................................ 16

III LITERATURE REVIEW..................................... 17

References for Chapter III............................. 33


4.1 Formulation of the Problem...................... 37

4.1.1 The conservation equations................ 37

4.1.2 Major assumptions for analysis........... 40

4.1.3 Boundary conditions...................... 43

4.2 Method of Solution .............................. 50

4.2.1 Initialization........................... 53

4.2.2 Details of computation................... 53

4.2.3 Step size and stability.................. 62

4.3 Chemical Kinetics............................... 67

4.4 Uncertainty in the Definition of Temperature.... 71

References for Chapter IV............................ 76


5.1 Introduction.................................... 78

5.2 Present Work..................................... 81

5.3 Major Assumptions for Analysis.................. 83



V 5.4 Theoretical Analysis............................

1. Calculations of conditions at beginning
of compression................. ...............

2. The compression process......................

3. Combustion................. ..................

4. Expansion stroke analysis....................

5.5 The Kinetics of Nitric Oxide Formation...........

5.6 Details of Calculations for Nitric Oxide
Formation ........................................

5.7 Injection of Water-Alcohol Mixtures..............

References for Chapter V..............................

VI RESULTS AND DISCUSSIONS................................

6.1 The Combustion Process............................

6.1.1 Concentration and temperature profiles....

6.1.2 Unity Lewis number approximations.........

6.1.3 Two-stage hydrocarbon combustion..........

6.1.4 Closing remarks.................... ........

6.2 Mathematical Model for the Otto Cycle............

6.2.1 The conibustion process....................

6.2.2 The expansion process.....................

6.2.3 Water-alcohol injection...................

References for Chapter VI.............................






























EQUILIBRIUM CALCULATIONS...........................

References for Appendix I..........................


References for Appendix II.........................

COMPUTER PROGRAMS..................................

References for Appendix III........................











Table: Page

4.1 Reactions considered in this analysis.................... 69

5.1 Reaction mechanism....................................... 103

6.1 Results for water-alcohol induction ...................... 151


Figure: Page

4.1a Staggered mesh scheme........... ....................... 51

4.1b Enlarged view of staggered mesh........................ 51

5.1 p-V diagram for the Otto cycle ......................... 79

5.2 Chemical reaction times for hydrogen-air and
hydrocarbon air mixtures at one atmosphere and unity
equivalence ratio...................................... 92

5.3 Schematic diagram illustrating the sequence of
processes occurring during combustion................... 94

5.4 Mass fraction of fuel-air mixture burnt as a func-
tion of crank angle..................................... 96

5.5 Experimental and theoretical nitric oxide formation
rates for combustion of hydrogen-air mixtures........... 112

5.6 Experimental and theoretical nitric oxide formation
rates for combustion of propane-air mixtures............ 113

6.1 Temperature profile at different instants of time
illustrating the development of the steady-state
profile................................................ 123

6.2 Concentration profiles for the species CEH, CO, CO2,
H 2 and 02 ............................................ 124

6.3 Concentration profiles for the species H2, OH, NO,
NO2, and N 20 .......................................... 125

6.4 Concentration profiles for the species H and 0......... 126

6.5 Concentration profiles for the species CH3, CHO, and N. 127

6.6 Test of unity Lewis number approximation for methane
and carbon monoxide in a methane-air flame.............. 133

6.7 Test of unity Lewis number approximations for methane
and carbon monoxide in 0.1 atm. methane-oxygen flame... 134



Figure: Page

6.8 Steady-state temperature profiles for unity and
non-unity Lewis numbers............................... 135

6.9 Pressure and temperature histories during combus-
tion process............................................ 139

6.10 Nitric oxide formation during combustion................ 140

6.11 Pressure-time and temperature-time plots for the
expansion stroke....................................... 142

6.12 Concentration-time histories of the species CO, H,
O, OH, and NO during the expansion stroke.............. 144

6.13 Concentration-time history of species N during the
expansion stroke ....................................... 145

6.14 Concentration-time history of species 02 during the
expansion stroke....................................... 146

6.15 Concentration-time histories of the species H2 and
NO2 during the expansion stroke........................ 147

6.16 Concentration-time histories of the species CO2 and
H20 during the expansion stroke........................ 148


a = Velocity of sound

A = Preexponential (frequency) factor in Arrhenius rate constant

C = Specific heat at constant pressure

D.. = Binary diffusion coefficient for species i and j

DT,i = Thermal diffusion coefficient for species i

E = Activation energy

f = Mass fraction of residual burnt gases in the engine cylinder

f. = External force per unit mass on species i

h = Heat transfer coefficient

h. = Specific enthalpy of species i

h = Standard heat of formation per unit mass of species i at
temperature To (= 298 degrees K)

k = Boltzmann's constant

k = Forward rate constant

kb = Backward rate constant

K = Equilibrium constant

H = Mach number; number of elementary reaction steps

n = Number of moles per unit volume

n. = Molar concentration of species i

N = Number of chemical species present

N. = Molar flux of species i-= n.V.

p = Thermodynamic pressure

p = Pressure tensor
q = Heat flux vector

Q = Heat transfer

r = Radius

R = Universal Gas Constant

s = Entropy

t = Time

T = Temperature

v = Mass average velocity of the gas mixture

V. = Diffusion velocity of species i

w. = Rate of production of species i by chemical reaction

W. = Molecular weight of species i
S= ole fraction of species i

X. = Masole fraction of species i

Y. = Mass fraction of species i

U = Internal energy

y = Shear viscosity

A = Thermal conductivity

V.i = Stoichiometric coefficient of species i in jth elementary
reaction step

0 = Collision integral

= Intermolecular potential function

p = Density

0 = Collision diameter

6 = Crank angle


b = Burnt gases

BDC = Bottom dead center

f = Flame element

TDC = Top dead center

u = Unburnt gases


o = Designates standard state

- = Written over a symbol designates an average quantity

P = Product species

R = Reactant species

Abstract of Dissertation Presented to the Graduate Council
of the University of Florida in Partial Fulfillment of the Requirements
for the Degree of Doctor of Philosophy



Shesh Krishna Srivatsa

August, 1974

Chairman: Dr. V. P. Roan
Major Department : Mechanical Engineering

An analytical study of the thermodynamic processes occurring

in an Otto cycle engine has been conducted with a view to obtaining a

better understanding of the combustion process and a realistic

prediction of the emission levels. Much of the work previously

reported in this area has been experimental, the reason being the

extreme complexity of any realistic solution to the governing equa-

tions. The availability of high-speed computing machines provides the

necessary tool for such an approach. With such machines, theoretical

calculations have assumed a new importance and in many cases are being

used to obtain estimates which would require extensive experimental

work. Theoretical analysis, by its very nature, involves the use of

simplifying assumptions to reduce any real problem to a stage where

it is solvable. Any unrealistic simplifying assumption will reduce

the solution to one of academic interest only and the literature on

any subject is replete with such examples. In this study, an analysis

with a minimum of simplifying assumptions has been carried out and the

results obtained are realistic and usable.


The study presented here has been divided into two parts.

Firstly, a study of the combustion process has been carried out by

numerically solving the set of partial differential equations govern-

ing one-dimensional unsteady laminar flame propagation in a homogeneous

fuel-air mixture. Secondly, a mathematical model for the Otto cycle

has been developed to predict the thermodynamic properties and the

species concentrations throughout the cycle based on a non-equilibrium

analysis and finite rate combustion.

For the first part, an explicit finite difference method has

been used to numerically solve the governing equations and thereby

obtain a detailed description of the flame structure in terms of the

distribution of the state properties and species concentrations through

the flame zone. Most of the previous analyses have ignored the varia-

tions in the thermodynamic and transport properties with temperature,

pressure and chemical composition and have also treated the chemical

kinetics as being governed by a single overall reaction. These

particular simplifying assumptions (particularly the last one which

is rather unrealistic for any system of practical interest) have been

removed in the computational scheme presented here. The scheme is

general and considers variable transport and thermodynamic properties

and a realistic reaction mechanism in terms of a set of elementary

reaction steps. The validity of the procedure has been demonstrated

by carrying out calculations for the methane-air flame and obtaining

good agreement with the limited experimental data available.

The mathematical model of the Otto cycle gives a detailed

description of the cycle in terms of the thermodynamic properties

and species concentrations at any stage in the cycle, including the


exhaust. The combustion process has been considered to proceed at

a finite rate instead of the commonly made assumption that this

process is instantaneous. Variations in the thermodynamic properties

with temperature and chemical composition have been allowed instead

of assigning them constant values. An accurate prediction of the

concentrations of the various species at exhaust was the main

objective in the development of the model.

The model has been used to study the effect of the intake

of an additive mixture of water and alcohols (methanol and ethanol)

in different proportions and amounts on the emission levels. The

model is quite flexible and permits the study of the influence of

other factors such as exhaust gas recirculation, air-fuel ratio,

compression ratio, speed, ignition timing, etc., on the performance

and emission levels.


The role of the internal combustion engine as a major

contributor of the hydrocarbon gases, carbon monoxide, and oxides

of nitrogen that contaminate the atmosphere has been recognized for

some years. The formation of these pollutants is closely related

to the combustion process and a realistic modeling of this process

would help the study of the methods to reduce the emissions. A

knowledge of the transient mechanism of high temperature oxidation

of hydrocarbon fuels would lead to a better understanding of the

combustion process and may ultimately help in the design of engines

which will emit lower concentrations of pollutants. Nhile the amount

of unburnt or partially burnt products in the exhaust is generally

insignificant from an energy point of view, it is significant, and

very much so, when considered from an environmentalist's point of

view. Thus these studies leading to new designs might or might not

improve the overall efficiency of the internal combustion engine (in

fact a reduction in efficiency might be achieved), but they might

help to reduce the emissions.

The high temperature oxidation mechanism of some simple fuels

like hydrogen and carbon monoxide have been widely studied and are

fairly well understood, but the state of development for the oxidation

of hydrocarbon fuels leaves much to be desired. Whatever little work

has been reported is for lower hydrocarbons (like methane, etc.) and

kinetic data and reaction mechanism postulations for some practical

fuels (like iso-octane) are virtually non-existent.

One of the methods available for the quantitative determination

of chemical kinetic constants is a detailed analysis of the flame

structure. A comparison of the flame structure predicted theoretically

on the basis of a postulated reaction mechanism, with the experimentally

measured data, provides a convenient way of checking the validity of

the reaction mechanism. The influence of varying the rate constant

parameters on the theoretical predictions permits the determination of

influence coefficients for the different input parameters. A parametric

study of this type is useful in arriving at some sort of a conclusion

without much trial and error.

To put this principle into practice, there remains the need

for a general theoretical computational procedure, and it is the aim

of this study to proceed in this direction by presenting a fairly

general computational scheme for a realistic reaction mechanism.

In addition to the study of the combustion process, a mathematical

model for the entire Otto cycle has also been presented. The model

can be used to study the effect of different parameters on the exhaust



2.1 Introduction

In this chapter, the formation of the various pollutants, the

factors affecting them, and the various control methods available are

briefly discussed. Finally, a description of the combustion process

together with some introductory remarks on laminar flames is provided.

There are three important sources of emissions from an auto-

mobile: the engine exhaust, the crankcase blowby, and evaporation

losses from the fuel system vents, the engine exhaust emitting by far

the largest amount of pollutants. Thus most of the research on

automobile pollution has been concerned with reducing the exhaust

emissions. The role of the other two sources of emissions is discussed

briefly for the sake of completeness.

2.2 Crankcase Blowby

This refers to the gases escaping from the engine cylinders

into the crankcase through the gaps between the sealing surfaces of the

piston and the cylinder walls. This leakage which occurs mainly during

the expansion and the compression strokes increases with increasing

load (increased pressures and mass flow) and typically consists of

about 85% unburnt fuel-air mixture and 15% combustion products. In

terms of hydrocarbon concentrations, this is about 6000-12000 parts

per million. On uncontrolled engines the crankcase emissions account

for about 20% of the total emissions from an automobile.

Crankcase emission control systems are at present in a well-

developed state, so that this source of emissions has been virtually

eliminated. All the control systems use the same basic approach:

recycling of the gases from the crankcase back into the combustion

chamber, and this is commonly referred to as 'positive crankcase

ventilation' (PCV).

2.3 Evaporation Losses

These losses occur from the carburetor and the fuel tank, and

on uncontrolled engines account for about 18-20% of the total hydro-

carbon emissions. Losses from the carburetor occur mainly after the

engine is shut off when the temperature of the fuel system rises due

to the heat stored in the engine. The losses depend upon the operating

conditions immediately before engine shutdown, the fuel volatility,

carburetor bowl capacity, and the ambient temperature.

The losses from the carburetor are generally small when the

engine is running, and the fuel tank accounts for the majority of the

evaporation losses. The losses from the fuel tank are dependent on

the fuel volatility, its temperature, what fraction of the tank is

filled with fuel, and the ambient temperature. Fuel evaporation loss

control devices store the vapors from the fuel tank and the carburetor.

The vapors are subsequently burnt in the engine combustion chamber. The

various controls in use are based on different modifications of this

basic principle. Evaporation loss controls have been in effect since


As mentioned before, these two losses (crankcase and evaporation)

contribute about 20% each to the total emissions of an uncontrolled

vehicle. The controls developed for minimizing the pollutants emitted

thereby are in a fairly good state of development so that these two

sources of emission have been virtually eliminated. All modern auto-

mobiles (1972 and later year models) are equipped with these controls

so that presently the engine exhaust accounts for almost all the


2.4 Formation of Pollutants

We shall now look at the formation of the various pollutants

which occur in the exhaust--our primary concern here.

2.4.1 Hydrocarbons

The presence of hydrocarbons in the engine exhaust shows that

combustion failed to occur in certain parts of the engine cylinder.

Calculations (chemical equilibrium or finite rate) reveal that the

concentrations of the hydrocarbons in high temperature combustion

products are negligibly small. The source of the unburnt hydro-

carbons is not this fairly homogeneous part of the combustion chamber.

The major source of the unburnt hydrocarbons is the quench layer

adjacent to the relatively cool combustion chamber walls, wherein the

oxidation of the fuel is only partially completed. Further, due to

improper mixing, the local mixture composition may not be within the

flammability limits and this would partly account for the unburnt

hydrocarbons. In an Otto cycle engine the quenching is by far the

more important of the two.

The quench layer may be considered to consist of two zones:

one known as the total quench zone which is immediately adjacent to

the wall and where virtually no chemical reactions occur, and the

other known as the partial quench zone where the temperature is high

enough for cool flame and slow oxidation reactions to take place, but

not at a fast enough rate to cause complete combustion.

After flame propagation, some of the hydrocarbons in the quench

layer mix with the burnt gases and experience partial oxidation. The

availability of oxygen and the degree of turbulence strongly influence

the extent of this oxidation. Also all the hydrocarbons present in

the quench layer are not expelled during exhaust and some remain as

residuals in the cylinder for the next cycle. Reactions continue in

the gases expelled into the exhaust pipe causing further oxidation of

the unburnt hydrocarbons. For these reasons, the amount of hydrocarbons

actually present in the exhaust gases is considerably less than the

hydrocarbons present in the initial quench layer.

Besides the quench layer the crevices between the piston rings

and other dead spaces and pockets in the combustion chamber which prevent

flame propagation also contribute to the exhaust hydrocarbons.

The hydrocarbons in the exhaust are a complex mixture of

oxygenated and partially oxidized compounds together with products

produced by the cracking of the fuel; thus the exhaust hydrocarbons

are different in composition from those produced by crankcase blowby

and evaporation losses which contain mainly the fuel itself.

2.4.2 Carbon monoxide

This is a product of incomplete combustion and results when

sufficient oxygen is not available for combustion, either on an overall

or a local basis. Significant amounts of carbon monoxide may be present

in the combustion chamber immediately after the combustion is complete

even when sufficient oxygen is available. This can be shown by a

chemical equilibrium analysis. However, as the gases cool during the

expansion stroke, much of this carbon monoxide gets converted to carbon

dioxide, so that except for fuel-rich mixtures, the amount of carbon

monoxide in the exhaust may be of the order of about 1% by volume or

about 10,000 parts per million. The conversion of carbon monoxide to

carbon dioxide during the expansion stroke cannot be predicted from

chemical equilibrium considerations since the chemical reactions

involved are not rapid enough to reach equilibrium. As such, the actual

concentration of carbon monoxide at exhaust is higher than that if

chemical equilibrium were to prevail.

2.4.3 Oxides of nitrogen

A brief discussion is provided here; details may be found in

Chapter V. The most important of the oxides of nitrogen of concern

here is nitric oxide (NO). The formation of nitric oxide, a high-

enthalpy species, occurs due to the existence of high temperatures

during the combustion process. The process of formation is slow

relative to the time scale of other processes occurring in a typical

internal combustion engine, and so has to be considered on a rate

dependent basis and not on a chemical equilibrium basis. During the

expansion process, the concentration of nitric oxide is essentially

frozen at an early stage and hence the concentration at exhaust is

higher than that if chemical equilibrium were to prevail.

2.5 Effects of Various Parameters on Engine Emissions

2.5.1 Hydrocarbons

The hydrocarbon emissions depend primarily on the air-fuel

ratio and engine load conditions. The emissions decrease with an

increase in the air-fuel ratio to a certain extent until the mixture

is so lean that regular flame propagation is hindered. This is commonly

referred to as the 'misfire' condition. With further increase of the

air-fuel ratio, the hydrocarbon emissions increase. By increasing the

turbulence and having heated inlet manifolds, regular flame propagation

can be extended to very lean mixtures. In terms of the operating condi-

tions the influence of the air-fuel ratio may be summed up as: hydrocarbon

emissions are a maximum during engine warming, idle and deceleration.

However, since the amount of air intake under these conditions is small,

the absolute amount (gm./hr.) of emissions is not necessarily a maximum.

The emissions also increase with increasing load and compression ratio

and early spark timing. With greater spark retard, combustion is

completed later during the expansion stroke and is spread over a longer

time. The exhaust temperatures are also higher, thus promoting oxidation

in the exhaust system. The result is lower emissions, but the performance

may be adversely affected.

The emissions decrease with speed to a certain extent. At very

high speeds, the ignition voltage available at the spark plugs decreases

due to the short time available for coil buildup. This causes only

partial combustion, increasing the concentrations of unburnt hydro-

carbons. The recently introduced capacitive-discharge ignition

systems maintain a high enough voltage irrespective of the operating

speed. Increasing the average operating temperature and freeing the

combustion chamber from deposits serve to decrease the hydrocarbon


The quenching effect and hence the emission of hydrocarbons

are decreased by a decrease in the surface area to volume ratio

within the combustion chamber, so that a larger part of the mixture

is in the flammable zone. This is achieved by increasing the

compression ratio, increasing the stroke to bore ratio and modifying

the geometry to obtain compact combustion chambers avoiding dead

spaces and pockets. Increasing the combustion chamber surface

temperature also serves to decrease the extent of the quenching effects.

2.5.2 Carbon monoxide

Like hydrocarbons, the carbon monoxide emissions also depend

on the air-fuel ratio and engine load conditions. The carbon monoxide

emissions are a maximum during idling and deceleration when the air-

fuel ratio is towards the fuel-rich side. As with hydrocarbons, the

absolute amount (gm./hr.) of the emissions under these conditions

need not necessarily be at its maximum value. The carbon monoxide

emissions also increase with load and during acceleration periods.

The concentration is generally at a minimum during moderate power

cruising. Other factors which affect the hydrocarbon emissions

also affect the carbon monoxide emissions in a similar manner.

2.5.3 Oxides of nitrogen

A detailed discussion is provided in Chapter V; here only

brief mention is made of the effect of the important parameters.

The nitric oxide emissions depend primarily on three factors:

(i) The emissions are a maximum at an air-fuel ratio corresponding

to about 10% excess air (as compared to stoichiometric), and decrease

for richer or leaner mixtures.

(ii) Higher combustion temperatures favor larger concentrations

of nitric oxide.

(iii) The emissions are greater, the greater the time available

for the formation of nitric oxide, this being a strongly time-dependent


2.6 Emission Standards

The standards set of 1976 vehicles are (based on constant

volume sample cold start test average with a constant volume sample

hot start test, both with the Federal 22-minute driving cycle):

0.41 gm./vehicle-mile for hydrocarbons,

3.40 gm./vehicle-mile for carbon monoxide, and

3.00 gm./vehicle-mile for oxides of nitrogen.

Oxides of nitrogen have to be cut to 0.4 gm./vehicle-mile for 1977


For the sake of comparison, the emission levels without any

exhaust controls for a typical automobile are given below:

Hydrocarbons: 850 ppm (approximately 11 gm./vehicle-mile)

Carbon Monoxide: 3.4% by volume (approximately 80 gm./vehicle-mile)

Oxides of Nitrogen: 1000 ppm (approximately 4 gm./vehicle-mile)

With control techniques in current use the emission levels are

approximately (1973 standards based on constant volume sample cold

start test):

3.4 gm./vehicle-mile for hydrocarbons,

39.0 gm./vehicle-mile for carbon monoxide, and

3.0 gm./vehicle-mile for oxides of nitrogen.

[These figures would be somewhat different (lower for hydrocarbons and car-

bon monoxide) with the testing procedure employed for the 1976 standards.]

2.7 Control Methods

A brief mention of the various methods used to control the exhaust

emissions is made here.

(1) Leaner air-fuel ratios tend to decrease both the hydrocarbons and

carbon monoxide emissions; oxides of nitrogen are also reduced after a

certain air-fuel ratio. However, at these air-fuel ratios if regular

flame propagation cannot occur, the former two emissions may actually

increase rather than decrease. Presently, research is being conducted

to study the extension of the lean limit and also how to cope with the

generally adverse effects on driveability and some other operational

characteristics of the engine at extremely lean air-fuel ratios.

(2) Secondary air pumped either into the engine cylinders directly

during the expansion stroke or later in the exhaust manifold brings

about a reduction in the hydrocarbon and carbon monoxide emissions by

promoting more nearly complete oxidation. The nitric oxide emissions

are not significantly affected, they may increase slightly.

(3) Design modifications: These include improved carburetor performance,

a fast acting and more accurate choke, an inductive or electronic igni-

tion system with modified timing, reduced compression ratios, air

preheating, etc. Some of these call for a compromise between the

operational characteristics and the emission levels. Previously, the

best design was considered to be the one that resulted in the best

operational characteristics, and this is generally different from the

condition for lowest emission levels.

(4) Post-cylinder exhaust treatment of the emissions prior to their

exit from the exhaust system can be done in several ways besides

secondary air injection mentioned above. The most important of these

is the use of catalytic converters to promote the reactions in the

exhaust system in a direction which favors the oxidation of the unburnt

or partially burnt products, and also to decompose nitric oxide into

molecular nitrogen and oxygen, the two processes generally being carried

out in different stages. Extensive research has been done in this area

but still some problems remain before the system can be put in operation

on an extensive commercial basis. The high initial cost involved, the

interference of the lead in the gasoline with the catalyst action, and

the associated reduction of effective lifetime are some of the problems

to be considered.

Other modes of exhaust treatment include flame afterburners, but

these do not result in the reduction of nitric oxide emissions, and

besides require considerable auxiliary equipment.

(5) Exhaust gas recirculation serves to reduce the nitric oxide emissions,

but the carbon monoxide and hydrocarbon emissions are adversely affected

as are some of the operational characteristics of the engine.

There are several other control techniques which have been

suggested and the reader is referred to the appropriate literature

on the subject.

2.8 Description of the Combustion Process in an Otto Cycle Engine

In a spark ignition engine, the high temperature of the electric

spark causes the ignition of a small portion of the fuel-air mixture

in the immediate vicinity of the spark plug gap. The result is a more

or less spherical volume of hot gas consisting of burnt products which

are rapidly equilibrated due to the reactions being very fast at the

high temperature involved. This causes a nearly spherical flame to

originate from the point of ignition. Spherical flames are truly one-

dimensional systems by virtue of their geometry, whereas flat flames are

quasi-one-dimensional systems, there being small but negligible changes in

the various state properties in the direction perpendicular to that of

propagation. The fuel-air mixture in the combustion chamber is exposed

to the flame, layer by layer, until combustion is complete as far as

practicable. In practice, combustion is never totally complete due to the

quenching effects of the walls and other reasons. The entire combustion

process is extremely rapid, e.g., for combustion spread over 40 degrees

of crank revolution at 2000 revolutions per minute, the process would be

completed in approximately 3.33 milliseconds. In spite of the rapidity of the

process, at any given instant of time, the fraction of the total

fuel-air mixture actually being burnt is not significant if normal

combustion is occurring. This would not be the case when detonation

and knocking occur. Therefore, under normal circumstances, the

combustion process occurs with uniform release of energy and may be

approximated in the initial stages by the model of a one-dimensional

laminar spherical flame propagating into a premixed fuel-air mixture.

The assumption of a laminar flame is made to render the problem

solvable (in general, turbulence would be present in varying degrees),

there being no convenient numerical techniques to handle a turbulent

flame for a general case. As the flame propagates, it would deviate

from its spherical shape, the deviation depending on the geometry of

the combustion chamber walls with respect to the point of ignition.

This deviation can be accounted for by considering the general three-

dimensional flame equations in their time-dependent form, and a solution

so obtained would give the flame shape and position as a function of

time. Even with modern electronic computers this problem is almost

insurmountable, and would hardly give more insight to the actual

combustion process than the solution of a one-dimensional problem.

2.9 Laminar Flames: Introductory Remarks

Laminar premixed flames represent a unique case of wave

propagation--the propagation of a self-sustaining exothermic reaction

wave into a mixture of combustible gases. Heat conduction, diffusion,

and a set of rapid exothermic chemical reactions all are significant

in the propagation of such a flame. Thus this is one of the earliest

combustion problems that has necessitated the simultaneous consideration

of both the fluid dynamic and chemical kinetics aspects for its solution.

Laminar flame velocity (also referred to variously as the burning

velocity or the normal combustion velocity), which is of primary

importance in the study of premixed flames, is defined as the velocity

at which unburnt gases move through a combustion wave in the direction

normal to the wave surface. This definition refers to the velocity

of unburnt gas relative to the flame front as determined at some point

far enough removed from the flame, such that the influence of the flame

itself at the point is negligible.

The propagation characteristics of steady, adiabatic, one-

dimensional laminar premixed flames are governed primarily by two

properties of the system--the activation energy of the rate controlling

step and the Lewis number, where the Lewis number represents the ratio

of the energy transported by conduction to that transported by diffusion.

The geometry is of secondary importance, and therefore, within reason-

able limits, the laminar flame velocity may be defined independent of

a particular experimental apparatus. It will be seen later that the

flame velocity is an eigenvalue of the conserv- ,.-- equations.

A laminar flame is made up of a series of fairly well-defined

regions. The thickness of the flame zone is rather arbitrary, similar

to the thickness of a boundary layer. In the first region, the unburnt

gas undergoes appreciable changes in composition and temperature due to

the transport processes of diffusion and thermal conduction, but the

reaction rates are negligible. This is referred to as the preheat

zone. After this there is a relatively thin primary reaction zone

wherein all the important rapid chemical reactions occur. This is the

luminous zone in hydrocarbon flames. In addition to this, there may

be one or more secondary reaction zones whose spatial dimensions are

an order of magnitude larger than those of the primary reaction zone,

and wherein slow reactions (e.g., carbon monoxide afterburning in

hydrocarbon flames) and radical recombination reactions occur. At

low enough pressures the radical recombination reactions and the

afterburning may occur in distinctly separate zones; on the other hand,

with increasing pressure, the secondary reaction zone may partly merge

into the primary zone. The existence of these zones and carbon monoxide

afterburning has been confirmed by experiments.

2.10 Closing Remarks

An introduction to the study of the combustion process in and

the emissions from an Otto cycle engine has been presented. These

ideas will be developed further in the succeeding chapters.

Chapter III presents a review of the extensive literature on

the methods of solution of the laminar flame propagation equations.

Chapter IV presents a mathematical analysis of the combustion

process--the formulation and solution of the laminar flame propagation


In Chapter V a mathematical model has been developed to obtain

a detailed thermodynamic analysis of the Otto cycle.

Results and discussions are presented in Chapter VI and finally

Chapter VII is used for some closing remarks.


The problem of laminar flame propagation in premixed gases has

been studied fairly extensively. The problem was first studied by

Mallard and LeChatelier (1883). Their theory has been referred to as

the thermal theory in the sense that they considered thermal effects to

be of primary importance and the chemical kinetics to be secondary.

Later on, Taffanel (1913) and Daniell (1930) demonstrated for a highly

simplified model that the flame propagation velocity is proportional

to the square root of the ratio of the thermal conductivity to the

specific heat at constant pressure and to the square root of the rate

of the reaction.

An excellent review of the early literature (up to about 1950)

on the steady state flame propagation in premixed gases has been

presented by Evans [1]. Williams [2, p. 95] has presented a review

with references to more recent literature (up to about 1960).

The assumption of Lewis number equal to unity (normal diffusion)

has been widely made in flame analysis. This assumption implies that

the energy fluxes due to diffusion and conduction balance each other,

and so the enthalpy is constant throughout the flame zone; hence the

enthalpy equation need not be integrated. This simplifies the problem

considerably and even permits closed form solutions for a few simple

cases. The assumption of unity Lewis number is, however, not valid for

many gases encountered in combustion problems, two notable examples

being hydrogen-air and propane-air mixtures (Le = 3.30, and Le = 0.49,


Assuming that the activation energy of the rate controlling

step and the heat of reaction are large and that all the reaction occurs

at the hot boundary, Zeldovich and Frank-Kamenetski [3] were able to

obtain a closed form solution for the burning rate eigenvalue. Semenov

[Sokolik, 4] has modified their solution for non-normal diffusion by

multiplying the flame velocity as calculated by the Zeldovich-Frank-

Kamenetski formulation by the factor 1/vLe for first order reactions

and by 1/Le for second order reactions. This analysis is also limited

by the requirement of large or infinite activation energies. Von Karman's

zeroth and first order approximations [5] also yield closed form solu-

tions, the former being a special case of the Zeldovich-Frank-Kamenetski

analysis wherein only the first term of a series expansion for the

eigenvalue is retained. Some of the other analytical solutions are

those of Boys and Corner [6, 7] and Adams [8]. Wilde [9] empirically

modified Adams' solution thereby obtaining better agreement with 'exact'

numerical solutions. Hirschfelder [10] has obtained two approximate

solutions in one of which the eigenvalue is obtained as a root of a

quadratic equation and in the other as a root of a third order equation.

The accuracy is comparable to that obtained with Wilde's solution.

Williams [2] has given a graphical comparison of the accuracy of these

approximate methods.

Spalding [11] has proposed a centroidd rule' which provides a

convenient way of estimating the flame speed to an accuracy of about 1%

for the case in which the reaction rate is an explicit function of

temperature. This is true for unity Lewis number of the reactants

and radicals for the following class of reactions: (a) single step

reactions, (b) reactions satisfying the steady state approximation,

and (c) explosive chain reactions. The method gives an expression

for the flame speed in terms of the centroid of the reaction rate

function without resorting to the solution of the governing dif-

ferential equations. Graphs are provided for estimating the flame

speed for non-unity Lewis number. These show that when the Lewis

number of the radicals is greater than unity the flame speed in-

creases, but only slightly. Small values of Lewis number (of the

radicals), on the other hand, greatly reduce the flame speed. For

the reactants, a Lewis number greater than unity decreases the

flame speed.

Rosen [12, 13] has suggested a variational method employing

the Rayleigh-Ritz method for obtaining successively better estimates

of the functional and eigenvalue. He has also given a plausible

derivation of the centroid rule proposed by Spalding [11]. The

centruid rule method has been extended for the case of non-normal

diffLusion with bimolecular reactions (A+B-C) in a stoichiometric

two-component system by Adler [14]. The flame speed has been shown to

be inversely proportional to the Lewis number of the reactants for

large activation energies.

Spalding [15] has compared several approximate solutions for

the case of unity Lewis number and temperature explicit reactions and

arranged them in order of merit as: Wilde (best), Adams, Von Karman-

Penner, Zeldovich-Frank-Kamenetski, and Corner. Wilde's method has

been extended to study the effect of the diffusion coefficient of

reactants and radicals on the flame speed and fairly good agreement

with 'exact' solutions is reported.

Adler [16] has considered temperature explicit reaction rate

functions and obtained limits on the burning rate eigenvalue independent

of the type of reaction considered for a fixed centroid of the reaction

rate function. For the limiting case of high activation energy the

upper and lower limits are close, so that the flame speed may be

obtained with good accuracy.

Favin and Westenberg [17] have considered one-dimensional

steady premixed spherical laminar flames and solved the equations for

the case of unity Lewis number and a unimolecular reaction process.

They have integrated the equations from the cold to the hot boundary

and considered the flameholder temperature, rather than the mass flow

rate, as the eigenvalue.

Adler and Spalding [18] have obtained analytical solutions for

the temperature gradient, radical concentrations, and the flame speed

for the case of non-branching chain reactions with negligible chain


De Sendagorta [19] has obtained a numerically exact solution

for one value of the reduced activation temperature E/RTb and used

this to show that good accuracy could be obtained by using an

exponential approximation for the species flux fraction in analytical

solutions. His solution covers both first and second order single

step reactions and is considerably more accurate than the other

analytical solutions cited above. Non-normal diffusion has been

considered by Wilde [9], de Sendagorta [19] and Von Karman-Penner

[20] utilizing Karman-Pohlhausen profile methods.

Millan and da Riva [21] have obtained a solution for the

distribution of radicals in laminar flames based on the steady state

approximation. A general criterion for the applicability of this

approximation is given and it is shown that this assumption is a first

order approximation to the actual solution, and this solution can be

improved by a method of successive approximations.

Spalding and Adler [22, 23] have considered one-dimensional

laminar flames with an enthalpy gradient (defined as positive for

heat loss at the cold boundary) which is useful when heat losses by

conduction or radiation are significant, and obtained numerical

solutions. For unity Lewis number, they have shown that the flame

speed decreases or increases as the enthalpy gradient decreases or

increases from zero, the effect being more, the higher the activation

energy. The flame speed decreases to zero for a certain critical

negative value of the enthalpy gradient.

Heat losses and flame propagation limits due to heat losses

have also been considered in the analyses presented by Mayer [24],

Spalding [25], Adler [26], and Idler and Kennerley [27]. From these

analyses it can be concluded that for a given amount of heat loss,

a given combustible mixture has two possible flame speeds, the lower

one of which is normally unstable. With increasing heat loss, the

two speeds approach each other and finally coincide for a critical

value of the heat loss. For still greater heat loss, the flame speed

is imaginary. The critical value of heat loss is identified with the

flammability limit. Some other studies of non-adiabatic flames are

given in References 28, 29, and 30.

With the advent of modern electronic computing machines, the

emphasis has gradually shifted to obtain more meaningful solutions with

less drastic assumptions, by using numerical methods. These studies

involve less stringent assumptions, are more accurate, and include

better estimates of thermodynamic and transport properties, and also

consider the effects of chain reactions instead of the classical

unimolecular process R P (R = reactant, P = product), considered

in most of the previous analytical solutions. Firstly, brief mention

is made of the early numerical solutions.

Von Karman's second approximation [5] involves an iterative

procedure with numerical integration. The approximate solution of

Boys and Corner [6] can be improved by numerical iterations. Some

of the other iterative procedures requiring numerical integration

are those due to Klein [31] and Johnson and Nachbar [32]. The latter

one is the most accurate of the simple approximate procedures and it

permits the determination of the upper and lower bounds for the

eigenvalue. Using their iteration scheme successive upper and lower

bounds can be obtained until the eigenvalue converges to any desired

accuracy. Generally the first approximation is accurate enough for

practical purposes so that it is seldom necessary to use the iterative


Friedman and Burke [33] have obtained numerical solutions for

the case of a hypothetical first order irreversible reaction with

exponential dependence of rate on temperature over a wide range of two

parameters: dimensionless activation temperature (E/RTb) and the

Lewis number.

The theoretical analysis presented by Giddings and Hirschfelder

[34] was the first one to consider branched chain reactions. The

assumption of a steady state condition was made to approximate the

radical concentrations and the governing equations were solved by the

integral equation approach devised by Klein [31].

Menkes [35] has solved the problem for a cylindrical flame

by transforming the differential equation into an integral equation

which is solved numerically by a method of successive approximations.

The cylindrical flame was investigated to get an insight into the

stability analysis. It was found that the flame radius decreases

and the specific mass flow increases rapidly with increasing inlet

temperature and at constant flame thickness. However, if the flame

radius is decreased, keeping the inlet temperature constant, the

flame thickness increases rapidly resulting in a decrease of the

specific mass flow rate. Similar results have been obtained by

Menkes [36] by an approximate analytical solution. Jain and Ebenezer

[37] have obtained analytical solutions for a cylindrical flame, wherein

the reaction rate is represented by a step function, and have studied

the accuracy of the thin-flame" approximation.

A thin flame is one whose thickness is small compared to its
radius of curvature, so that the mass rate of burning per unit area
can be taken as equal to that in a plane flame.

Spalding [38] has obtained numerical solutions for the case

of a steady laminar spherical flame by stepwise integration from the

cold boundary. The model used consists of a combustible gas emerging

from a point source into an atmosphere of adiabatic combustion products.

The results have been compared with those obtained from the thin-flame

theory and the two are shown to be asymptotically equal for large mass

flow rates and high reaction rates. For the same model Spalding and

Jain [39] have obtained analytical solutions for the cases of (a) step

function reaction rate curves and (b) Adams-type reaction rate curves

based on the assumption of a constant temperature gradient in the

reaction zone. It is shown that the radius of the flame depends more

on the position of the centroid of the reaction rate curve rather

than on its specific shape. The same results have been confirmed by

means of a resistance network analog computer by Spalding, Jain and

Samain [40].

Jain and Kumar [41] have employed the method of matched

asymptotic expansions and obtained analytical solutions for the

limiting case of large activation temperatures. They have also

obtained numerical solutions for a wide range of values of both the

reduced activation temperature and Lewis number, for both first and

second order reactions. Based on the numerical solutions, a reason-

ably accurate rule for the flame speed eigenvalue in terms of Lewis

number and the centroid of the reaction rate function is suggested;

it is linear in Lewis number for first order reactions and quadratic

in Lewis number for second order reactions. There is good agreement

of these results with those of Spalding [11, 15] and de Sendagorta [19].

However, the present solution agrees with Adler's [14] solution only

for Lewis numbers around unity, and there is considerable disagreement

when the Lewis number differs significantly from unity for which a

possible explanation has been given.

Besides the various methods mentioned above, other recent

methods have sought not just to determine the laminar flame speed

but also the detailed flame structure regarding the distribution of

temperature and the various chemical species over the thickness of

the flame zone. A detailed analysis of the flame structure gives

a better insight into the mechanism of flame propagation and also

permits the study of chemical reactions from the viewpoint of

mechanism and kinetics. The difficulties involved in such a detailed

analysis of the flame zone may be summed up into three categories:

(a) a reliable postulation of the reaction mechanism for the oxidation

of the fuel; (b) a reliable estimate of the transport properties,

especially at temperatures encountered in flames and for free

radicals; and (c) the absence of suitable numerical techniques to

handle the complicated system of equations in a convenient way.

Basically two mathematical formulations have been widely used

in the detailed analysis of the flame zone. The first one which is

a trial and error method is due to Hirschfelder [42]. This method

consists of setting up the differential equations for the conservation

of mass and energy and the diffusion equations in a coordinate system

which is at rest with respect to the flame front. The boundary conditions

are specified at the burnt and unburnt ends. In practice, integrating

these equations to obtain the concentrations of species and state

properties in the flame zone requires a knowledge of the flame speed.

Since this is not known a priori, the problem reduces to an eigenvalue

problem; a solution to the equations exists for only one unique value

of the flame speed, known as the eigenvalue of the conservation

equations or the eigenvalue of the two-point boundary value problem.

Assuming various flame speeds, the equations are numerically integrated

from the burnt to the unburnt end. The solution which agrees best with

the boundary conditions is the one assumed to be correct. Convergence

criteria for this method are not available. For more complicated cases,

where the chemical kinetics is governed by a set of reaction steps, the

assumption of flame speed alone is not sufficient to start the integra-

tion of the conservation equations, and one is confronted with a multi-

eigenvalue problem in non-separable equations [1, 43, 44]. The number

of eigenvalues is equal to the number of linearly independent species

flux fractions in the flame equations. This multi-eigenvalue problem

has, however, been reduced to the determination of a single eigenvalue,

without loss of generality, in a method introduced by Campbell, Heinen

and Schalit [44]. They have obtained numerical solutions for a one

dimensional flame with idealized kinetics and transport properties.

Three examples using this trial and error method have been given

by Hirschfelder [42]: (a) the unimolecular decomposition of hydrazine,

(b) the bimolecular decomposition of nitric oxide, and (c) the two step

chain mechanism describing the decomposition of ozone. Agreement with

experimental results is fair* for the first and third cases; no

It should be noted that the solution method itself, even if
extremely powerful, would not give results in good agreement with the
experimental data if the accuracy of the reaction kinetics and transport
property data is poor. Thus the judgment of accuracy should rather be
based on the solution of the same problem with the same input data by
different methods and comparison with experiments.

experimental data were available for comparison in the second case.

The second method is due to Spalding [45] and consists of

setting up the differential equations for conservation of mass and

energy and the diffusion equations in their time-dependent or unsteady

form. Assuming arbitrary initial profiles for the state properties

and the species concentrations, the equations are solved by forward

marching finite difference methods until the various profiles of the

species concentrations and the state properties converge to a steady

state. A big advantage of this method is that almost any arbitrary

initial profile can be used for starting the solution; however, by

using profiles which closely resemble their steady state shape, the

computation can be cut down considerably. For the temperature and

the major species, S-shaped profiles asymptotically approaching the

boundary values are commonly used. The radicals are usually assumed

to be at zero concentration initially. The assumed profiles will-

develop into the steady case if enough energy is stored in them.

If this initial energy is insufficient then the gradients of temper-

ature and species concentrations will decrease monotonically and the

flame is extinguished. Thus this method can also be used to obtain

an estimation of the minimum ignition energy, if required. No

iterations are required with this approach. The computational

advantage or disadvantage of either method has not been conclusively

established as yet.

This method has been demonstrated by Spalding for the case of

the hydrazine decomposition flame. The conservation equations reduced to

two partial differential equations and these were solved simultaneously

by Schmidt's graphical method for the solution of the unsteady heat

conduction equation. Calculations were carried out for two different

flame temperatures and fairly good agreement between theoretical and

experimental values were obtained.

Spalding's method has been used by various researchers for

solving different flames. The later techniques used are basically

refinements of the finite difference solutions and the inclusion of

chain-branching reactions and other improvements. A review of some

of the recent work using this method is presented below.

Zeldovich and Barenblatt [46] have reviewed previous work

on the steady state propagation of flat flames. They have also

considered the effects of the Lewis number, heat loss and chain

reactions on the velocity of propagation of flames using Spalding's

time-dependent approach. Using a perturbation technique the stability

of the steady state solution has been examined. They have shown that

for non-unity Lewis number, the total enthalpy reaches a maximum or

a minimum in the flame zone, depending on whether the Lewis number is

greater or smaller than unity. The same result was also obtained

earlier by Lewis and Von Elbe [47]. Heat loss was found to decrease

the velocity of propagation and there was found to exist a certain

critical value of the heat loss parameter 'b' [q = b(t-t )], above

which a steady state solution did not exist.

Adams and Cook [48] have used the time-dependent method to

study the hydrazine decomposition flame which has been fairly

extensively studied. As before, the problem was reduced to the

simultaneous solution of two partial differential equations: one

involving the temperature and the other the mass fraction of hydrazine.

A central difference formula was used for approximating the space

derivatives and a forward difference formula for the time derivatives.

The equations were then numerically solved by explicit forward marching

steps in time until the profiles of temperature and mass fraction reached

a steady state. Considering chain breaking reactions of both second and

third order, it was shown that the dependence of flame speed on temper-

ature is insensitive to changes in pressure and the mechanism of the

chain breaking reactions, while the dependence on pressure is sensitive

to the mechanism assumed.

Dixon-Lewis [49] has implemented Spalding's time-dependent

method for the hydrogen-oxygen flame. To reduce the computational

work, some simplifications regarding the transport and thermodynamic

properties were made. The chemical kinetics was treated in terms of

a set of elementary reaction steps and thus this analysis avoids some

of the misleading results of previous analyses which considered the

chemical kinetics to be governed by a single reaction or a simple

unrealistic mechanism. Important conclusions regarding the influence

of some parameters on the hydrogen-oxygen flame have been made in

this analysis. For the sake of computational stability, the standard
2 1
stability condition of At/(A)2 < 1 for the heat conduction equation

had to be modified to At/(AI) 2 1/30 for this problem.

Detailed transport property calculations have been made [50]

and considered in the analysis of the hydrogen-oxygen-nitrogen flame

by Dixon-Lewis [51] using the Runge-Kutta-llerson procedure (Runge-Kutta

method with variable step size).

Spalding et al. [52] have utilized the mathematical similarity

of the equations governing steady two-dimensional boundary layer flows

with the one-dimensional unsteady flame equations. The computational

procedure developed by Patankar and Spalding [53] for the study of the

hydrodynamical and heat transfer aspects of steady two-dimensional

boundary layer flows has been adapted to the study of the one-dimensional

unsteady hydrazine decomposition flame. Both these phenomena are governed

by sets of simultaneous parabolic partial differential equations. The

assumptions of unity Lewis number and constant and uniform specific

heats were made in the solution method, although neither was necessary

for the computational procedure. Initial concentration profiles of the

form yA = Au (10w 15w + 6w ) were used for the major species, and

the radicals were assumed to be initially at zero concentration. Here,

w = (-' )/(u b), where k = stream function, and the subscripts b

and u refer to the burnt and unburnt ends, respectively. A variable

step size for the space dimension was used and the equations were solved

by an implicit method after linearization at each step. In the absence

of formal stability criteria, the largest time step size was found by

trial and error.

Wilde [54] has treated this problem as a boundary value problem

and formulated it in three different ways: (a) a set of first order

time-independent ordinary differential equations for species conserva-

tion, diffusion and thermal flux; (b) a set of second order time-

independent ordinary differential equations; and (c) a set of time-

dependent parabolic partial differential equations. The first for-

mulation was solved by quasi-linearization, also known as Newton-

Raphson-Kantorovich linearization, and then integrating by using

algorithms suitable for stiff equations. The second formulation

was also linearized and then the equations were written in finite

difference form, the result being a set of linear algebraic equation,

block tridiagonal in form. For the last formulation, explicit finite

difference methods were found to require a very small time step size

due to the presence of the term involving the rate of production of

species by chemical reactions. Considerable time saving was obtained

by using forward difference formulas in time, implicit representation

of the space derivatives and an averaging of the chemical rate term

at the present and future time step. The second formulation was

faster per iteration than the first one, but required closer initial

estimates and a greater number of iterations for convergence. The

last method was longer in running time than the other two but was

surer in convergence. Three flames were studied using this method:

ozone decomposition, hydrogen-bromine flame, and the hydrogen-oxygen

flame. The calculated flame velocities were about twice the cor-

responding experimental values. The discrepancy was explained by

blaming the unreliable kinetic and transport data.

Bledjian [55] has provided a general formulation which allows

for the consideration of spherical and cylindrical flames in addition

to plane flames. An important difference between plane and non-planar

flames should be mentioned here. For non-planar flames the spatial

variable (radius) appears explicitly in the gradient terms while for

plane flames the spatial variables appear only in the derivatives.

Thus for the simple cases, the space variable can be eliminated and

the governing equations can be reduced to a single equation with

temperature as the independent variable and the mass flux as the

dependent variable. An added advantage of this is that the range

of integration is now finite in the independent variable.

The method of lines was used to solve the governing equations

in Bledjian's formulation [55]. This consists of converting the set

of partial differential equations to a set of ordinary differential

equations by discretizing all but one of the independent variables

at a time, the resulting equations being then integrated using

Runge-Kutta's fourth order method. Difficulties were found when the

set of equations was stiff. The step size for the next integration

cycle was estimated based on the local truncation error, as calculated

by using Richardson's extrapolation formulas. Calculations were carried

out for two flames: the hydrazine decomposition flame and the ozone

decomposition flame. The results for the first case agreed well with

Spalding's calculations [45] and for the second case good agreement

was obtained with experimental results. Estimations of the minimum

ignition energy for the two cases have also been obtained. For the

sake of computational stability, the time step size was restricted by

the condition (X/pC )At/(A,)2 1/30 instead of the standard stability

condition (X/pC )At(ALp)2 1/2 for the heat conduction equation.


1. Evans, M. W., "Current theoretical concepts of steady-state
flame propagation," Chem. Rev., 51, 363 (1952).

2. Williams, F. A., "Combustion Theory," Addison-Wesley (1965).

3. Zeldovich, Y. B. and Frank-Kamenetski, D. A., "Theory of uniform
flame propagation," J. Phys. Chem. (U.S.S.R.), 12, 100 (1938).

4. Sokolik, A. S., "Self-ignition, flame and detonation in gases"
(Translated from Russian), Jerusalem (1963).

5. Von Karman, T., "The present status of the theory of laminar
flame propagation," Sixth (Int'l) Symp. on Combustion, Reinhold
Publishing Corp., New York (1956).

6. Boys, S. F. and Corner, J., "The structure of the reaction zone
in a flame," Proc. Roy. Soc. London, A197, 90 (1949).

7. Corner, J., "The effect of diffusion of the main reactants on
flame speeds in gases," Proc. Roy. Soc. London, A198, 388 (1949).

8. Adams, E. N., Rept. CF-957, Univ. of Wisconsin Naval Research
Lab. (1948).

9. Wilde, K. A., "Improved approximate solutions of flame equations
governed by simple chemical reactions," J. Chem. Phys., 22, 1788

10. Hirschfelder, J. O., "Theory of flames produced by unimolecular
reactions," Phys. Fluids, 2, 565 (1959).

11. Spalding, D. B., "I. Predicting the laminar flame speed in gases
with temperature explicit reaction rates," Combustion and Flame,
1, 287 (1957).

12. Rosen, G., "An action principle for the laminar flame," Seventh
(Int'l) Symp. on Combustion, Butterworths Scientific Publications,
London (1959).

13. Rosen, G., "Generalization of the laminar flame action principle
for Arrhenius-type rate functions," J. Chem. Phys., 32, 311

14. Adler, J., "The prediction of laminar flame speeds in stoichiometric
mixtures with non-normal diffusion," Combustion and Flame, 9, 273

15. Spalding, D. B., "II. One-dimensional laminar flame theory for
temperature explicit reaction rates," Combustion and Flame,
1, 296 (1957).

16. Adler, J., "The limits of the eigenvalue of the laminar flame
equation in terms of the reaction rate-temperature centroid,"
Combustion and Flame, 3, 389 (1959).

17. Favin, S. and Westenberg, A. A., "The theory of a spherical
premixed flame," Combustion and Flame, 4, 161 (1960).

18. Spalding, D. B. and Adler, J., "Flame propagation by non-
branching chain reactions," Combustion and Flame, 5, 123

19. De Sendagorta, J. M., "Propagation velocity and structure of
single step reaction flames," Combustion and Flame, 5, 305

20. Von Karman, T. and Penner, S. S., "Selected Combustion Problems,"
AGARD, Vol. I, Butterworths, London (1954).

21. Millan, G. and Da Riva, I., "Distribution of radicals in laminar
flames," Eighth (Int'l) Symp. on Combustion, Williams and Wilkins,
Baltimore (1962).

22. Spalding, D. B. and Adler, J., "One-dimensional laminar flame
propagation with an enthalpy gradient," Combustion and Flame,
4, 94 (1960).

23. Adler, J. and Spalding, D. B., "One-dimensional laminar flame
propagation with an enthalpy gradient," Proc. Roy. Soc. London,
A261, 53 (1961).

24. Mayer, E., "A theory of flame propagation limits due to heat
loss," Combustion and Flame, 1, 438 (1957).

25. Spalding, D. B., "A theory of inflammability limits and flame
quenching," Proc. Roy. Soc. London, A240, 83 (1957).

26. Adler, J., "One-dimensional laminar flame propagation with
distributed heat losses: Thin-flame theory," Combustion and
Flame, 7, 39 (1963).

27. Adler, J. and Kennerley, J. A., "Steady laminar propagation
with conductive heat loss," Combustion and Flame, 10, 191

28. Chen, T. and Toong, T., "Structure and propagation of laminar
flames near a heat sink," Combustion and Flame, 4, 313 (1960).

29. Yang, C. H., "Burning velocity and structure of flames near
extinction limits," Combustion and Flame, 5, 163 (1961).

30. Adler, J., "Downstream boundary conditions of non-adiabatic
one-dimensional laminar flames," Combustion and Flame, 11,
85 (1967).

31. Klein, G., "A contribution to flame theory," Phil. Trans. Roy.
Soc. London, A249, 389 (1957).

32. Johnson, W. E. and Nachbar, W., "Deflagration limits in the
steady linear burning of a monopropellant with application to
ammonium perchlorate," Eighth (Int'l) Symp. on Combustion,
Williams and Wilkins, Baltimore (1962).

33. Friedman, R. and Burke, E., "A theoretical model of a gaseous
combustion wave governed by a first order reaction," J. Chem.
Phys., 21, 710 (1953).

34. Giddings, J. C. and Hirschfelder, J. 0., "Flame properties and
the kinetics of chain-branching reactions," Sixth (Int'l) Symp.
on Combustion, Reinhold Publishing Corp., New York (1956).

35. Menkes, J., "A contribution to the theory of laminar flames
with radial symmetry," Combustion and Flame, 4, 1 (1960).

36. Menkes, J., "Analytical solution of a flame with cylindrical
symmetry," Eighth (Int'l) Symp. on Combustion, Williams and
Wilkins, Baltimore (1962).

37. Jain, V. K. and Ebenezer, J. S., "Analytical solutions of a
steady cylindrical flame," Combustion and Flame, 10, 391 (1966).

38. Spalding, D. B., "The theory of steady laminar spherical flame
propagation: equations and numerical solution," Combustion
and Flame, 4, 51 (1960).

39. Spalding, D. B. and Jain, V. K., "The theory of steady laminar
spherical flame propagation: analytical solutions," Combustion
and Flame, 5, 11 (1961).

40. Spalding, D. B., Jain, V. K., and Samain, M. D., "The theory of
steady laminar spherical flame propagation: analog solutions,"
Combustion and Flame, 5, 19 (1961).

41. Jain, V. K. and Kumar, R. N., "Thegry of laminar flame propagation
with non-normal diffusion," Combustion and Flame, 13, 285 (1969).

42. Hirschfelder, J. O., Curtiss, C. F., and Campbell, D. E., "The
theory of laminar flame propagation IV," J. Phys. Chem., 57,
403 (1953).

43. Campbell, E. S., "Theoretical study of the hydrogen-bromine
flame," Sixth (Int'l) Symp. on Combustion, Reinhold Publishing
Corp., New York (1956).

44. Campbell, E. S., Heinen, F. J., and Schalit, L. M., "A theoretical
study of some properties of laminar steady state flames as a
function of properties of their chemical components," Ninth (Int'l)
Symp. on Combustion, Academic Press, New York (1963).

45. Spalding, D. B., "The theory of flame phenomena with a chain
reaction," Phil. Trans. Roy. Soc. London, A249, 1 (1956).

46. Zeldovich, Y. B. and Barenblatt, G. I., "Theory of flame
propagation," Combustion and Flame, 3, 61 (1959).

47. Lewis, B. and Von Elbe, G., "Combustion, flames and explosions
in gases," 2nd Ed., Academic Press, New York (1961).

48. Adams, G. K. and Cook, G. B., "The effect of pressure on the
mechanism and speed of the hydrazine decomposition flame,"
Combustion and Flame, 4, 9 (1960).

49. Dixon-Lewis, G., "Flame structure and flame reaction kinetics
I. Solutions of conservation equations and applications to
rich hydrogen-oxygen flames," Proc. Roy. Soc. London, A298,
495 (1967).

50. Dixon-Lewis, G., "Flame structure and flame reaction kinetics
II. Transport phenomena in multicomponent systems," Proc.
Roy. Soc. London, A307, 111 (1968).

51. Dixon-Lewis, G., "Flame structure and flame reaction kinetics
V. Investigation of reaction mechanism in a rich hydrogen+
nitrogen+oxygen flame by solution of conservation equations,"
Proc. Roy. Soc. London, A317, 235 (1970).

52. Spalding, D. B., Stephenson, P. L., and Taylor, R. G., A
calculation procedure for the prediction of laminar flame
speeds," Combustion and Flame, 17, 55 (1971).

53. Patankar, S. V. and Spalding, D. B., "Heat and mass transfer in
boundary layers," 2nd Ed., Intertext Books, London (1970).

54. Wilde, K. A., "Boundary-value solutions of the one-dimensional
laminar flame propagation equations," Combustion and Flame, 18,
43 (1972).

55. Bledjian, L., "Computation of time-dependent laminar flame
structure," Combustion and Flame, 20, 5 (1973).


4.1 Formulation of the Problem

In this chapter the case of laminar premixed flame propagation

in a stoichiometric mixture of methane and air, initially at a pressure

of 10 atmospheres and at a temperature of 628 degrees K, which approx-

imately is the condition in an Otto cycle engine at the time of ignition,

has been considered. The time-dependent approach of Spalding [1] has

been used here. The reason methane has been considered is because of

the availability of the reaction kinetics data. The method is general,

and if the reaction mechanism and kinetic data for the oxidation of any

other fuel are known, the same method can be used to obtain the detailed

flame structure. The variations of the transport and thermodynamic

properties with temperature, pressure and chemical composition have

been accounted for as far as the available data permit.

4.1.1 The conservation equations

The flame zone is generally a few hundred mean free paths

thick, and so the concepts of continuum analysis may be used here.

The general conservation equations for multi-component reacting

ideal gas mixtures are given in various references [2,3]. These

equations are;

(1) Overall continuity:

-t + V. () = 0 (1)



+ (vV)v = (V-p)/p + g Y.f.


aU -- -> -V
p + p*vVu = Vqj p: (V) + p Z Y .iV.

Continuity of species:

5Y. w.
i i -
-+ vVY. [V(pYV)]/ ,
1t i p ( )]

i=1,2,..... N

In addition to the above conservation equations, other equations

describing the fluxes of mass and energy, the chemical kinetics of the

system and some constitutive relations are necessary. These are (again

in their general form):

(1) Diffusion in multi-component mixtures is governed by:

N X.X-ij ij +
VX. .(V V.) + (Y.- X.)E + i j i j
1 j=l Di. j jD i 1 1 p j A =1

ff (ff D% (VT]
P- Di Y Y
J=l j 3

i=1,2,......... .N

Y.V = 0

The derivation of equation (5) assumes that the medium under

consideration can be treated as a dilute gas. However, there is


experimental evidence that this result also applies quite well in

dense gases and liquids [3, p. 718].

(2) Heat flux vector:

N N N rX.D
q = VT h p h.Y.V. + R T (Vi V.) (7)
I I i u j=l WiDij

The same comment as above applies to this equation also.

(3) Chemical kinetics for a set of M elementary reactions between

N species:
M R N tXpn Jk N X p k
k=l ( 2 j=l kN kbk N l u '
k=1 u (8)

(4) Thermally perfect gas:

P = pR T/W (9)

S= X X.W. (10)

(5) Caloric equatiai of state:

h. = h + C .dT i=1,2,........N (11)
S 1i p, '

Implicit in these equations is the assumption that every point

in the flow field has a thermodynamically definable state, i.e., it is

possible to attach significance to the concepts of local temperature,

pressure, density and chemical composition (see Sec. 4.4 for un-

certainties in the definition of temperature).

4.1.2 Major assumptions for analysis

The general equations are now simplified for the case of one-

dimensional laminar flame propagation in spherical coordinates using

the following assumptions.

(1) The kinetic energy associated with the motion of the gases is

small compared to the heat of the reaction and may be neglected. Also,

since the velocities and velocity gradients are small, viscous effects

may be neglected. Thus the momentum equation can be replaced by the

assumption that the pressure is constant through the flame. In fact,

solving the momentum equation for the pressure change across the

flame yields:

Ap = P, P = p M ; M = (12)
u u

From experimental observations it is seen that in general

M < 0.02, V /V < 15, and y < 1.667, which yields Ap 0.01 p .
u b u u
Hence for ordinary flames the process may be considered essentially

isobaric. This assumption would not hold for detonation waves.

(2) Radiative heat transfer is negligible. Although relatively

high temperatures (- 25000K) are encountered in flames, radiation is

not important since the gases have low emissivity and absorptivity and

so the radiative heat transfer during the reaction time is small compared

to the total energy flux. With flames containing solid carbon

this is generally not the case.

(3) The Soret and Dufour effects are negligible.

(4) Any diffusion caused by pressure gradients is negligible. Since

assumption (1) already considers the process to be isobaric, this

particular assumption is not necessary. This term may be significant

in strong detonation waves wherein tremendous pressure gradients

may be established.

(5) External body forces are neglected. For very hot flames wherein

ionization occurs to a significant degree, this term may have to be

considered if an external electric field is also present; the external

force on an ion is equal to the product of the ionic charge and the

local electric field strength. This gives rise to what is termed

as 'forced diffusion'.

(6) The chemical kinetics is assumed to follow Arrhenius equation.


kk A= kTb exp(-E /RT) (13)

The backward reaction rate constant kbk is then obtained from

a knowledge of kfk and the equilibrium constant KCk:

K k (14)
Ck kbk

There is considerable controversy [4,5] regarding the validity

of equation (14); it will nevertheless be used in the absence of a

more accurate relation.

(7) The gas mixture is thermally perfect, i.e., it obeys the ideal

gas equation. Nitrogen is a major constituent (> 70%) of the

mixture and it closely approximates thermally perfect behavior for

the conditions considered in this study. Other constituents are

present in relatively smaller concentrations, so as not to affect

the overall thermally perfect behavior.

The equations finally simplify to:

(1) Overall continuity:

D- + (r2pv ) = 0 (15)
-t 2 Br r

subscriptt r denotes the radial direction)

(2) Energy equation:

9T DT 1 1 3 2 DT 1 N h i
t r Dr pC 2 3r r2 r C i ir ri
p r p i=l

1 C h.w, (16)
pC iJ i i
p i=l

(3) Species continuity:

BY. BY. w.
S+Vr rv 1 2 (r ) i=1,2,..... N (17)
9t r r p p 2 r i ir

(4) Diffusion equation:

ax. N X.X
( V. ) i=1,2,......N (18)
Dr jjr ir
J=l I j ; I

and Y.V. = 0 (19)
j jr
j=1 3 Jr

(18) represents a set of (N-l) linearly independent equations, since

N X.

The unknowns in these equations are: Y., Vir, T, v X. and

P are determined once Y. and T are known respectively, using the

following two equations:

X = 1 1 i=1,2,.......N (20)
I (Y./W )

and p = p-' (the ideal gas equation) (21)

4.1.3 Boundary conditions

The next step is to set up the boundary conditions at the hot

and cold ends for the set of partial differential equations given


a. lne cold boundary difficulty. This has been discussed

considerably in the literature. Basically the difficulty arises from

the incompatibility of the energy conservation equation with the

boundary condition it is subjected to at the cold end. The model for

the flame which assumes both 'steady state' and 'adiabatic' conditions

does not represent truly the physical processes occurring and these

conditions contradict each other for a realistic reaction rate law,

wherein the reaction rate at the cold boundary is small but finite

[6]. To satisfy the conditions of both adiabatic and steady state

requires any one of several assumptions. The assumptions are

mathematically equivalent but have different physical meanings.

Also the validity of these assumptions on physical grounds is open

to question, although these assumptions have been used widely to

solve the flame equations. The assumptions are:

(1) This relaxes the adiabatic condition at the cold boundary

by allowing the existence of a small but finite heat sink at the

boundary in the form of a flame holder. It was introduced originally

by Hirschfelder et al. [3]. This was done because a solution does

not exist when the very small heat loss at the cold boundary is

neglected; the neglect of heat transfer (due to lateral conduction

and radiation) at other points in the flame zone does not, however,

affect the existence of a solution. Mathematically, this is equivalent

to relaxing the condition (T/Dr) = 0 at the cold boundary, and instead

allowing for (3T/r) = E, where E is a small positive quantity. The

value of E determines the 'strength' of the flameholder in terms of

the heat abstracted by it and has to be greater than a certain critical

value to permit the existence of a solution. This assumption is physically

unrealistic for a freely propagating flame which does not have a flame-

holder or heat sink propagating with it.

(2) The second assumption compromises another physical aspect

by allowing for the existence of a temperature (greater than the initial

temperature of the reactants) below which the heat release due to

chemical reactions is identically zero. This temperature has widely

been referred to as the 'ignition temperature' [7]. The choice of

this temperature is again not arbitrary but depends on the particular

flame being studied and it can be determined from the governing

equations. It permits the adiabatic assumption (3T/3r) = 0 at

T = T to be satisfied; however, now this occurs at r = o, as

required for a strict mathematical model [6].

The use of the ignition temperature concept can be explained

by relaxing the adiabatic condition. It is assumed that the heat

loss from the cold boundary and the heat released by chemical

reactions (not identically zero) at the cold boundary reach an

equilibrium at the temperature T Further, assuming the rate of

the chemical reaction to be extremely slow (the initial mixture is

in a retastable state), the change in concentration with time at the

cold boundary can be neglected. Thus it is as though no reaction has

taken place, the heat released at the cold boundary having been

conducted away.

(3) The third assumption involves the use of truncated

kinetics--wherein the reaction rate goes to zero as the temperature

approaches its cold boundary value. Tlis has been shown to be

equivalent to the use of an ignition temperature which is very close

to, but not equal to the cold boundary temperature. Again, this

assumption is not physically possible. However, for most systems

of interest, the reaction rate at the cold boundary is so small as

to be immeasurable, and so it can be considered to be essentially

zero. The reaction rate law is then written as

kf = ATn exp(-E/RTbT) (22)

T- T
where T = u (23)
b u

as compared to the actual law,

k = ATn exp(-E/RT) (24)

The factor T varies from zero at the cold boundary to unity at

the hot boundary. The change in magnitude of the reaction rate due to

this modification becomes smaller as the temperature approaches its

hot boundary value and this is where most of the chemical reaction

and its associated heat release occur.

In the calculations carried out here for the methane-air

flame with an initial temperature of 628 degrees K, it was found

that the chemical reaction rate at this temperature was lower than

the smallest number that can be handled by the IBM 370/165 computer

[i.e., exp(-176)]. In fact, the reaction rates were ignored at

temperatures below 650 degrees K for all the elementary reaction

steps, and above 650 degrees K, the reactions for which E/RTbT was

greater than 120.0 were ignored. The condition E/RTb T 120.0 was

used instead of E/RTbT 2 176.0 to avoid the occurrence of a number

smaller than exp(-176) later during the computations.

Again, when integrating from the hot to the cold boundary,

the integration was stopped arbitrarily for (Fk- Fk-l)/Fk j5 0.005;

(F = dependent variable; temperature and mass fractions) the subscript

'k' refers to the grid point in the r-direction. This had to be done

since, when carrying out numerical calculations, the condition

(8F/ar) = 0 can strictly be achieved only at r = c.

b. Hot boundary. When only the initial (unreacted) and final

(reacted) states of the gas are considered, it is seen that the gas

has passed through the flame zone without any net gain or loss of

heat. All that has occurred is a change in composition due to

chemical reaction and all the heat of the reaction has gone to raise

the temperature of the gas. This adiabatic condition is not fulfilled

for flames which lose significant amounts of heat by other processes,

e.g., radiation, heat losses to flameholders, etc. However, within

the framework of the assumptions used in this analysis, the adiabatic

assumption is satisfied. Again, the adiabatic condition is not

satisfied for any element of gas within the flame zone, since it

loses or gains heat by conduction and diffusion of species. Therefore,

for a constant pressure flame, since the heat transfer equals the

change in enthalpy, the enthalpies (per unit mass of mixture) of only

the initial and the final states are equal.

For the special case of unity Lewis number, the constant

enthalpy condition is satisfied throughout the flame zone. Since

the heat fluxes due to diffusion and thermal conduction are equal

in such a case, any element of gas within the flame zone does not

suffer any net gain or loss of heat and so goes through the flame

zone at constant enthalpy.

The condition of equal enthalpy at the hot and cold ends is

utilized to obtain the boundary condition at the hot end. Since the

gas mixture is in chemical equilibrium at the hot end, the concentra-

tions of the various species are obtained by simply calculating the

equilibrium concentrations corresponding to an isenthalpic combustion

process, starting with a mixture of given concentrations and state

properties at the cold end, i.e.,

h(Tb,p) = h(Tu,p) (25)

which is known from the initial conditions at the cold end. Thus

h(Tu,p) = C (hXiu )TT (26)
i=l u

where the h. are obtained from the JANAF Tables.
Considering a Taylor series expansion about Tbk, the kth

estimate for the value of the burnt gas temperature, and neglecting

higher order terms,

h(Tp) = h(Tbk) + (Tb Tbk) (-h (27)

where the partial derivative is evaluated at (Tbk,p) and

h(Tbk,P) = I (hXi )T=T (28)
i=l b b bk

h(Tb,p) h(Tbk,p)
T = p(Tbk) (29)
P k

Tbk+l = Tbk + ATbk


A first estimation of the temperature at the hot end is

obtained from:

Ah fuel
T = T + (31)
b u C (T) (31)
P u

where Ah = heat of combustion (kcal/gm.mole)

C (T ) = specific heat at constant pressure of the mixture at
temperature T (kcal/gm.mole OK)

= Cpi (Tu)Xi (32)
i=l u

The Ci are again obtained from the JANAF Tables.

The equilibrium concentrations of the species are calculated

at the estimated value of the burnt gas temperature by the method

given in Appendix I. Then equation (30) is used to obtain the next

estimate of this temperature and the process is repeated till

convergence is obtained to within a predetermined accuracy:

A11 T (33)
Tbk -- (33)

= 0.001 was used in the calculations.

This determines the boundary conditions at the hot end:

T = T

X. = X. i=1,2,.......N
i ib

4.2 Method of Solution

An explicit finite difference method was used for numerically

solving the partial differential equations (15-19), and obtain the

profiles of the state properties and the species concentrations

through the flame zone. A scheme involving a symmetrically staggered

mesh in the r-t plane was used. This method has been used by Nash [8]

for the solution of the three-dimensional boundary layer equations.

The staggered mesh scheme uses a set of secondary points in

addition to the primary collocation points, the former lying at the

centers of the rectangles formed by the latter, as shown in Fig. 4.1.

In this method, values of the dependent variables are at first

calculated at the secondary points [at time (t+At/2)] using the values

associated with the primary points at time t Then using values

associated with both the primary points at time t and the secondary

points at time (t+At/2), values of the dependent variables at the

primary points at time (t+At) can be obtained. This is more accurate

than proceeding directly from time level t to level (t+At).

The difference equations for the two steps are readily obtained

by double Taylor series expansions about the appropriate points, and

are [8]:

c 2 (F + F ) + {Ar + T AtM + B + e
c 2 A B [3A B





FAr- Ar/2 t -

Fig. 4.1(a). Staggered mesh scheme, r-t plane


At --


o primary points

o secondary points


1 1
C: (i + k + )
2D: (i 2
i 1
D: (i k + -)
2' 2

Fig. 4.1(b). Enlarged view of staggered mesh


A(i,k) B(i+l,k)

I-----I. I -- I ^ --I--- --C- -I----- -py

~--- --

- .~ -- -- ----- 1--I --------- --- ----

FE = FA +I (1-p)(FD + FC) + X(l-p)Ar L -
SD D l cJ

+ (l+p)At D+ J + e2 (36)

Here F denotes any one of the dependent variables. The quantities

(8F/)t) are obtained from the governing partial differential equations

(16) and (17), the r-derivatives in these equations and in equations (35)

and (36) being expressed in finite difference form. A fourth order dif-

ference formula is used for this purpose:

l-ar1 2
F 1 (F F ) ( F
ri, 12Ar i-2,k i+2,k 3Ar i-l,k i+l,k

1 4 5F
+ I (Ar)4 (37)

In equations (35) and (36), el and e2 denote the discretization

or truncation errors, and contain higher order terms of the double

Taylor series expansions.

e. = ( X )(Ar)2 2 (At)2 t + At(Ar)2 + (38)

1 2 d2F 1 2 D2F
e = (l- X) (Ar) -- (l-0)(At)

1 2 8 F
+ (4X-4Xh + p)At(Ar)2 + ...... (39)
8 tr2

X and p are constants which govern the accuracy and the

stability of the computation. The pair of values X = 1/8 and p = 0

serve to minimize e2 and also assure stability [8]. With these values

the two-step scheme reduces to a one-step scheme:

FC = (FA + FB) + Ar B + At (40)

FE = (FD + FC) + Ar D- C + At D (41)

F(r) = Fu r > ru
4.2.1 Initialization

S-shaped profiles for the dependent variables Y. and T were

assumed to start the integration. Second order curves were chosen and

these were of the form:

F(r) =F r r
u u

(r r)2
F(r) = F + (F Fu) u r r rb +6
b2 (42)
F(r) = F (Fb- Fu) rb+ 6 r > rb

F(r) = F rb r

An arbitrary value of 6 = 2 (r rb) = 0.005 cm. was chosen.

It is immaterial what value of 6 is chosen since the final profiles are

independent of the starting ones.

4.2.2 Details of computation

The initial profiles having been fixed, the integration process

can now be started. Expressing (DX./Or) in finite difference form in


accordance with equation (37), the diffusion equations (18) and (19)

reduce to a set of N simultaneous algebraic equations with the N

diffusion velocities V. (j=l,2,....N) as the unknowns. These

equations were solved using Gauss-Jordan elimination with maximum

pivot strategy. Thus the diffusion velocities were obtained at each

grid point for a given time level t.

Before proceeding further, two other methods for the determina-

tion of diffusion velocities will be mentioned. These methods are

approximate, but involve lesser computation than the method mentioned

above, and may therefore be used where conditions warrant their use.

(1) In flames involving air as oxidant, one component (namely, nitrogen)

is present in great excess so that it is possible to simplify the problem

by regarding each of the other species, one at a time, as forming a

binary mixture with the main component. In general, the effective

binary diffusion coefficient D. for the diffusion of species i in a
mixture m is given by [9]:

S(1/nD)(X.N. X.N.)
1 _= 1 1 .
1 1--- (43)
nD. N
mN. X. ) N.
1 j=l

For the diffusion of trace component i=1,2,3,...N (ifj), in nearly pure

species j,

D. z D.
im li

and so the molar fluxes and hence the diffusion velocities are obtained

by solving the set of algebraic equations:

X. N
N. = nD. X. N (44)
1 im ar 1 (44

Ni = n Vi (45)

(2) This approximation may be used when species i is still a trace,

but there is no single component which may be considered to be in

excess compared to all the others. For such a case the diffusion

velocity of the trace species is obtained from:

1 i
ri dr
Vi = 1dr- (46)
X. Xj/D.
Sj i J j

which shows that the quantity X./Dij is an effective diffusion
coefficient for the species i and the mixture. The above equation is

not valid for the diffusion of a species present in significant amounts.

Such species may be treated in two ways. Often a species which is

present in excess is not significantly affected by the chemical reaction

(e.g., nitrogen) so that its concentration is more or less constant

through the flame zone and hence its diffusion velocity may be neglected.

On the other hand, for species which are chemically reactive, and not

present as a trace, the diffusion velocity may be computed by difference

from equation (19), which is

SYiVi = 0

Getting back to the main problem, it is seen that the mass-

average velocity v is as yet unknown, and a suitable profile for it

was assumed subject to vru = 0, i.e., the unburnt gases are at rest.

This need be done only for integration over the first time step; for

succeeding time steps, the values of vr from the preceding time step

were used for a first approximation. With this assumed profile for

Vr, (3T/at) and (3Yi/3t) (i=1,2,....N) were calculated at each grid

point for a given time level t, using equations (16) and (17). The

derivatives in r occurring in these equations were approximated, again

by the fourth order difference formula (37); values of Vjr (j=1,2,....N)

as calculated above were used, and the JANAF Tables were used to obtain

the values of the state properties h and C The w. (i=1,2,....N) were
p 1
calculated according to equation (8) and modified by truncating the

reaction kinetics as discussed in connection with the cold boundary


Using the ideal gas equation, it may be shown that:

SP, N ti (47)
P W 1 (47)
Si=i i

and thus the overall continuity equation may be rewritten in the form:

Sr 1 Y. 1 2
v -2 lp + T tJ r dr (48)
pr i=l i
(for given time t)

satisfying r = 0.

The values of (DT/at) and (gY./8t) calculated above with an

assumed profile for v were now used in equation (48) to obtain better

values for vr, which were then used to obtain improved values of

(8T/at) and (OY./8t). The iteration can be terminated when the

vr values converge to within a predetermined accuracy.

Actually the values of (DT/8t) between two successive iterations

were compared for checking convergence. This was done since (T/8t) is

used to determine the temperature at the next time level (equations 40,

41) and so convergence of (8T/t) is of greater importance than that

of v which is not used anywhere else in the calculations. The con-

vergence of (8Y./ t) (i=1,2,....N) may also be tested, but these were

not found to vary significantly and the convergence of (aT/9t) generally

assures the convergence of (gY.//t) for all i. A convergence criterion

of the type

At {[-T] (49)
( k _k-1

was used with E = 0.0005. This assures the errors in the calculated

temperatures due to uncertainty in the term (3T/8t) to be less than

E/2 each time either equation (40) or (41) is used. The subscripts

in the above convergence condition refer to the iteration number and

this condition is to be satisfied at each grid point in the flame

zone before proceeding further.

Equation (48) involves quadrature; this was done using the

trapezoidal rule with the same space interval Ar as used for the

other computations. Simpson's rule along with Newton's 3/8 rule may

also be used (IBM subroutine QSF). The difference between the two

was found to be negligible since the space interval Ar is small.

Hence the trapezoidal rule which needs less computer running time

was used.

The term involving (8Y.i/t) in equation (48) has only a small

influence on the final profile of v and may be neglected to a first

approximation. This would amount to making the assumption that the

average molecular weight of the gas mixture is constant through the

flame zone. This is approximately true for methane-air flames in

which a single constituent, nitrogen, accounts for more than 70% of

the concentration at any point in the flame zone. For the sake of

completeness, this term was included in the calculations carried out

here, although its effect is small.

The iterative scheme used here for the solution of equation (48)

is similar to Picard's method of successive approximations for the

solution of a Volterra integral equation [10]. This is seen from the

analysis given below.

Equations (16) and (17) may be written as:

-= -r + G(r) (50)
ot r or

v + H. (r) i=1,2,.....N (51)
at r ar r


( ) 1 1 @ r 2 X BTj
pC 2 ar (r2
G(r) pr
p r


YiV i ar
1 r ar

i- I h.w.
pCp 1 1
p i=1

H.(r) = i (r2 pY.Vi ) i=1,2.....N
1 p p 2 Dr 1 r


At a given instant of time the right hand sides of (52) and (53)

are functions of r only.

Equation (48) now reduces to:

1 f pr2
r vpr

N Y.

i=l 1


N H.(r)
+\-- --
i=l 1

+ (r) dr


Again, for a given instant of time,

N aY
V- N 1 _
_W. Dr
i=l 1

N H.(r)
i=1 1

1 aT
T ar



= A (r)

SA (r)

= A3(r)


where A A2 and A3 are functions of r only. So (54) finally reduces



vr(r) = A3r) A3(r) [vr(r)Al(r) + A2(r)]dr




This is nothing but a Volterra integral equation of the second

kind, and the iterative technique described above is Picard's method of

successive approximations. The first approximation to v (r) in this

method is given by:


rv(r) A3(r) A3(r)A2(r)dr

1 0K { 1 i 1 (r YiVi}
Pr Ii=l pr
r -

N N 3h N
+ 1 1 1 V i 1 2
p r p i=l p i=l


Having finalized the values of the time derivatives (DT/Dt) and

(OYi/t) at a given time instant t, equation (40) was used to proceed

to the secondary points at the time level (t+At/2) and by a similar

process the values of the dependent variables at time level (t+At) were

obtained using equation (41). The mole fractions were then calculated

from the mass fractions using equation (20) and the mean gas density

was obtained from a knowledge of the temperature and the mean molecular

weight using the ideal gas equation. Calculations in the r-direction

at the cold boundary were stopped at any given time level when it was

found that the relative changes in the temperature and the species

mass fractions were less than a predetermined value, i.e.,

F Fk
k F (58)

E = 0.005 was used in the calculations. This process was repeated

until the dependent variables converged to their steady state values.

The transport properties in the conservation equations were determined

(Appendix II) at each grid point every 75 time steps; the error

introduced by not calculating the new values of the transport

properties at every time step would cause small errors in the transient

state solutions, but would not affect the steady state solution. However,

since the estimation of the transport coefficients is quite often approx-

imate, the error introduced by not calculating these coefficients at

every time step, might in fact be lower than the error involved in the

calculations based on the approximate theories.

The last step involves the computation of the laminar flame

velocity based on the unburnt gas density. This was done by considering

a control volume which is large compared to the thickness of the flame

zone, and equating the net rate of production of species i due to

chemical reactions to the rate of outflow of species i. The result is:


f w.dr
u = (59)
uu Y. Y.
lu ib

The values of the flame velocity, S as calculated from the

above equation converged as the state properties and species concentra-

tions approached their steady state values only when species i was

considered to be one of the major species (i.e., CO2, H20, CH4, 02);

however, for each one of the species a different value of the flame

velocity was obtained. The reason for this are the unreliable kinetic

data. Further, this type of behavior has been observed before [11].

The flame velocity calculations based on the other species did not

converge to a steady value, presumably because the term (Yiu Yib)

involved in the denominator of equation (59) is small for these

species and small uncertainties in this term could cause large errors

in the calculation of S These values might, however, converge if

computations are carried on further. In the calculations reported

in Ref. 12, it was found that the flame velocity based on the unstable

species needed more time steps to converge to a steady value.

4.2.3 Step size and stability

The S-shaped profiles used initially were at least qualitatively

correct for the major species (CO2, CH 02, and H20). However, this

was not so for the intermediate species which have a maximum concentra-

tion somewhere in the flame zone and decrease towards both the hot and

the cold boundaries. For this reason initial time steps were restricted

to a very small value until the profiles of the intermediate species

approached their correct form, after which larger step sizes could be

used. A starting step size of At = 1.OE-08 sec. was used and this was

doubled every 50 steps until it reached a maximum value of 32.0E-08 sec.

The step size had to be restricted to this maximum value to satisfy the
stability condition (X/pC )At/(Ar) < 1/2. Since X, p, and C vary
from point to point, the combination of values of and Cp which
from point to point, the combination of values of A, p, and C which

would give the least upper bound for At was used to show that

At < 1.28(Ar)2 which gives At < 32.0E-08 sec. with Ar = 5.0E-04 cm.
max max
Thus this explicit scheme is conditionally stable.

The use of these values of At, however, caused instability in

the integration of the term involving the chemical rate of production

of species. So this term was integrated separately with an appropriate

value of step size. The integration can be done by Runge-Kutta's fourth

order method for a set of first order differential equations. The


i. dn.
dt -fi(nn2'.....N) i=1,2,.....N (60)

form a set of what is termed as 'stiff'* equations, and such equations

need very small step sizes and hence prohibitively long computer running

times with Runge-Kutta's method.

Gear [13] has developed an algorithm for dealing with stiff equa-

tions, and this was used in the integration of the reaction rate equations

(details of the method in Appendix III). Instead of starting the integra-

tion by a first order method, as proposed by Gear, a higher method was

used. This necessitated the calculations of up to six derivatives of ni,

since the maximum order for which the algorithm is devised is six.

Generally, it was found in this problem that a starting fourth order

method resulted in the maximum value of the step size. The order for

subsequent time steps was decided by the program itself. The following

equations were used to obtain the derivatives.

So called because this behavior was first noticed in connection
with servomechanisms having a tight coupling between the driving and the
driven member.

dn. M
di (i i)[Pj P2j]
dt 1 i 1i 1i 2j

d n. M
2 (v.. -v R)[P P -P .]
dt2 j=l 1j Ij3j P2jP4j
dt j=1

d ( P R +2 p 2 P
t j (1 i ii 3j lj 5j P2j4j 2j6j
dt j=1

d n. M R 3
SVP R p3
-7-= (.. .)[P .P + 3P .P .P + P .P
4 13 13 ij )[Plj 3j lj3j5j + 7
dt j-1

- P p3
2j 4j

- 3P 2P .P
2j 43 6j





- P2jPj ]

d n. 4 2 2
__ ( P )[p P + 6P P P + 3P P + 4P 2P
dt5 j=1 i j i lj3j l+ 3j 5j + 53I j + lj 3jP7j
dt j=1

4 2 2
+ .P -p P.4 6P.P -3 2 P 4P2.P.P.
-lj j P2j 4j 6P2j 46j 3P2j 6j 2jP4jP8j

- P2jP j I
2j 10j


d n. M
i P 5 3 2
1t 16 j 1 3j 1j 3j 5 j j3j 5j
dt j=l

+ o10P .P2 + 10P .P .P. + 5P- .P. + P .P-
1 P 3j 7j + 1P 5j7j jj + PjP 1j

5 3 2
P.P 10P.P .P 15P .P .P
-2j 4j 2j 4j 6j 52jP4j 6j

O1P p2 P -5P P(60f)
2j 4j 8j 2jP6j 8j 5P2j lOj -2j 12j

Pj = kfj i=

P2 = kb 1
2j i=l

3j r=1

P =
P5j E =l


7j =l

P9j = zj Q4Z

Plj= -1

kj li

Vj Q2Z

j Q3Z

(60i) P4j = 1 Q

(60k) P6j = Z j Q2

(60m) P8 =
OJ R=1

V j Q3

(60o) P10j= k 1 j Q4

k JQ5e

1 dn
QI n dt

(60q) Pl2j = I
3 =i

(60s) Q2Z


vj Q5Z

1 d2n
n dt2

dn 3 3 dn
dt J 2dt

dn] 2
d- j

Idn 4
dt .

S1 d 4
"n dt4


[n ij

V. .
[ni] j




(60( )







Q3~, 3

- -

1 d3n

Sn dt3


d t2

2 2

dt 2


d 3t

4 dn
2 dt


dn} 2
dt )

n k

24 dn 5 60 dnl 3 d2n 30 dn d n 2 20 d 2 d3n
S0 + +20
5 5 dt 4 dt J 2 3 d dt 2 3 dt 3
d 2 n d nn d dt n dt n dt

10 dt d3n] 5 dn d4n 1 dn
2 1d J [dt P 2 dt + (60w)
n2 3 2 dt 4 n dt5
n dt dt n dt dt

The calculation of these derivatives requires a knowledge of

the temperature to determine the values of the rate constants involved

in these equations. The temperature was assigned its value at the

beginning of a time step of the main integration and regarded as

constant over this main time step. The computations might be repeated,

however, using a suitable average value once an estimate has been

obtained for the temperature at the end of the time step. Instead of

using a simple arithmetic average T = 1 (T + T2), an average of the
nature exp(-E/TmR) = [exp(-E/RT1) + exp(-E/RT2)], which in fact

averages the rate constant as opposed to the temperature might be used.

(E = activation energy of the rate controlling step.) This iteration

across a time step would not significantly alter the results, and so

this was not done in this study.

The term h.w. in the energy equation was treated similarly.

Since the temperature was assigned its value at the beginning of a

main time step and regarded as constant over the main time step, the

h. (i=1,2,......N) are also constant over this main step and so,

h.w. dt = h w. dk (61)
i=l x i=l

So the result of the previous integration was used here instead of

resorting to a separate integration of the term X h.w.. Again, a

mean value for the h. can be used for improved accuracy.

4.3 Chemical Kinetics

The conservation equations contain terms involving the rate of

production of species due to chemical reactions. These terms must be

evaluated before the equations can be solved. The evaluation of these

terms requires a knowledge of the reaction mechanism and the rate

constant parameters for the different elementary reaction steps.

Rate constants for the elementary reactions occurring in

the C-O-H-N system have been tabulated by several authors. However,

there is considerable disagreement over the values of the rate

constants from author to author. Again, there is disagreement over

the reaction mechanism proposed for the oxidation of methane. In

this analysis a total of 16 species were considered. It was found

that it is possible to write a total of 84 (42 forward and 42 reverse)

bimolecular elementary reaction steps for these 16 species. In addition

to these, several third order reactions involving a third body are also

possible. The kinetic data for all these reactions are not available;

moreover it is not necessary to consider all the reactions for which

data are available. A detailed analysis of all the bimolecular

reactions for which data are available was made to determine the

contribution of each reaction towards the total rate of production

of each of the species. This was done over the temperature range

600-24000K. On the basis of this, reactions which contributed less

CO2, CO, N2, H20, H, OH, O, NO, NO2, N20, N, 02, H2, CH3,

than 0.1% to the total rate of production of any of the species over

the entire temperature range were eliminated from further consideration.

The contribution of each reaction from an energy standpoint was not

tested; it being assumed that the reactions which contribute less

than 0.1% to the total rate of production of any of the species

would not cause any significant contribution to the energy equation

also, unless their heats of reaction are at least an order of

magnitude more than those of the reactions which were not eliminated

from consideration. A similar analysis was also done for the third

order reactions. Finally, a mechanism consisting of 25 elementary

reaction steps (Table 4.1) was used for the analysis. The first

14 reactions are the ones which Seery and Bowman [14] have suggested

in their proposed scheme for the oxidation of methane. The remaining

reactions are concerned with the N-0 system.

The kinetic data used in this analysis have by no means been

firmly established, and suggested error limits for some of the rate

constants are as high as 25%. The error is greater the higher the

temperature, since most of the experimental data are gathered in

the temperature range 300-1000 degrees K and involve extrapolation

to higher temperatures. Also the oxidation mechanism of methane

has not been firmly established and much work remains to be done on

the effects of pressure, temperature, mixture composition, and the

presence of inerts on the oxidation mechanism. In the absence of

better data and a more definite reaction mechanism, the data referred

to above were used. The method developed however is general, and as


Table 4.1

Reactions considered in this analysis

k = ATn exp(-E/RT)

Forward rate constant



CH4 + M = CH3 + H + M
4 3
OH + CHO = CO + H20
0 + CH = HO + CHO
CHO + M = CO + H + M
0 + CH4 = OH + CH3
4 3
CH4 + H = H2 + CH3
OH + CH4 = H20 + CH3
4 2 3
OH + CO = CO + H
H + 0 = OH + 0
0 + H2 = H + OH
0 + H20 = OH + OH
H + H20 = H2 + OH
M + H20 = H + OH + M
0 + CH = HCO + H

N + 0 = NO + N
NO + 0 = 0 + N
N + 0 = NO + NO
NO + 0 = 0 + NO
M + NO = 0 + NO + M
0 + N20 = NO + NO
NO + NO = NO + 02
NO + N2 + 02 = N20 + NO2
OH + M = 0 + H + M
H + NO = OH + N
N2 + OH = H + N20




Units: E/R in degrees K N-
rate constants in (cc/g.mole) /sec.
where N = order of the reaction







Table 4.1 (Continued)

Reactions which were not considered because of negligible contribution

Forward rate constant

Reaction A n E/R Ref.

N + NO = NO + N 0

N2 + 02 = 0 + N20
NO + H = H20 + NO
OH + NO = H + NO2
0 + H = OH + OH
0 +N20= NO2 + N
N + M = N + N +M
NO2 + NO = NO + NO + 02
NO + N + 0 = N20 + NO2
0 + 0 + M = 0 + M
CO + NO = NO + CO2
NO + NO = NO + N
NO + M = 0 + N + M
NO + M = N +N 02 + M

N2 + 0 = NO2 + N
2 2 2

1. 0Oxl14







Reactions for which kinetic data were not available

CO2 + N2 = CO + N20
CO2 + OH = 02 + CHO
CO2 + N CO + NO
CO + OH = 0 + CHO
CO + CHi = CH3 + CHO
N2 + NO = N20 + N
H20 + N = NO + H2

CO2 + H = O + CHO
CO2 + 0 = CO + 02
CO2 + 2 = CO + H 0
CO + H2 = H + CHO
N2 + H20 = N20 + H2

H20 + 0 = 02 + H2

Most of these reactions are not of much significance from a

practical standpoint.

Table 4.1


more is known about the oxidation mechanism and the rate constants,

the new information can be readily incorporated to obtain more

realistic solutions. In the computer program developed any number

of species and elementary reaction steps may be considered, the only

limitation being the computer time and storage available.

4.4 Uncertainty in the Definition of Temperature

It is now generally realized that there does exist a distinct

possibility that the various chemical species in a flame zone may not

be in equilibrium with respect to their various degrees of freedom.

Thus the concept of a 'temperature' to define the thermodynamic

state of the gases in the flame zone is rendered difficult. This

problem is briefly discussed here; details may be found in the

appropriate literature.

The statistical definition of temperature is based on the

Maxwell-Boltzmann equilibrium energy distribution function. It is

defined as a parameter characteristic of a system at equilibrium

and it specifies the distribution of energy between the particles

comprising the system. Normally the definition is based on the

distribution of the translational kinetic energy. Since polyatomic

molecules have in general other modes of energy storage, it is

possible to analogously define various temperatures, e.g., Tib'

Tro T elec. For a system at complete equilibrium all these
rot' elec
definitions give a unique value of temperature, the various degrees

of freedom being at equilibrium with each other.

The kinetic theory definition of temperature for a system

composed of different chemical species is expressed in terms of the

mean particle energy, averaged over all the species [3, p. 455]:

3 kT= 1 1 n. ( mV2
2 jl 2 i i

Temperature for the other degrees of freedom may be defined

analogously. For a multi-component system not at equilibrium, it is

possible that not only are the various degrees of freedom of a

particular species not in equilibrium with each other, but that there

is non-equilibrium between the different chemical species, i.e., it

may not be possible to define a unique translational (or rotational

or vibrational) temperature for the entire system.

In a laminar premixed flame, the temperature of the gas is

raised extremely rapidly (more so in detonations and shock waves).

When this happens the different degrees of freedom would in general

come to equilibrium at different rates. Translational equilibrium

is attained the fastest--in about five mean free paths, as demonstrated

by shock-wave experiments in gases with no internal degrees of freedom

(Ne, Ar, Kr). The temperature of primary concern here is the trans-

lational one, since this is used in the equation of state for

transport property calculations and in the evaluation of reaction

rate constants. Since the assumption of a continuum has already

been made, i.e., changes in macroscopic properties occur over

distances, large compared to a mean free path, the concept of a

local translational temperature should be valid.

The rotational degrees of freedom in diatomic and polyatomic

molecules also come to equilibrium at a rapid rate though not as fast

as the translational mode. It is the vibrational degrees) of freedom

which poses problems. Equilibration or relaxation times for the

vibrational degrees of freedom are in general much larger than for

the translational or rotational degrees. The reason for this is that

energy differences between quantized vibrational energy levels are

much greater than for the translational or rotational levels. As

a result, when the temperature is being raised rapidly, the number

of collisions needed to excite the next higher vibrational energy

level may be significant, of the order of 5000, compared to 50 for

translational and rotational equilibrium. For triatomic molecules

(N20, CO2) it may be more, about 2500-50,000. Thus the translational

and rotational temperatures of a gas being heated rapidly follow

their equilibrium values closely while the vibrational temperature

may still be its cold-gas value. In such cases a strict unique

definition of temperature is rendered impossible. It may be neces-

sary to define two different temperatures--one translational and the

other vibrational.

Thus the definition of temperature is linked to the requirement

of equilibrium. If there is non-equilibrium between the different

chemical species or the different degrees of freedom of a particular

species, the Maxwell-Boltzmann equilibrium distribution function

ceases to hold, and so, strictly speaking, a temperature cannot be

defined for the system. To use the concept of temperature for a

non-equilibrium system, the assumption generally made is that even

though the system as a whole is not at equilibrium, local equilibrium

exists at any particular point in the system. This boils down to

assuming that in any elementary volume (whose dimensions are 'small'

compared to those of the system) the particles of the different species

are in equilibrium with respect to each other and their individual

degrees of freedom, so that the Maxwell-Boltzmann energy distribution

may be used to define a local temperature. This is how the use of

temperature in flame analysis is justified.

In flames where the transfer of energy between the vibrational

and translational energy modes is extremely slow, the need for separate

definitions of translational and vibrational temperatures may destroy

the concept of a local temperature also. Such cases have been shown

to exist by spectroscopic means [20]. The chemical reaction rates

and the transport properties are evaluated at values of temperature

based on a static system where local equilibrium has been achieved.

For the case just mentioned these evaluations may be somewhat different

due to non-equilibrium effects; it is not known whether the difference

is significant in all cases. Some work on the high temperature

oxidation of hydrogen and the dissociation of oxygen has been

reported. The results show that the slowness of vibrational relaxa-

tion does change the processes to a certain extent; the chemical

reactions occur at 'temperatures' different from their local

equilibrium values.

For chemical reactions which are extremely rapid (half-lives

of 10 microseconds or less) the classical assumption of equipartition


of energy is invalid. Observations for some rapid reactions have

revealed that some of the reaction products are formed almost

entirely in higher vibrational levels (e.g., 14th or 15th) and the

transition to the equilibrium ground state occurs rather slowly [21].


1. Spalding, D. B., "The theory of flame phenomena with a chain
reaction," Phil. Trans. Roy. Soc. London, A249, 1 (1956).

2. Williams, F. A., "Combustion Theory," Addison-Wesley (1965).

3. Hirschfelder, J. 0., Curtiss, C. F., and Bird, R. B., "Molecular
theory of gases and liquids," John Wiley and Sons, New York

4. Pritchard, H. 0., "The kinetics of dissociation of a diatomic
gas," J. Phys. Chem., 65, 504 (1961).

5. Rice, 0. K., "On the relation between an equilibrium constant
and the non-equilibrium rate constants of direct and reverse
reactions," J. Phys. Chem., 65, 1972 (1961).

6. Berlad, A. L. and Yang, C. H., "On the existence of steady state
flames," Combustion and Flame, 3, 447 (1959).

7. Von Karman, T. and Penner, S. S., "Selected combustion problems,"
AGARD, Vol. I, Butterworths, London (1954).

8. Nash, J. F., "An explicit scheme for the calculation of three-
dimensional turbulent boundary layers," Trans. of ASIE, Journal
of Basic Eng., 94, 131 (1972).

9. Bird, R. B., Stewart, W. E., and Lightfoot, E. N., "Transport
Phenomena," John Wiley and Sons, New York (1960).

10. Tricomi, F. G., "Integral Equations," Interscience Publishers,
New York (1970).

11. Bledjian, L., "Conputation of time-dependent laminar flame
structure," Combustion and Flame, 20, 5 (1973).

12. Spalding, D. B., Stephenson, P. L., and Taylor, R. G., "A
calculation procedure for the prediction of laminar flame
speeds," Combustion and Flame, 17, 55 (1971).

13. Gear, C. W., "The automatic integration of ordinary differential
equations," Comm. of ACH, 14, 176 (1971).

14. Seery, D. J. and Bowman, C. T., "An experimental and analytical
study of methane oxidation behind shock waves," Combustion and
Flame, 14, 37 (1970).

15. Baulch, D. L., Drysdale, D. D., and Lloyd, A. C., "Critical
evaluation of rate data for homogeneous gas phase reactions of
interest in high temperature systems," Dept. of Physical
Chemistry, The University, Leeds, England 4lay 1968-July 1970).

16. Bahn, G. S., "Reaction rate compilation for the H-O-N system,"
Gordon and Breach Science Publishers, New York (1968).

17. Bowman, C. T., "An experimental and analytical investigation of
the high temperature oxidation mechanisms of hydrocarbon fuels,"
Combustion Sci. Tech., 2, 161 (1970).

18. Newhall, H. K., "Kinetics of engine generated nitrogen oxide and
carbon monoxide," Twelfth Symp. (Int'l) on Combustion, The
Combustion Institute, Pittsburgh, Pa. (1969).

19. Ellis, G. and Bahn, G., "Engineering selection of reaction rate
constants for gaseous chemical species at high temperatures,"
The Combustion Institute, Western States Section, Paper 62-27

20. Broida, H. P., "Temperature: Its measurement and control in
science and industry," Vol. II, Reinhold Publishing Corp.,
New York (1955).

21. Strehlow, R. A., "Fundamentals of Combustion," International
Textbook Co., Scranton, Pa. (1968).


5.1 Introduction

The Otto cycle may be envisioned as a sequence of thermodynamic

processes, during which the thermodynamic properties and the chemical

composition of the working fluid undergo changes. The cycle starts

with the intake of a fuel-air mixture at state 1 (Fig. 5.1). The

mixture is compressed until state 2, where it is ignited with an

electric spark. Process 2-3 represents the combustion and 3-4 is the

expansion process at the end of which the gases are expelled from the

engine cylinder and a fresh charge drawn in.

The early theoretical analyses made several assumptions to

simplify the problem and make it solvable without a great deal of

computational work. For example, it is common to assume that

processes 1-2 and 3-4 are isentropic and that the combustion process

2-3 and exhaust process 4-1 occur at constant volume. Further

simplifications are obtained through the assumption that air is the

working fluid throughout the cycle and its properties (specific heats)

are those of atmospheric air (known as the air standard cycle). Recent

analyses have removed several of these simplifying assumptions and

modeled the processes occurring in the Otto cycle engine on a more

realistic basis. This has largely been possible due to the availability

of high-speed electronic computing machines. A review of the early

1 1



Fig. 5.1. p-V diagram for the Otto cycle

- -r~--------------- I'--- ---f


~-1Ur*LII ---^ II_--_ 1__1_

work is given in Ref. 1. Among some of the recent work considering

the complex sequence of processes to some detail are those of Edson

[2], Spadaccini [31, and others.

Edson's analysis was made on the assumption of chemical

equilibrium at every step in the cycle; the main aim, however, was

to study the influence of compression ratio and dissociation on the

thermal efficiency. This assumption of chemical equilibrium has

been made in several analyses due to the lack of chemical kinetics


Spadaccini's analysis [3] has been one of the most detailed

analyses made so far. In this analysis the compression process is

assumed to be isentropic and the specific heat is assumed to be a

linear function of temperature for this process. Combustion is

assumed to occur at constant volume and the composition and

thermodynamic properties of the mixture after combustion are found

on the basis of a chemical equilibrium analysis. For the expansion

process, a non-equilibrium finite rate kinetics analysis with

allowances for heat transfer is made considering a total of 13

chemical species to be present. It will be seen later that an

equilibrium analysis of the expansion process does not convey the

true state of affairs with respect to the concentrations of the

species at exhaust. Since the present interest is to study the

exhaust emissions, a non-equilibrium analysis is a must to obtain a

realistic estimate of the emission levels. Other work in this area

has been referred to at the appropriate places in this chapter.

5.2 Present Work

Here again the compression process has been assumed to be

isentropic. The combustion process is initiated sometime before top

dead center is reached and completed shortly after top dead center.

Instead of assuming the combustion process to be instantaneous as

done in many of the previous analyses, here the process has been

assumed to proceed at a finite rate. Using an assumed burning rate

law which gives the mass fraction of the total fuel-air mixture burnt

as a function of the crank angle, a step-by-step analysis of the

combustion process has been made to obtain the state properties and

the chemical composition of the working fluid at the point where

combustion is complete. After this, for the expansion process analysis,

the same reaction mechanism as used by Spadaccini has been used here;

however, different mathematical techniques have been used to solve the

governing equations.

The concentrations of the various species at exhaust are thus

obtained and so this model can be used to study the influence of several

parameters such as air-fuel ratio, exhaust gas recirculation, compression

ratio, speed, ignition timing, etc., on the exhaust emissions. Due to

the approximations involved in the theoretical analysis, and the fact

that the processes occurring in the exhaust pipe have not been considered

in this analysis, agreement with the experimentally measured values may

not be exact. However, the trend in the results will definitely show

up in this theoretical analysis and so experimental runs need be made

to obtain results only for a few combinations of the different

parameters which seem promising from the emissions viewpoint, after

a preliminary detailed theoretical analysis.

The only empirical input here is the burning rate law. The

validity of the burning rate law can be checked by comparing the

calculated thermodynamic properties (i.e., the pressure and the

temperature) with the experimentally measured values. The law may

then be modified, if necessary, to obtain better agreement with

experimental results. The assumptions made here agree qualitatively

with the experimental results of Refs. 4 and 5 and other published

literature; quantitative agreement of the mass fraction burnt as a

function of the crank angle could not be checked due to differences

in the operating conditions.

A rigorous theoretical analysis involving the solution of the

general three-dimensional time-dependent conservation equations for

a flame propagating into a premixed fuel-air mixture would yield the

mass fraction burnt as a function of the crank angle. However, with

present day computational techniques, this type of an analysis of

the problem is almost impossible. Moreover, the analysis is beset

by difficulties due to the limited knowledge of the complex nature

of the processes occurring in the engine: the reaction mechanism

and the rate constant values, degree of turbulence in the combustion

chamber, problems arising from the fuel-air mixture not being homo-

geneous and the fuel not being completely vaporized, heat transfer

at the combustion chamber walls and the associated quenching effects

on the flame, together with a host of other difficulties.

The effect of the injection of alcohol-water mixtures on the

emissions has been studied using the model. In general, many of the

control techniques which have been suggested involve some treatment

of the exhaust gases. The injection of water-alcohol mixtures strives

to reduce the formation of the pollutants in the engine itself rather

than offer some treatment of the exhaust gases.

The hydrocarbon emissions result mainly from the heterogeneous

processes occurring in the quench zone and are influenced by several

other factors. Due to the lack of knowledge of the complex processes

occurring, it was not possible to predict the hydrocarbon emissions

using this model. However, it should be remembered that control

techniques which result in the reduction of carbon monoxide also

generally result in the reduction of hydrocarbons. So the latter need

not be studied separately, but can be determined by subsequent experiments

after the preliminary theoretical analysis.

5.3 Major Assumptions for Analysis

Before describing the mathematical model of the Otto cycle,

a discussion of the assumptions made in the analysis is in order.

In this section the assumptions relating to the entire cycle are

mentioned. Assumptions for the specific processes are mentioned in

the appropriate places in this presentation.

(a) The fuel is assumed to be completely vaporized and mixed homo-

geneously with air. If the fuel-air mixture is not homogeneous, it

will be locally fuel-rich at some points and locally fuel-lean at

others, though it may be overall stoichiometric. Compared to

stoichiometric combustion, the fuel-rich parts would yield larger

amounts of unburnt hydrocarbons and carbon monoxide while the

fuel-lean parts would yield a larger amount of the oxides of

nitrogen. Also the energy release would be less since a lesser

amount of the fuel is actually burnt. The result would be (as

compared to the homogeneously mixed case):

(i) relatively low mean exhaust temperature, and

(ii) higher concentrations of carbon monoxide, unburnt

hydrocarbons, and oxides of nitrogen.

Similar conclusions would hold when the mixture was overall

fuel-lean or fuel-rich, and not homogeneous. Thus the theoretical

calculations assuming perfect mixing represent the upper limit of

performance and the lower limit of the emissions.

(b) The gas mixture is assumed to be thermally perfect. Nitrogen

is a major constituent (> 70%) of the mixture and it closely

approximates thermally perfect behavior for the conditions

encountered in a typical Otto cycle engine. Other constituents are

present in relatively smaller concentrations, so as not to affect the

overall thermally perfect behavior. The only constituent likely to

be considerably thermally imperfect would be a fuel like iso-octane.

However, since its concentration in the mixture is small (mole fraction

< 0.02), the results would be affected only slightly.

(c) Reactions in the unburnt gas region are neglected, i.e., there

are no precombustion reactions. The maximum temperature in the

unburnt gas region is generally less than 1000 degrees K for normal

combustion, so this is a fairly good assumption. Some reactions

(i.e., low temperature pyrolysis, slow oxidation, cool flames) do

occur in this region in practice, but the rates of these reactions

are significant only at temperatures above about 1400 degrees K.

(d) The quench layer is assumed to occupy a negligible fraction of

the total volume. The quench layer is a major source of the hydro-

carbons and products of incomplete combustion present in the engine

exhaust. This assumption implies that only the homogeneous chemical

and physical processes occurring in the bulk of the gases are

considered and the heterogeneous processes occurring near the combustion

chamber walls are neglected.

(e) There is no heat flow or mixing between adjacent elements of gas

during combustion, i.e., transport phenomena of viscosity, thermal

conductivity (also radiation and convection) are neglected. Transport

phenomena are of importance only in the flame element where appreciable

gradients of temperature and species concentrations exist, and these

have been considered in the analysis of the flame zone presented

earlier in this study. In the unburnt gas region, gradients in these

properties are generally absent, since this region is fairly homogeneous.

In the burnt gas region, gradients in these properties are present, but

are not appreciable enough to cause significant fluxes of heat and mass.

Besides, in this study a bulk analysis has been made considering the

entire burnt gas region to be at a uniform temperature and concentration.

(f) No pressure gradients exist in the combustion chamber. Although

the pressure changes with time, it is spatially uniform within the

combustion chamber; any pressure gradient introduced would be

Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd