NBSIR 76-994

A Numerical Solution of the

Time Dependent Partial

Differential Equations which

describe a one-dimensional.

Laminar, Premixed Flame

R. L. Brown

Institute for Materials Research

National Bureau of Standards

Washington, D. C. 20234

January 1976

Final

U. S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS

NBSIR 76-994

A NUMERICAL SOLUTION OF THE

TIME DEPENDENT PARTIAL

DIFFERENTIAL EQUATIONS WHICH

DESCRIBE A ONE-DIMENSIONAL

LAMINAR, PREMIXED FLAME

R. L. Brown

Institute for Materials Research

National Bureau of Standards

Washington, D. C. 20234

January 1976

Final

U.S. DEPARTMENT OF COMMERCE, Elliot L. Richardson, Secretary

Edward O. Vetter, Under Secretary

Dr. Betsy Ancker-Johnson, Assistartt Secretary for Science and Technology

NATIONAL BUREAU OF STANDARDS. Ernest Ambler, Acting Director

A NO^imCAL SOLUriON OF THE TIME • DEPENDENT PARTIAL DIFFERENTIAL

EQUATIONS WHICH DESCRIBE A ONE-DIMENSIONAL, LAMINM, PREMIXED FLAME

by

R. L. Brown

Institute for Materials Research

National Bureau of Standards

Washington, D. C, 2023-4

ABSTRACT

The set of time dependent, parabolic differential equations, which

describe the physical and chemical processes in a one-dimensional, laminar, premixed flame is solved by adapting a solution procedure originally developed to solve the two-dimensional steady state boundary layer

equations. The flame equations are integrated by an implicit method

Until the steady state is reached. This corresponds to a flame propagating steadily through a mixture of combustible gases. By a suitable

choice of boundary conditions, it is also possible to model a flame

which is stabilized on a burner. Solution of the flame equations yields

the concentration profiles of the different chemical species as well as

the temperature profile. From these one can also calculate the production rates of each species, the rate of each chemical reaction, and

the heat release rate at each point in the flame. The velocity of the

freely propagatinK flame can be calculated from the integrals over the

whole flame zone of uny of bho .''pecies prodU'itlon rates. The model

incorporates roylistic t?iorinodynamic data and transport property data

that are functions of both temperature and concentration. A complete

documentation of the computer program which accomplishes the integration

is presented.

TABLE OF CONTENTS

I. INTRODUCTION 1

IE. METHOD USED TO SOLVE THE FLAME EQUATIONS 3

III. DETAILED DESCRIPTION OF THE PROGRAM 9

A. Description of the Main Program SPALD/3 9

B. Description of Subroutine SWITCH( KRATE) 29

C. Description of Subroutine DIFUSE(PRE3S) 31

D. Description of Subroutine OlMEGA 33

E. Description of Subroutine INITL/l

F. Description of Subroutine TRANS( ICALC)/6 38

G. Description of Subroutine RATCN/3 LH

H. Description of Subroutine ENTRN( J, IFOW) 51

I. Description of Subroutine SCHMl/l 52

J. Description of Subroutine SCHM2(J)/7 53

K. Description of Subroutine CALC(J) 58

L. Description of Subroutine SENTP/l 60

M. Description of Subroutine ZCALC/l 62

N. Description of Subroutine OUTPUr(L,l[EADNG,RUliID,X, INPUT,

LMAX) 63

IV. TEST CASE 7/^

LISTING OF THE PROGRAM 84.

R.EFEREi^JCES 131

APPENDIX A. DERIVATION OF THE ONE-DIMENSIONAL LAMINAR FLAME

EQUATIONS A1

a) Conservation of Total Mass A1

b) Spocior. Conservation Equations A2

c) Energy t'onservation Equ-ition A3
APPEimX B. DETERMINATION OF THE VELOCITY OF A STEADILI

PROPAGATING FLAME kU

APPENDIX 0. CALCULATION OF DIFFUSION VELOCITIES A6

APPENDIX D. CALCULATION OF THE THERMAL CONDUCTIVITY OF THE FLAME. A8

APPENDIX E. VON MISES TRANSFORMATION OF THE FLAME EQUATIONS A9

APPENDIX E. SPALDING'S TRANSFORMATION OF THE FLAME EQUATIONS All

APPENDIX G. U.\[ITY LEWIS NUl^ER APPROXIMATION A1 3

APPENDIX H. FINITE-DIFFERENCE FORMS OF THE FLAME EQU7ITIONS A1

5

APPENDIX I. CALGUIJITION OF THE SOUltCE TERMS A20

a) Calculation of the chemical source term A21

b) Calculation of the energy source tertri A24.

c) Calculation of source terms arising from diffusion of

molecular weight k2U

APPENDIX J. CALCULATION OF THE ENTRAINMENT RATES A26

APPENDIX K. BOU^IDARY CONDITIONS A28

a) Boundary conditions for a freely propagating flame A28

b) Boundary conditions for a flame stabilized on a burner. . . A28

APPENDIX L. TRANSFORMATION FROM THE au,t COORDINATE SYSTEM BACK

TO THE y,t SYSTEM A33

APPENDIX M. CALCUUTION OF THERMAL AND TRANSPORT PROPERTIES FOR

PURE SPECIES A3

5

a) Binary diffusion coefficients A35

b) Thermal conductivities A36

c) Heat Capacities A37

d) Enthalpies A37

APPENDIX N, INITIAL PROFILES A38
1

I. INTRODUCTION

* The factors which determine how fast a flame will move through a

mixtiJre of combustible gases have been known for many years. Flames of

practical interest are generated by chemical reactions involving free

radicals. The radical reactions usually have high activation energies

and thus are rapid only at elevated temperatures. A flame propagates

because the cold gases in its path are raised in temperature by thermal

conduction and infused with radicals by mass diffusion from the flame

front. It is a simple matter to set up eqiiations for the conservation

of energy and mass which give the time dependence of the temperature

and concentration of chemical species at different positions in the

flame. These equations contain transport parameters and chemical rate

constants. The transport parameters consist of mass diffusion coefficients for each chemical species in the mixture and the thermal conductivity of the mixture. The kinetic theory of gases provides the

functional dependence of these transport properties on the composition

and temperature of the mixture. In addition, it also requires a knowledge of binary diffusion coefficients between pairs of species and the

thermal conductivities of the pure species. Although the mathematical

expressions for the transport properties are very complex, this constitutes no impediment to the solution of the flame equations; their changes

with time are slow and do not require calculation at every integration

step. Furthermore, the kinetic theory is sufficiently accurate, and our

experimental and theoretical knowledge of diffusion coefficients and

thermal conductivities sufficiently good that the transport phenomena

impose no major barrier to a flame modeling calculation.

The significant barrier arises from a lack of accurate rate constants

and until the last few years in efficient numerical techniques for solving the flame equations. While the rate constant problem remains, it is

now possible to solve the flame equations on a high speed digital computer in times ranging from a few seconds to several minutes depending

on the complexity of the flame. The method used in the present approach

1 was originally developed by Patankar and Spalding to solve the twodimensionnl steady -state boundary layer e(|Uatlons. One of the two
2

spatial variables in that problem becomes time in the flame equations,

while the other represents the distance through the flame fromt. The

equations for the two systems are otherwise the same.

In reading this report, it is desirable to begin with Section II.

This gives a brief discussion of the method used to solve the flame

equations. Appendices A and B should then be read. These give a simple

derivation of the one-dimensional flame equations and show how the

steady-state flame velocity is calculated from the solutions of the

equations. These equations contain the diffusion velocities of the

various species. The expressions used for these are derived in Appendix

C. We have considered here only the diffusion arising from a concentra*-

tion gradient. That caused by a thermal gradient has been neglected.

Appendix D, which shows how the thermal conductivity of the flame was

calciilated can be skipped on a first reading. Appendices E and F show

how the flame equations are transformed to a new spatial variable.

This is done in two steps. The first is the Von Mises transformation

which eliminates the equation for overall mass conservation. This removes the explicit dependence on the flame velocity in the remaining

equations. Next, a transformation devised by Patankar and Spalding is

Used to convert the equations into the form of the boundary layer equations. Appendix G discusses the assumptions underlying the unity Lewis

n'Jmber approximation and can be skipped. Appendices H through K show

how the equations are written in finite-difference form, how the various

source terms are calculated, the way the grid size is controlled during

the integration of the equations, and how the boundary conditions are

incorporated into the calculation. The last three Appendices can be

skipped.

After this, one can begin with Section III, which gives a detailed

discussion of how the program works. Throughout this section, references

are made to relevent parts of the Appendices for detailed discussions

of the various calculations.

A test case is presented in Section IV. For this a H^-Br^ flame

was used. This shows the form of the input data and the type of output

thrit the program yields.
3

II. METHOD USED TO SOI.TO THE FLAME EQUATIONS

In this section we present an outline of the method used to solve

the flame equations. The details are presented in the Appendices.

The equations which describe a one -dimensional pre-mixed laminar

flame can be vrritten in the form of the classical diffusion equation.

at

1, (1)

X

C. is the concentration of the j'th species, t . is a diffusion

J J parameter and S. is the chemical source term for 1. Both t. and S. .J J J are functions (presumed to be known) of the concentrations of the various

species. The independent variables are the time t and a distance x.

X is actually a transformed spatial variable which is roughly proportional to the distance y which would be measured in the laboratory,

(x = \p , see Appendix E. ) There is one of these equations for each

species in the flame, and also a similar equation for the temperature

or the enthalpy. The derivation of these equations is given in Appendices

A and E.

To solve these equations we set up a two-dimensional grid for the

independent variables as shown in Fig. 1. The subscript on C in the

figure now denotes a particular

grid point on the x axis and is

not a species index. The primes

denote the values of the C at

m

the time t' . The differential

equation (l) gives the value of

ac/at in terms of the derivatives of the spatial variable.

Thus, at a particular grid point

m we have,

X

ra+1

m

m-1

C

m

m-1

T

6t.

Figure 1 . Grid Structure

m

(2)
4

To simplify the discussion, let us assume that ^t/^x = 0; i.e., t is

a constant independent of x. We also will assume that t is independent of time. This gives

§1 - * ^TM ,

To approximate the spatial derivative, assume that C can be represented

by a quadratic function in x in the region from x^_^ to • Thus, let

o

C = a + bx +cx'' for x . ^ x -s: x m-1 - m+1

From this we get

= b + 2cx; pCj = b + 2cx^

5x \§x/m

a?C = 2c = ih^G )

To obtain the parameter c, we must solve the three equations

Vl = ^ ^^m+1 °^m+1

2

C = a + bx + cx

m m m

C . = a + bx ^ + cx^ . m-1 m-1 m-1

A little manipulation of these equations yields

G = ^(C , - 2C + C ^J/(6x)^ m-1 m m+1 " ' '

where 6x is the spacing between the grid points on the x axis. This

gives Us an approKirnyte expression for |^^^/^-^^jui terms of the values

of C at the thi-ee f':rid points m-1, in, and mM.

j^^C \ ^ -A (C -2C + C U)

From (3) and (4.) wo have

j^C) ^ ^

(C^ ^ - 2C^ + C^.J + S (5)

\§t!m m-1 m m+r m '

'
5

Equation (5) gives us an approximate value of which could be

Used to estimate C from C : m m

C «^ C + m m ill" (7)

This is not the best estimate for C^. If varied vrith time as shown

in Fig. 2, one could get a value of closer to the true value by using the average of ^C/^t and

To be general, let us take a weighted average of the two time

derivatives,

< I > = 4Ui ^ HI}:m

Our estimate for C will now be

m

/Eq.(7)

/Eq.(8)

C ^ C + < >6t

m m ^ dt

t t'

Figure 2.

6x m m+1 X 6tS'

m

(10)

We are now confronted with the problem of determining S^, the value

of the source term at t' . For a small change in S^, we have

S' f« S +- ^^m(C' - C, ) +

m m — Im Im"^

Im

+ miu' - C,, ) Nm Nm

Nm

where the first subscript on C is a species index. The difficulty here

is that we don't have values for the C until an equation like (lO) jm ^ is evaluated for each species j. To get around this problem, we assumie

that all of the partial derivatives ^^jn/^^jm negligible except

in the case of the particular species under consideration in (lO).
6

We then have, (dropping the species index),

S' « S + ^^ni(C' - C ) (11)

m m m m

m

This is obviously not a good assumption. Nevertheless, we need only

a rough estimate of at t' . As t increases, we eventually reach a

steady-state in which C = C . Therefore S also stops changing with '' mm m

^ time. One could probably use the value of S at time t for S' and still

reach the same steady-state. Using (11 ) to estimate S^, however, allows

one to take somewhat larger time steps and consequently speeds up the

integration.

Inserting (ll) into (10) and rearranging, we get,

Z = 1 + + (12) m m m-1 m m+1 m

where

K = A(0/DJX)

A(\) = X5tT/(6x)^

Bix) = (1 - A)6tT(C + C ) + [l - 2(1 - X)6tT - X6t + S 6t

m -—7 m-1 m+1 L T^T ::ip~'J ni m

6x'^ 6x2 ac^j

\{X) = 1 + 2A6tT - \6t ""jn

We are free to choose any value between 0 and 1 for the weighting parametei-. If wo take X ~ 0, then

C' = (0)/D (0) •• ^(C , - 2C + G + S 6t

m m " m ^^'^ '"-1 ni m+1 m

This is an example of a so-called explicit difference equation, in which

the value of C is calculated from the concentrations at time t only, m

When X ^ 0 ve have an implicit difference equation. This gives an

equation connecting three grid points at time t' . This results in a set
7

of simultaneous algebraic equations (l2) which must be solved. This set

has the form

+ AC'

1 2

''2

A c

3 2

2 3

C

3

-B.

= -B,

= -B,

0

A c - d' +A c = -R

M-1 M-2 M-1 M

where M is the maximun number of grid points along the x axis.

For this particular system involving three grid points, the solution of

the above system of equations can be obtained without having to invert

a matrix. Equation (12) can be written in the simpler form,

C = A C' , + B

m m m+1 m

(13)

where

^•1 = m

B

m

'm

Am ^m-1 Bm

(U)

= A

1

B* - A,C' + B,

1 10 1

The values of at the first and last grid points and are always

specified by the boundary conditions for the problem. We can solve (12)

by first calculating A^ and B* for all the grid points starting from

m = 1. It is then easy to obtain the values successively with (13)

starting from (^_^ and working down to CJj

.
3

This implicit method is what we have used to solve the flame

equations. The actual form of the difference equations for the species

and enthalpy used in the program is derived an Appendix H. While the

methods used there to approximate the spatial derivatives are somewhat

different from those used in this section, the basic solution method

is the same.
9

III. DETAILED DESCRIPTION OF THE PROGRAM

The following is a step-by-step description of the main program

and each of its subroutines. Flow diagrams for most of the routines

are provided at the end of each section.

A. Description of the Main Program SPALD/3

The function of SPALD is to control the calculation. It uses

various subroutines to calculate the parameters contained in the finitedifference coefficients in Eq. A40, Appendix H, p. A18. After calculating the coefficients, it then calls a routine CALC which solves the

difference equations and yields values of the concentrations and

enthalpy a time t + 6t. It continues for as many time steps as the

User specifies. No test is made for convergence.

SPALD begins by reading the following data:

1) RIJNID(M), M=1,12

Format(l2A6)

This is the run identification. It consists of 12 words, in A6

format, making 72 characters available for whatever identifying

remarks the user wishes to employ.

2) HEADNG(M), M=1,126

Format(63Al)

This provides headings for the output profiles. There are 126

characters available. Two cards must be used. The output profiles

are listed in 14- columns. The first 12 are for species concentrations

and are in E9.5 format. Column 13 contains the enthalpy profile and

is in El 0.5 format; this extra space is used for the sign of the enthalpy

which can be negative. The last column is for the temperature profile

anjhas E9.5 format. The headings should be spaced within the 126

character range available so that they are aligned at the top of the

appropriate column.

3) LH,L0,K)H,LH02,LH202,LX,LHX,LK2,L02,LH20,LH2,LN2

Forma t(UI2)
10

These are species indices for the H^-O^-N^ flame. They should

run from 1 to 12. In earlier versions of the program they were used

in the subroutine which calculated the chemical source term. Now,

however, they are needed only in a portion of the routine which calculates

transport properties. They must be included although they are used only

for the H^-O^ flame.

U) LD,LE,LC

Format (312)

These three numbers determine the number of integration steps for

which certain subroutines are called at each step. For example, if

LD = 20, then TRANS(l ) is called for the first 20 integration steps.

This means that a complete ualculation of the transport properties is

made at each of the first 20 steps. The second number controls the

calling of SCHM1 , which calculates the part of the chemical source term

arising from the diffusion of the mean molecular weight. SCHMI is

called at each integration step for the first LE steps. The last niimber

controls the reference to SENTP, the routine which calculates the

enthalpy source term. This will be called at every step for the first

LG steps.

5) IDELA,LDELB,LDELC

Format (312)

These numbers control the number of time steps for which TRANS(l),

SGHMI, and SENIP are not called after the first LD, IE, and LC steps,

respectively. For example, after calculating the transport properties

for the first LD steps, the program, from this point on, will calculate

them only at intervals of LDELA steps.

There are two special cases which can arise;

a) Suppose one does not wish to call SCHMI at any time during

the calculation. In most flames, the mean molecular weight changes

little throughout the flame so that this part of the chemical source

term can be neglected. To do this, let LE = 0 and LDELB =0.

b) If one is considering a constant enthalpy flame, SENTP need not

be called. To prevent this, let LC = 0 and LDELC = 0.
n

6) N,NII,NEE,OMR

Forma t( 312, D5. 1

)

N is the number of slices into which the u; axis from 0 to 1 is

divided into to form a grid. It can have a maximum value of /^.0. The

g-rid points start atl = 2 and run to N + 2. (See Appendix H, p. A19.)

Slices between the grid points Nil and NEE have equal widths. There

is a non-uniform grid spacing near the boundaries with the spacing

increasing from I = Nil to 2 and from NEE to N + 2 so that consecutive

slices are in the fixed ratio OMR. This non-uniform spacing is used

to increase the number of grid points in the region which has the

largest concentration gradients. The spacing of the central grid

points is given by the formiila,

[ qM^^^KNII - 1) ^- OMR-'^^KN - nee + 3) - 2.^'-0MR ^ 1

"1

L OMR - 1 + NEE - Nil

J

An example of this grid spacing is shown in Fig. 3 for the case,

N = 15, Nil = NEE = U, and OMR = 1.$.

1=22 Nil NEE 17

i . i . . .

'

0)

Figure 3

Sample Grid Spacing

The grid points Nil and NEE also denote the positions near the

hot and cold boundaries, respectively, used in the grid width control

procedure. They correspond to the points HH and CO (see Appendix J)

.

7 ) JH, JRAD, JM, JMR, JT , JTR, ITER, ITERI , JENTRN, JBODYA, JBODIB

Format(l1l2)

JH is the index of the enthalpy variable F(JH,l). It can have

a maximum value of 13. It must be larger by one than the total number

of species in the flame. The index of the last species will then be

JhMl = JH - 1 . In the output, the number of the column will correspond
12

to the first index in F(J,l). The enthalpy profile, however, will

always appear in column 13.

JRAD is the index of the last minor species. Minor species are

always given the lowest index values. Major species start with the

index JRADPl = JRAD + 1

.

JM and JMR are indices for the mean molecular weight and its

reciprocal which are contained in the array FS. The mean moleciilar

weight is FS(JTyi,l) and its reciprocal is FS(JMl,l). We use the values

JM = 1 , and JMR = 2.

JT and JTR are indices for the temperature and its reciprocal

which are contained in the array FS. The temperature is F3(JT,l)

and its reciprocal is FS(JTR,l). ¥e use JT = 3 and JTR = ^.

ITER is the number of iterations used in the Newton-Raphson •

procedure to calculate the temperature from the enthalpy value at a grid point. ITER = 2 has been found to give satisfactory results.

(See Appendix M, p. A39.)

ITER1 is not used. Set it e(|ual to zero.

JENTRN is the index of the species which is used to calculate the

entrainment rate. A major species is always used for this purpose.

(See Appendix J, p. A27.)

JBODYA and JBODYB are indices for third body concentrations which

are stored in the array FS. FS(JBODYA,l) is the third body concentration

to be Used for the reaction H 4- OH + M ii H^O + M, and F3( JBODYB, I)

that for the reaction H + 0^ + M if HO^ + M in the H^-O^ flame calculations. These concentrations depend on the mixture composition and . are calculated in TRANS. Normally one uses the reciprocal of the mean

molecular weight for the third body concentration in the mole/kg concentration units. We let JBODYA = 5 and JBODYB = 6.

8) XI(J), J=JRADP1 , JHM1

XE( J), J=JRADPl ,JM1

Form'at(5D7.3)
13

These are the mole fractions of all the major species at the hot

and cold boundaries, respectively, at the start of the calculation.

For example, consider a H^-O^-N^ flame having on the cold side

XE(L02) = 0.063, XE(LH20) = 0.0, XE(LH2) = 0.700, and XE(LM2) = 0.237.

Assuming that the reaction goes to completion, we get hot side values

XI(L02) = 0,0

Xl(lJi20) = 2.^>XE(L02)/(1 . - XE(L02)) = 0.13^5

XI(LH2) = (XE(LH2) - 2. ^XE( L02) )/(l . -XE(K)2)) = 0.6126

XI(LN2) = XE(LN2)/(1. - XE(L02)) = 0.2529

9) IPRFL(J), J^JMDPI, JHM1

Format ($11)

If IPRFL(J) - 0, the initial profile for the major species J

decays from the cold to the hot side of the flame. If IPRFL(J) = 1,

then we have a growth profile. For the above H^-O^-N^ example,

IPRFL(L02) = 0

IPRFL(LH20) =1

IPRFL(LH2) = 0

IPRFL(LN2) = 1

10) TOOLD, PRESS, YWIDIH

Format(D6.0,D6.2,D5.1

)

TGOLD is the temperature of the cold boundary in degrees K.

PRESS is the pressure of the flame in atmospheres. YWIDIH is the

approximate width of the flame front in meters. The program needs a value for this qu-antity to start the calculation. It uses this to

calculate a starting value for the grid width parameter PEI = T] - ^q~^

(See Appendix F, p. All.) Although the grid width is automatically

adjusted at each step in the calculation, it is best to start with a

'

value as close as possible to the final value. We have not experiment

ed with the effect of the choice of YWIDTH (and thus the initial PEl)

on the course of the calculation.
14

1 1 ) HaEF( J )

, DGP( J )

, ECP( J ), FCP( j) , Wr( J) , J^l , JHMl

Forraat(/^D12.0,D10.2)

HREF(J) is the reference enthalpy value h? for the species J

in jo-oles/kg. (See Appendix M, Eqs. (A78a) and.(A78b), p. A37 for

their definition.) DCP(J), ECP(J), and FCP(J) are the parameters

d^. , e^. , and f^ in the quadratic expression used for the heat capacities

of the individiJal species. The units are joules/kg-deg. (See Appendix

M, Eq. (A76), p. A37.)

It is worth noting here that one needs values of these parameters

only for the major species. The enthalpy, h =X^.h,Y., so that contriJ J J butions from species having small Y. are not important in its calculation.

The h. (actually h'i" - M.h.) also appear in the source term for the

J J J J enthalpy as the products h*^cp./Scju. Here again the contributions from

J J the trace species will be negligible compared to that from the major

sioecies. This is so because the values of Scp./^ut) will be small when

J cp. is small even though these gradients could be large on a percentage

J basis. Calculation time could be saved by not summing over the trace

species in the enthalpy source term. The effect of neglecting these

terms should be examined in future work.

WT(J) is the molecular weight of each species in kg/mole.

12) DMN,DXMAX, STEPS

Format(2D5.0,D6.0)

The size of the time step DX is determined from the formiila

DX = DXMIN + ((DXKAX - DXMIN ) /STEPS )^'-L, where L is the number of the

step and STEPS is a number comparable to the maximum number of steps

to be Used. The step size increases linearly from DXMIN to DXMAX when

L - STEPS. We normally have DXMIN « DXMAX. Thus DX is small at the

beginning of the calculation and increases as the steady-state is

approached. Other step size control formulas could be investigated.

Perhaps the step size could be made inversely proportional to the rate

of change of the grid width PEI. As the steady-state is reached, PEI

changes little from step-to-step justifying larger values of DX. At

the beginning of the calculation, when the profiles are changing rapidly.
15

PEI will also be changing and the step size would be kept small.

13) LMAXjLPRINT

Format (213)

IMAX is the maximum number of steps to be used in a given calculation. LPRINT specifies the step interval at which profiles are to be

printed. One page of profiles is printed every LPRINT steps. When

L reaches IMAX, the complete output is printed.

U) ALPHA

Format(F3.2)

ALPHA is the parameter v which appears in the grid control formulas

Eqs. A58-60. (See Appendix J, p. A26.) We have been setting ALPHA

equal to 0.1 . ,

15) KRAT

Forma t( 12)

This is the number of chemical reactions. For a particular

reaction, the foCTfjard and reverse reactions count as one reaction.

(Note that MATE = KRAT. ) KRAT can have a maximum value of 30.

16) FRQ(K),BETA(K),EAGT(K),BFLAG(K),EFLAG(K),K=1,KRATE

Format(D6.2,D7.2,D8.2,2Ll)

These are the Arrehenius parameters for the forward rate constant

of reaction K. If we express the temperature dependence of a rate

constant k as k = AT^ exp ( -E/RT ) then

FRQ(K) ^ A

BETA(K) =3

EACT(K) = E/R (in °K)

k must be in units of m /mole and m /mole for two and three body reactions, respectively.

BFLAG(K) = .TRUE, if BEIA(K) ^ 0, and EFLAG(K) = .TRUE, if

EACT(K) ^ 0.
16

17) EQA(K),E(3B(K),EQC(K),B0D1(K),B0D2(K),K=1,KRATE

Format(3D11 .5,2L1

)

The first three quantities are the coefficients in the expansion

of the free energy of each reaction as a power series in the

temperature. We use

-AF° =EQA(K) + EQB(K)*T + EQC(K)^^T^

The values of these coefficients must be such that AF° will be in units

of kcal/mole.

If B0D1(K) = .TRUE., then either the forward or reverse reaction

is third order. If B0D2(K) = .TRUE., then the reverse reaction is

third order.

18) ITEST(M),M=1,10

Format (1011)

When set to unity these flags produce various outputs used for

checking the operation of the program. They are normally set to zero.

An examination of the listing of SPALD will show what quantities are

printed. They should be used only for flames with JH ^ 6. With some

modification of SPALD they could be used for larger values of JH.

19) INPITT

Format (12)

When INPUT = 1, input profiles are supplied from a data file.

This file must contain the concentration and enthalpy profiles and the

grid width PEL This file is taken from logical unit 4-5 by the statement KEABU5) ((F(J,I),J=1,JH),I=1,NP3),PEI. These can be output

profiles from a previous calculation. The output profiles are placed

in logical unit by the statement WRITE(4.6) ( (F( J,I) , J=1 , JH) ,1=1 ,NP3),

PEI. (NP3 equals N + 3) If the user is generating his own input profile

they must be stored in the above order and be double precision numbers.

Note that F(J,l) and F(J,NP3) are not actually used and should be given

zero values. Remember that F(JH,l) is the enthalpy, and that small

values of I refer to the hot side of the flame.

When INPUT = 0, the initial profiles are generated by INITL.
17

20) L0RDER(J),J=1,JHM2

Format(UI2)

These numbers specify the order in which one wishes the species

equations to be solved. Spalding has found that the efficiency of the

calculation is improved by solving the species equations in order of

increasing species concentration. The trace species equations are

solved first, followed by the major species equations. (Note that

JHK2 = JH - 2; this is the index of the next to the last major species.)

The species JHM1 is calculated from the relation^j:p^.M^. = 1 , and so

LORDER(JHMl) need not be specified.

As an example, consider the H^-Br^, flame (JH = 6) where the species

have been given the indices Br = 1 , H = 2, -- 3, Br^ = and HBr = 5.

If we wanted to solve the species equations in the order H, Br, Br^,

H^, HBr, we would have

LORDERd ) = 2

L0RDKH(2) 1 L0RDER(3) =

LORDER(/i) = 3

The equation for J = 5 will automatically be solved last.

(Note that in the test case to be presented later for which this flame

was Used, we solved the species equations in the order of species index,

i.e., we used K)RDER(J) = 1, 2, 3, ^.)

21) SPECIE ( J ),J=1,JM2

Format (UL1)

If SPECIR(J) = .FALSE., then the equation for that J is not solved.

22) REACT (K),K=1,KRATE , Format(30Ll)

If REACT(K) = .FALSE., then both the forward and reverse rate

constants for reaction K are set to zero.

23) FFLAG(K) ,K^1 ,KRATE

RFLAG( K) , K^l , KRATE

Forma t(30Ll)
18

If FFLAG(K) = .FALSE., then the forward rate constant for reaction

K is set to zero; if RFLAG(K) = .FALSE., then the reverse rate constant

is set to zero.

2U) LA(K),1£(K),LRA(K),IJIB(K),K=1,KRATE

Format(4.l2)

These numbers specify the reaction mechanism to be used. Consider

the H^-Br^ flame as an example. We use the reactions,

K - 1 Br^ + M il Br + Br + M

K=2 H^+MftH + H+ M

K = 3 Br + ?i HBr + H

K = ^ H + Br^ ^ HBr + Br

For a species index assignment Br = 1 , H = 2, = 3, Br^ = U, HBr = 5,

we have the following assignment for LA, LB, LRA, and LRB:

K LA(K) LB(K) UIA(K) LRB(K)

1 U 0 1 1

2 3 0 2 2

3 1 3 5 2

2 U 5 1

We see that the indices LA(K) and LB(K) refer to the species entering

into the forward reaction, while LRA(K) and LRB(K) pertain to those

contained in the reverse reaction. There is space for two species on

each side of a reaction. The third body M is not counted as a species.

In cases where there is only one species aside from M involved in the

reaction, one must . always set either LB(K) or LRB(K) equal to zero. The

index, of the one species present should be assigned to LA(K) or LRA(K).

For example, if we had written reaction 1 as

Br + Br + M ^» Br^ + M

then the assignment would have been,

LA(i) =r 1, IB(1) = 1, LRA(l) = 4., LRB(i) = 0.
19

25) FLAG(I),I=1,20

Forma t( 30 L1)

These flags provide a series of options regarding the calculation

of transport properties, source terms, and the nature of the output.

FLAG(1 ) = .FALSE. The complete calculation of the transport parameters

is made.

= .TRUE. Constant and equal transport parameters are used.

These are contained in the array PREF and are assigned values in a DATA

statement in SPALD, line SPD00350. Constant, but non-equal values could

be used by changing this data statement.

FUG(2) = .FALSE. If the flame contains a buffer gas like in a low

temperature H^-O^-N^ flame, its mass fraction can be considered to remain

constant throughout the flame. This buffer will be assigned the species

index JHM1 . The relation between the mass fractions gives

WHM3y ^ ^J=1 J JHM2 JHM1

When Yji^.^ is the buffer mass fraction and is kept constant,

be calculated from this equation and its species equation need not be

solved. Therefore when this option is used one must also have SPECIE(JHM2)

= .FALSE..

= .TRUE. The buffer concentration is allowed to vary. This

option must be used for flames containing no buffer species.

FIAG(3) - .FALSE. The pressure equals 1 atmosphere.

= .TRUE. The pressure is different from 1 atmosphere.

FIAQU) (Not used)

FIAG(5) = .FALSE. The enthalpy source term is not calculated. This option

can also be achieved by using LC = 0 and LDELC = 0.

= .TRUE. The enthalpy source term is calculated.

FLAG(6) (Not used)
20

1i'LAG(7) = .FALSK. The transpoi-t propertied for the ozone decomposition

flame are not calculated.

= .TRUE. The transport properties for the ozone flame are

calculated.

FLAG(8) = .FALSE. The transport properties for the hydrogen-bromine

flame are not calculated.

= .TRUE. The transport properties for the hydrogen-bromine

flame are calculated.

FLAG(9) = .FALSE. The rate constants as a function of distance through

the flame are not printed in the final output.

- .TRUE. The rate constants are printed.

FLAG(IO) = .FALSE. The rates of each reaction at each point in the

flame are not printed in the final output.

= .TRUE. The reaction rates are printed.

FLAG(11) = .FALSE. Heat release rates for individual reactions at each

point in the flame are not printed in the final output.

= .TRUE. Heat release rates for the reactions are printed.

FLAG(12 through 20 ) (Not used)

This completes the input data requirements for SPALD.

After reading this data SPALD calls the subroutine S1«TGH. Its

function is to set three flags, G1 , G2, and G3, which describe the

nature of each reaction. These flags are then used in the calculation

of the chemical source term by subroutine SCHM2.

FUG(l) is then tested. If false the routine DIFUSE is called.

This calculates binary diffusion coefficients for the ^2'^2~^2 ^^^'^^ from Lennard -Jones parameters. The values of these parameters are

contained in a DATA statement in DIFUSE.
21

The index JSTOP is set equal to JHM3. FLAG(2) is then tested;

if true, then the b^offer gas will vary and JSTOP is equated to JHM2.

The flag INPUT Is tested. If INPUT = 1, profiles are read from

a file located in logical unit 4-5.

The quantities L, LCOUMT, X, and DX are initialized. L is the

number of the step, LCOUNT is the number of steps taken since the

last profile print-out, X is the time variable, and DX is the time

step size.

The grid structure is calculated by a call to OMEGA.

The flag INPUT is again tested and if equal to zero INITL is

called to calculate the initial profiles.

After these preliminaries, the point 706 is reached to which the

program returns after every integration step. This therefore represents

the beginning of each step. TRANS(l) is then called if a complete

transport calculation is to be made, otherwise TRANS(O) is called.

TRANS(O) calculates the temperature at each grid point from the enthalpy and stores it in FS(JT,l). TRANS(l ) in addition calculates the

transport parameter array PREF(J,l).

The flag INPUT is tested once again and, if zero, an initial value

of PEI is determined by calling ZCALC with PEI = 1.0. This routine

produces a value of ^(l), the laboratory spatial variable, at the cold

boundary grid point NP2. PEI is then taken to be the ratio YWIDrH/Y(NP2)

.

To see the rationale behind this, refer to Eq. (A70), Appendix L. This

formula gives for Y(NP2)

EN+2

i-2 p . -

Pi

The quantities oj. wore calculated in OMEGA, and the p. were calculated

in TRANS(l). We want to have a value of 'I] (^PEI) such that Y(NP2) YWIDTH.

By calling ZCALC with PEI = 1.0, we get the quantity

vN+2

2 ^+1

" ^i

i=

Dividing YWIDTH by this gives a value of PEI which makes Y(NP2) = YWIDTH.

i-2 P.+^ - p.
22

SPALD next calculates the quantities PEID2 ^ ^T], PEIOMD(l) =

Tl(ci)^^^ - 0).), and PEI0M2(I) = " \_^) • RATGN is then called to calculate the values of the rate constants

at each point in the flame.

To calculate the entrainment rates and (see Appendix J, p. A26)

,

it is necessary to have a value of R( JE.NTRM, l) ; this is the produ.ction

rate of species JENTRN in .Tioles/kg-s . Since this is normally calculated

by SCHM2 later in the program, on the first integration step it is.

necessary to call SC1M2( JENTRN) to calculate the production rate for

this species.

SPALD then evaluates the finite -difference coefficients A(J,l),

B(J,I), and D(J,l) (minus the source term). These are the quantities

A', B', and D + 23^ of Eq. (A40), Appendix H, p. A18.

If desired, the chemical source term arising from the diffusion of

the mean molecular wieght is then calciolated by a call to SCHM1

.

The finite-difference coefficient C(J,l) (minus the source term)

is then calculated. This is the quantity C - 2Sp of Eq. (A4.0).

The chemical source terms Sp and S^ are then calculated for each

species by a call to SCM2(J). Also calculated by this routine are

the production rates R(J,l). Following this, C(J,l) = C and D(J,l) =

D are calculated.

New species concentrations for time t + 6t (i.e., X ^ DX) are then

calculated by calling CALC(J) in the order specified by the numbers

LORDER(J).

Next, the hot boundary condition {^cp/<^uo[|^ = 0 is set by equating

the concentration at the hot boundary point to the value at the adjacent

grid point; i.e., we set F(J,2) = F(J,3).

The value of the concentration of the JMI species is then calculated

from the relation^.cp .M . = l ; if the buffer concentration is being

J J J kept constant, then the concentration of the JHM2 species is calculated

from this relation.

If the enthalpy is allowed to vary, its source term is calculated

by a call to SENTP and its new value at t + 6t is calculated by calling

CALC(JH).
23

The step counters L and LCOUNT are then incremented. If an output

point has been reached, flame velocities are calculated for each species

by a call to ENTRN(J,l). These velocities should be the same for each

species. In practice, this will happen only for the major species

(see Eq. (A25) Appendix F, p. A12). The laboratory distances Y(l) are

then calculated by calling ZCALC and the partial results are printed

out by a call to OUIPUT. When L = LMAX, the call to OUTPUT produces

a complete print out. In addition, the profiles for concentrations

and enthalpy are stored in logical unit 4-6.

If L is not equal to IMAX, then the program returns to point 706

and the calculations are repeated.
24

SPALD/3

READ DATA

JST0P^JHM3

JST0PtJM2 3

l^AD INPUr PROFILES

I

"

L=0

LcouFr=o

X-0.0

DX^DXMIN

CALL INPRFL
25

SPALD/3 (cont.)

J »—^CALL TRANS (1)^'

LTABA^LTABA+1

LTABA=0

PEI=1 .

0

( CALL ZCALC

PEI^YWIDTH/Y(NP2)

CALGiaATFJ:

PEID?

PEIOMD(I)

PEI0M2(I)

CALL RATCN
26

CALL ENTM( JENTfm,0)

CALCUUTE:

A(J,I)

B(J,I)

D(J,l) minus

source term

CALCUIJITE:

C(J,I) minus

chemical

source term
27

CALCULATE:

C(J,I)

D(J,I)

r^LL CALC(J)J

F(J,2)=F(J,3)

SPALD/3 (cont.)

T CALCULATE:

F(JHM1 ,1)

buffer varies

CALCULATE:

F(JHM2,I)

buffer const.
28

SPALD/3 (cont.)

^CALL CALC(JH)^

L=L+1

LC0UNT=LC0UNT+1

DX=DXMIN+FAC2^L

X=X+DX
29

B. Description of Subroutine SWITCH( KEIATE)

This subroutine sets three flags G1 , G2, and G3 which denote the

presence and location of a particular species J in a particular reaction K. These flags are initially set to the .FALSE, position. SWITCH

sets G3(K,J) to .TRUE, if the species J appears twice on either side of

the reaction equation. If it appears on the left side of the equation,

Gl(K,J) is set to .TRUE. J if on the right side, both Gl(K,j) and G2(J,K)

are made .TRUE..

As an example, consider the following reactions which are tested

for the presence and location of species A.

Values of the Flags after Calling SWITCH

Reaction G1 G2 G3

A + B C + D .TRUE. .FALSE. .FALSE.

C + D A + B .TRUE. .TRUE. .FALSE.

A + A C ^ D .TRUE. .FALSE. .TRUE.

C + D A + A .TRUE. .TRUE. .TRUE.

B. + D, E -f- F .FALSE. .FALSE. .FALSE.

These flag values are used in SCHIyl2 in the calculation of the

chemical source terms.
30

t

SUBROUTINE SWITCH(KRATE)

DO: J^1,JHM2

K=1 , KRATE

G3(K, J) = .TRUE.

Gl (K, J)=.TRtJE.

G3(K, J) = .TRUE.

G1 (K, J)^.TRUE.

G2(K,J) = .TRUE.

ZZD
31

C. Degcription of Subroutine DIFU3S(PRBS3)

This routine calculates the temperature independent portion of

the binary diffusion coefficients for the H^-O^-N^ flame from LennardJones parameters. These parameters are assigned values in the DATA

statement starting at line DIF00100. M is an array containing the

molecular weights in gm/mole; EK is the Lennard-Jones potential

o

° parameter e ./k in K; and S is the parameter a. in A. These are listed

J J in order of species index. The quantities calculated are DA(I,J) and

DB(I,J) which are the temperature independent portion of Eqs. (A72) and

(A73), Appendix M, p. A36. Also calculated is the cut-off temperature

below which DA(I,J) is used rather than DB(l,J). This is TDIV(I,J) =

3e^^./k. I and J are species indices. The index I runs from 1 to JHM2,

while J goes from JRADPl to JHMl . This means that DIFUSE calculates

these quantities for interactions between minor and major species,

major and major, but not minor and minor.

If FLAG(3) = .TRUE., then DIFUSE divides DA and D3 by the pressure

PRESS.
32

SUBROUTINE DIFUSE(PRESS)

DO: I=1,JHM2

J-JRADPI , JHM1

CALCULATE:

(e../k)°-"

2

a. .

T C Ur=3e

^ ^./k=TDI V(I , j)

FAC^( (M. ^^4 . )/M.M . ) V( Pa? .

)

DA(l,j)=1.26x10"'^(FAC) X (e.^./k)"°-^

DB(l,J)=1.66x10"'^(FAC) X (e../k)

^ J

-0.17

DA(I,J)=DA(I, J)/PRESS

DB(I,J)=DB(I, J)/PRESS

I
33

D. Description of Subroutine OMEGA

This routine takes the values of N, Nil, NEE, and OMR and calculates the grid structure. (See Sectionlll-A, p.n.) The quantities

determined are

0M(I) = (D.1

OMD(I) = iu.^^ - CO. 1+1 1

ROMD(I) = 1./0MD(I)

OMP(I) = u; . + OD.^^

BOM(I) = 0).^^ - aj._^

B0MT3(I) = 3.'^B0M(I)

Note that lines OMGOO3OO and 310 can be removed since the variables

OMI and OME which are calculated there are not used.
34

E. Description of Subroutine INITL/I

This program begins by defining the function lOu;^ - 15u/^ + 6cd^

which is used to generate S-shaped profiles. It then tests the variable

INPRFL(J), where J is the index of a major species. If .FALSE., then

the initial profile for J is a decay profile and is calculated from the

formula

XI(J) + (XE(J) - XI(J))(lOaa^ - + 6(ju^)

If .TRUE., then J has a growth profile which is calculated by interchanging XI and XE and inserting 1 - cju for ou in the above formula. This

yields mole fraction profiles. To convert them to units of moles/kg,

INITL calculates the average molecular weight at a particular grid

point from the formula (M> =^.x.M.. The desired profiles are gotten

J J «] by dividing the mole fraction profiles by (M).

At this point a small subroutine INPRFL is called to set the initial

profiles for the trace species. Normally, the trace species are assumed to have ngeligible concentrations at the cold boundary grid points.

In setting up INPRFL, do not use equal concentrations at the two boundary grid points, 2 and NP2. This would result in a division by zero

on the first call to ENTRN. (See ENTRN, lines ENT00190, 210, and 230.)

INITL next calculates ET(J,JP2), the enthalpy at the cold boundary

grid point for all the major species. This will be in units of joules/kg.

From these quantities and the composition at the cold boiJndary, the

total enthalpy value h^ = F(JH,NP2) is calculated. The enthalpy at the

hot boundary, F(JH,2), is then set equal to h^-, and the temperature

there is calculated from this and the initial hot side composition by

the Newton-Raphson method discussed in Appendix N, p. A39. From the

two temperatures TCOLD and THOT an S-shaped temperature profile is

calculated. Finally, the enthalpy at each grid point is given the cold

boundary grid point value.
35

DEFINE: i ^ FNC(U)=iolr-1$lr+6Ir

I

DO: J- JRADP1 , JHM1

DIFF^XE(J)-XI(J)

DO: 1=2, NP2

U-OM(I)

F( J, I ) =XI ( J) +DIFF^^FNC( U)

..£77TTTI

Convert F(J,l) to mol/kg

DO: 1=2, NP2

AVM^

JHM1

JRADP1

F(J,l)^^Wr(J)

F(J,I)=F(J,I)/AVM

f*GALL INPRFL

J

SUBROUTINE INITL/i

_ ^

DIFF=XI(J)-XE(J)

DO: 1=2, NP2

U=1 .0-OM(I)

F ( J , I ) =XE ( J ) +DIFF^'^ FNC ( U)
36

T=TCOLD

INITL/1 (cont.)

CALCUUTE: ^ ET ( J , NP2 ) =h° +d . T+i-e . T'^+if

J J J tJ

F ( JH, NP2 ) = ^

J .TMJjp^ WT ( J ) ^^ET ( J , NP2 ) ^ F( J , NP2

)

F(JH,2)=F(JH,NP2)

A=

JHM1

:JRADP1 h°WT(j)^^F(J,2)

J

BJHM1

:JRADP1 d.WT(j)^^F(J,2)

J

C=E.JHM1

:JRADP1

e.l/fr(j)^^F(J,2)

JHM1

=JRADP1

f .WT(j)^fF(J,2)

T:=(F(JH,2)-A)/B

ICOUNT-0

I

T ^T^( JH, 2) -( A+B'iM-J-, GT"^ f^DT^]l/(B+GT fDT'^)

.2>

I

IG0UNT=IC0UNT+1

THOT-T
37

INITL/I (cont.)

DO: 1=2, NP2

.0-OM(I)

FS( JT , I ) = ( THOT-TCOID) ^^FNC( U) H-TCOLD

DO: 1=3, NP2

F(JH,I)=F(JH,2)
38

F. Description of Subroutine TRANS( ICALC)/6

This subroutine calcu3.ates the temperature, the mean molecular

weight, the density, and the transport parameters for three flames;

ozone decomposition, H^-Br^j and H^-O^-N^. It begins by calculating

the temperature and its reciprocal from the enthalpy by the NewtonRaphson procedure described in Appendix N, p. A39. At the same time

it also calculates the mean molecular weight and its reciprocal at

each grid point.

The quantity ICALC is then tested; if it equals zero, then TRANS

jumps to the next grid point. Under this condition the transport

properties are not calculated. If ICALC equals one, then TRANS continues

by calculating the density RHO(l) at each grid point. FLAG(3) is tested,

and if .TRUE., the density is multiplied by PRESS which, in this case,

will be different from unity; if .FALSE., then this multiplication is

not made.

FLAG(7) is tested and if .TRUE., transport properties are calculat2 ed for the ozone decomposition flame. The desired parameters are A.p

for diffusion, and ^P/^^ for thermal conductivity. The average heat

capacity C (=CPMIX) is first calculated and from this and the density,

P -3 i

Xp/G (=PREF( JH,I)) is calculated from the expression 2.2026x10 ^T^p/C . P P

For this flame a Lewis number of 0.94- was used so that the quantities

A.p^ (=PREF(J,I)) are all given the values 0.9/;^^PREF( JH, l) . The program

then jumps to statement 57.

If FLAG(7) was .FALSE., FLAG(8) is tested; if .TRUE., transport

properties are calculated for the H^-Br^ flame. We used the following

expressions to calculate the diffusion parameters;

Ayp^ ^ 1 .05x10"^(T/TC0LD)'' •^'^p^

P^' ^HBrP^ " 0.101xlO"^(T/TCOLD)''-^^p^

The thermal conductivity of the mixture was taken to be

XP/Op ^ 3.349Ux10"^(T/TC0LD)°*^'^p/CPMIX.

The program then jumps to statement 57.
39

If FLAG(8) was .FALSE., FLAG(l) is tested; if .TRUE., the program

jumps to the next grid point. This is the situation where we are using

constant values for the transport parameters. If FLAG(l ) = .FALSE.,

then TRANS proceeds to calculate the transport properties for the

H^-O^-N^ flame. It begins by calculating the binary diffusion coefficients using the values of DA(I,J) and DB(I,J) determined in DIFUSE.

1 .9/ Since each species at a particular grid point uses either T * or

T

* as the factor multiplying DA or DB, it is desirable to avoid

doing these exponentiations more than once per grid point. This feat

is accoaiplished by assigning .TRUE, values to the variables TAGA and

TAGS (initially .FALSE.) if either t'' or T^*^'^ has already been calculated for a particular species. Testing these flags on subsequent

passes through the DO loop allows one to avoid repeated evaluation of

T

' or T*

. The temperature dependent binary diffusion coefficients

are stored in the array DD(KI,KJ), where KE and KJ are species indices.

The diffusion coefficients A. are then calculated from the DD(KI,KJ)

J values and the concentrations F(J,l) by means of Eq. (A9), Appendix C,

p. A7. Diffusion parameters PREF(J,l) = A.p are then calculated.

J TRANS next calculates the thermal conductivity of the mixture

from the formula given in Appendix D, p. A8. For this, it is necessary

to have values of the temperature dependent coefficients A^^. defined by

Eq. (AII), Appendix. D, p. AS. These were calculated from this formula

by a separate program, and the results fit by a least squares procedure

to a power series in the temperature. These expressions for the A^^.

are given by lines TRN01510 through 01620. Note that the values are

stored in the DD array since the binary diffusion coefficients previously

stored there are no longer needed. Also given just before these are

the power series expansions of the thermal conductivities of the pure

major species; lines TRJNfOU20 through 0U50. These are placed in CON(J).

Note that it is for these arrays that the indices L02, LH20, LH2, and LM2

are required. From the A . and CON(J) values, the mixture thermal conJ ductivity is calculated. This is then converted to the thermal conductivity parameter PREF(JH,l) = Xp/C^.
40

FLAG(5) is then tested. If .FALSE., TRANS goes to the next grid

point. If .TRUE., the enthalpy source term is being calculated and

values of ENT(J,l), the enthalpy of each species at each grid point

are determined.

After the calculations have been completed for all the grid points,

ICALC is tested. If equal to zero, a return to SPALD occurs. If equal

to one, TRANS proceeds to calculate the values of the reciprocal of the

mean molecular weight at the control volume boundaries. This is taken

to be the average value of FS(JM.R, I) at two adjacent grid points, and

is stored in RMB(l). The same thing is also done for the diffusion

and thermal conductivity parameters, and the results placed in PREF(J,l)

and PREF(JH,I). Thus, these arrays will contain the control volume

boundary values and not the grid point values.

If FLAG(5) = .TRUE., several quantities appearing in the enthalpy

source term are evaluated. Note that both parts of the enthalpy source

term are calculated together and not separately as in the case of the

chemical source terra.
41

SUBROUTINE TRANS(ICALC) /

1

T=(F(JH,I)-A)/B

AB=

JHM1

JHM1

JIMI

d WT(j)-"-F(j,l)

e.Wr(J)^F(J,l)

f .WT(J)^'^F(J,I)

fs(jmr,i)=2];j4^^ f(j,i)

FS(JM,I)=1.0/FS(JMI,I)

ICOUNT-^O

i

T^T+ [(F( JH, I ) -( A+BT+i-CT^+iDT^ )] / (B+CT+DT^)

I

ICOUNT^ICOUNT+1
42

0

TRANS/6 Icont.)

RHO ( I ) =1 . 21 867'^ FS ( JM, I ) ^FS( JTR, I

)

RH0(I)=RH0(I)^;-PRESS

CPMIX=B+^CT+iDT'

PREF( JH, I ) = ( 2 . 2026x1 0 ^T^ ) ^^RHO (I ) /CPMIX

©

I

DO: J=1,JHM1

PREF( J, I ) ^0 . 9^^=-PREF( JH, I

)

I

FAC^(T/TCOLD)^"M .67

I

DO: J^3,5

, PREF(J,l)--0. 101x10 "FAG

I

PREF( 1, 1 ) :=0 . 1 5 5x1 0'^^FAC

PREF(2,I)^1 .05x10"^FAG

PREF( JH, I ) ^3 . 3^9Ux1 0"^( ( T/TGOLD) ^=-^0 . 67

)

CPMIX-^]3^4GT+iDT^

PliEF ( JH , I ) =PREF ( JH , I ) ^^RHO ( I )/CPMIX

R2DM=RH0 ( I ) ^^RHO( I ) ^^FS( JMl, I

)

DO: J=1,5

PREF(J,I)^PREF(J,I)TM0(I)^"^2
43

)

TAGA=. FALSE.

TAGB=. FALSE.

DO: KI=1,JHM2

KJ=JRADP1 , JHM1

DD(KI,KJ)=DA(KI,KJ)^fTA

TRANS/6 (cont.)

TAGA=.TRIIE.

DD( KT , KJ) =DB ( KI , KJ) ^^TB

TB=T-!^^^1 .67

TAGB=.TRUE.

)
44

I

TRANS/6 (cont.)

DO: J=1,JFM1

PREF(J,l)=rO.O

I

DO: K^JRADP1,TM1

DD(J,K)^DD(K, J)

1

PREF(J,I)=PREF(J,I)+F(K,I)/DD(J,K)

-13

PREF( J , I

) =R2DM/PREF ( J , I

)

(b

J > jrad: PREF(J,I) = (1 .0-F(J,l)^^FS(JM,l))^PREF(J,l)
45

TRANS/6 (cont.)

CALCULATE:

Thermal conductivity

of major species in

H^-O^-N^ flame.

Power series in T is

Used.

CON(J)

PREF(JH,I)=0.0

DO: KI=JRADP1 , JHM1

zn

SUM(KI)^0.0

DO: KJ^JRADP1 , JIMI

S IM( KI ) =SLH( KI ) -f-DD( KI , KJ) •^F( KJ, I

)

CALCULATE:

Array A(I,J) for major

species in ^2~^2~^2

flame.

Power series in T is

used.

This array is stored

in DD(J,I)

CPMIX=B+|-CT^DT^

4

SLM(KI)-1 ,0+SUM(KI)/F(KI,l)

PREF( JH, I

) -PREF( JH, I

) +CON( KI )/SUM( Kl)

4

PREF( JH, I ) ^PREF( JH, I ) -"-RHO( I )/CB4IX
46

9

I

I

I r

k I

I I

I L

TRANS/6 (cont.)

€)

ENT ( J , I ) =WT ( J) h?+d . TH^e . T^-fif .

)

J J J J

I

1 DO: 1^2, NP1

I

i

I

RMB ( I ) FS( JMR, I ) +FS( JMR, 1+1 )

)

PREF( JH, I )=!-( PREF( JH, I ) +PREF( JH, 1+1 )

)

zzA

DO: 1-^2, m

J^l , JM1

I

PREF( J , I ) -^k ( FUEYiJ , I ) +PIIEF( J, I +1 )

)

ENTLPZ=^(ENT(J,I)+ENT(J,I+1 ))

HRDM( J , I ) =^ENT LPY^ PREF( J , I ) ^^RMB ( I

)

TRCP( J, I) -ENTLPY^^PREF( JH, I

)

DR2M( J, I ) =PREF( J , I ) ^^RMB ( I

)
47

G. Description of Subroutine RATCN/3

This subroutine takes the Arrhenius parameters and calculates the

forward rate constants for each reaction; from these, and' the free

energies of reaction, it then calculates the rate constants for the

reverse reactions.

RATCN starts by calculating the quantities (RT) at each grid

point. R is the gas constant and has the value 8.2057x10 atm-m /mole-de

These are stored in RGT(l) and are used later in the calculation of the

reverse rate constants.

The calculation of the forward rate constants starts by the testing

of FFLAG(K); if .FALSE., RATCN sets RATE(K,l) to zero and goes to the

next K value. This option allows one to examine the effect of neglecting the forward rate without having to change the Arrhenius parameters

which may have been used in a previous calculation.

Next, BFLAG(K) and EFLAG(K) are tested together; if both are .FALSE,

then both 3 and E/R are zero for this reaction and so exponentiation

can be avoided. In this case RATE(K, l) is given the value FRQ(K), the

Arrhenius frequency factor, and RATCN moves to the next K. If either

BFLAG(K) or EFLAG(K) or both are .TRUE., BFLAG(K) is tested and if

.FALSE., RATE(K,I) is given the value FRQ(K)-'^DEXP(-EACT(K)/FS( JT,l))

.

If BFLAG(K) = .TRUE., then EFLAG( K) is tested and if .FALSE., RATE(K,l)

takes the value FRQ( K)^KFS( JT, l)^^^fDETA(K) ) . If EFLAG(K) = .TRUE., then,

RATE(K,I) is given the value FRQ( K)^^FS( JT, l)^!-^'^riETA( K) )^^DEXP( -EACT( K)/

FS(JT,I)). Only in this last case is it necessary to calculate the complete Arrhenius expression.

After calculating all of the forward rate constants, RATCN evaluates

all of the equilibrium constants. Consider the general reaction,

k

R.12+ Ro+- + R i^P, +P^ + + P

^ k' ^ ^ P

R^ and P^ represent one molecule of a species, so two symbols may

represent the same molecule. The equilibrium constant in terms of the

partial pressures of the various species is

Kp ^n^^l^P. = exp(-AF7RT)

n?=i Pr.

1
48

We want the equilibrium constant in terms of concentrations; this is

related to by the expression = Kp(RT) ^, where n = p - r and R

is the gas constant in atm-m'^/inole-deg. Since K = k/k' , we have for

c the reverse rate constant k' = k(RT) /K . / p

To calculate k' RATCN firsts tests RFLAG(K); if .FALSE., it sets

the reverse rate constant k' (=RA.TER( K, I)) to zero and goes to the next

K value. If this flag is .TRUE., then it calculates -AF°/RT from the

expression -EQA(K)/T - EqiB(K) - EQC(K)^^T, evaluates exp(-AF7RT), and

storesthe results in ECON(l). (Note that EQA(K), EQB(K), and EQC(K)

as read by SPALD give AF° in kcal/mole. SPALD immediately multiplies

them by 1000/R, where R = 1.9869 cal/mole-deg to give AF°/R. ) B0D1(K)

is tested next; if .FALSE., then the number of molecules does not change

in the reaction; n = 0 and RATER(K,l) is given the value RATE(K,l)/

ECON(I). If B0D1(K) is .TRUE., then we have a 3-body reaction with a unity change in the number of molecules. B0D2(K) is tested to find out

the sign of n. If B0D2(K) = .FALSE., then n = -1 and RATER(K,l) is given

the value RATE( K, l)^RGT( I)/EC0N( l) ; if .TRUE., then n = +1 and it has the

value RATE(K,I)/(EC0N(I)-=^RGT(I)). (Remember that RGT(l) = (RT)"''.)
49

SUBROUTINE

RATCN/3

DO: 1=2. NP2 . RGT ( I ) =1 . 21 867x1 ( JT , I

)

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

IL

I

DO: K-=1,KRATE

DO: 1=2, NP2

RATE(K,I)=0.0

DO: 1=2, NP2

RATE(K,I)=FRQ(K)

DO: 1=2, NP2

RATE( K, I ) =FRQ( K) ->^DEXP( -EACT ( K)/FS( JT , I )

)

DO: 1=2, NP2

RATE ( K, I ) =FRQ( K) FS ( JT , I ) ^^^-^BET A( K)

)

DO: 1=2, NP2

RATE( K, I )=FRQ( K) ( FS( JT , l) ->"^BET a( K) ) ^^^DEXP( -EACT ( K)/FS( JT , I )

)

"(5
50

RATCN/3 (cont.)

DO: 1=2, NP2

RATER(K,I)=0.0

DO: 1=2, NP2

ECON(l)=DEXP(-EQA(K)/FS(JT,l)-EQB(K)-EQC(K)^-FS(JT,l))

DO: 1=2, NP2

RATER ( K, I )=RATE( K, I)/(EC0N( I )^HiGT(l )

)
51

H. Description of Subroutine ENTRM( J, IFOW)

The purpose of KNTRN is to calciolate the so-called entrainment

rates m^^ (=RME) and m^^ (=RMl). RME is the mass flow rate in kg/m -s

across the cold boundary and RMI that across the hot boundary. When

a steady-state is attained, RME and RMI become equal. Their values

depend on the concentrations of the species JENTRN at the hot and cold

boundaries and at the points NEE and NIX. The grid width PEI is changed if RME and RMI have different values. By this means the grid can

be kept centered about the region containing the largest concentration

gradients. (See Appendix J, p. A26.) This routine also calculates

the flame velocity.

During the course of the integration ENTRN is called with J =

JENTRN and IFOW = 0. The routine begins by evaluating the integral

where R. is the mass production rate of species j in kg/m -s. The

J quantity R./pM. has units of moles/kg-s and is the quantity stored in

J J the array R(J,l). An approximation with three grid points is used for

this integral. This is given by Eq. (A6l), Appendix J, p. A27.

The flame velocity is gotten by dividing the value of the above

integral by the density at the cold boundary grid point, RH0(NP2).

IFOW is then tested; if equal to zero, ENTRN continues by calculating RME and RMI from Eqs. (A58), (A59), and (A60), Appendix J, p. A26.

The grid width PEI is then incremented by the amount (RMI - RME)^;-DX.

The condition IFOW = 1 is used at an output point. Then ENTRN

is called for J - I through JHM1 to calciolate the flame velocity from

the production rater, of each species. In principle, these velocities

should be the same in the steady-state flame. When IFOW = 1, RMI and

RME are not calculated.
52

I. Description of the Subroutine SCHM1/1

This short routine calculates the chemical source term arising

from the diffusion of mean molecular weight. (See Eq. (A37), Appendix

I, p. k2U') Concentration values at the control volume boundaries are

required; these are calculated from the average of two adjacent grid

points and stored in BF(J,l). The other quantities required are all

evaluated elsewhere. Values of this source term are placed in the

array SCI ( J, I).
53

J. Description of Subroutine SGHM2(J)/7

This subroutine calculates the part of the chemical source term

arising from the chemical reactions.

SCHM2 begins by testing Gl(K,j); if .FALSE., then species J does

not appear in reaction K and the program moves to the next K value.

If .TRUE., then K is tested to see if it equals 5 or 10. These are the

reactions H + OH + M H^O + M and H f 0^ + M HO^ + M in the H^-O^-N^

flame which have the special third-body concentrations FS(JBODIA, I) and

FS( JBODYB, I) . This routine can also be used for the ozone decomposition

and the H^-Br^ flames because both require fewer than 5 reactions. To

Use SCH^ for other flames, it will be necessary to remove the following lines:

SC200260

SC200270

SC200/;20

SC200/^80

If K is not equal to $ or 10, then B0D1(K) is tested; if .FALSE.,

K does -not contain a three-body reaction and SCHM2 proceeds to calculate

the forward and reverse rates of the reaction. These are put into FOR(l)

and REV(l). At this point, these rates are in Units of moles-m /kg -s

and will later be multiplied by the density to give them in units of

moles/kg-s. If BODl(K) = .TRUE., then K contains a three-body reaction

and B0D2(K) is tested to find out whether it is the forward or the reverse reaction. If B0D2(K) = .FALSE., then it in the forward reaction

and FOR(I) and REV(I) are calculated accordingly. If B0D2(K) = .TRUE.,

then the opposite situation holds. FOR(l) and REV(l) are now combined

in various ways depending on the values of the flags G1 , G2, and G3.

The results are placed in SU(J,l) and SD(J,l). This process is repeated

for all the K values, and the results a dded to SU(J,l) and SD(J,l).

After multiplying SU(j,l) and SD(J,l) by the density RHO(l), the

production rate of species J, R(J,l), is obtained by adding them together. The final desired quantities are gotten by dividing the value in
54

SD(J,l) by the concentration of J, F(J,l), and multiplying both SD(J,l)

and SU(J,I) by ^Tja

For an example of how this routine functions consider the reaction

mechanism discussed in Appendix I, p. A21 . There are two reactions

A + B C + D

A + A+ MiiE^M

v

2 ,

We want to .calculate the source term parameters Sp(A) and Sj,(A) for the

species A defined in Eq. (A53), p. A23.

For K = 1, we have BODl(l) = .FALSE.

B0D2(1) = .FALSE.

G1(1,LA) = .TRUE.

G2(1,LA) = .FALSE.

G.3(1,LA) = .FALSE.

For K - 2, we have B0Dl(2) = .TRUE.

B0D2(2) = .FALSE.

G1(2,IA) = .TRUE.

G2(2,LA) = .FALSE.

G3(2,LA) = .TRUE.

When K = 1, SCHM2 places in SU(M,l) and SD(LA,l), respectively, the

quantities kJjcpQcpjj and - k^cp^"^* Moving to K ~ 2, it adds to SU(LA,l)

and SD(LA, I), the quantities 2 ( k^cp^cp^cpj^p +

^2'^e'%^

^^"^ "

Note that the concentration of the buffer species, is taken to be

the reciprocal of the mean molecular weight, FS(J]yiR, l). SU(LA,l) and

SD(LA,I) are then both multiplied by the density. At this point they

have the values

SU(LA,I) kjjpcpccpp + Sk^p^cp^cpj^ + 2k^pcpgcf^

SD(LA,I) - k^^pcp^cpg - ^'^k^p^cp^cp^cp;^

The production rate R(LA,l) of species A is the sum of SU(LA,l) and

SD(LA,I).
55

SU(LA,I) and SD(LA,I) are not yet equal to the parameters Sp(A)

and Sj,(A). These are obtained by dividing SD(LA,l) by cp^ and multiplying both SD(LA,I) and SU(LA,l) by =

^'^^^i+i " ^i-i ) = PEI0M2(l).

The result is

SU(U,I) = Sp(A) = |-T]n(k^P9^cpQ + Sk^p^cp^cp^ + 2k^pcpgCp^)

SD(LA,I) = Sp(A) = i-Tlfi( -k^ pcpg - ^-k^p^cp^cp^^)

This is the same as that shown in Eqs. (A5l) and (A52), Appendix I,

p. A23.
56

SUBROUTINE SCHM2(j)/7

DO: 1=2, NP2

FOR( I ) =RATE( K, I

) ^^F( LA( K) , I

) -"-F( LB ( K) , I

) ^'^FS( JBODYA, I

) ^^RHO ( I

)

REV( I

) =RATER( K, I

) -:^F( LRA( K) , I

) ^^FS ( JBODYA, I

)

DO: 1=2, NP2

F0R(I)=RATE(K,I)^^F(LA(K),I)^^F(LB(K),I)^FS(JB0DYB,I)TM0(I)

REV( I ) =RATER( K, I ) F( LRA( K) , I ) ^FS ( JBODYB , I

)

DO: 1=2, NP2

FOR(l)=RATE(K,l)'^^F(LA(K),l)^fF(LB(K),l)

REV( I ) =RATER( K, I

) ^^F( LRA( K) , I

) ^^r( LRB ( K) , I

)

DO: 1=2, NP2

F0R(I)=RATE(K,I)^^F(LA(K),I)^^F(LB(K),I)^FS(JMR,I)^TM(I)

REV( I ) =RATER( K, I

) -;^F( LRA( K) , I ) ^^FS ( JMl, I)

DO: I-2,NP2

FOR( I

) ^RATE( K, I) ^^F( LA( K) , I) ^^FS( JMR, I

)

REV( I ) =RATER( K, I

) ^^^F( IJU( K) , I

) ^^F( LRB( K) , I

) -'^FS ( Jlffi, I ) ^RHO ( I

)

6 &

I
57

G2(K,J)

SCHM2/7 (cont.)

DO: 1^2, NP2

SU(J,I)=SU(J,I)+F0R(I)

SD(J,I)=SD(J,I)-REV(I)

DO: 1=2, NP2

SU(J,I)=SU(J-,I)+2.0^KF0R(I)+REV(I))

SD(J,l)=SD(J,l)-^.0^:-REV(l)

DO: 1-^2, NP2

SU(J,I)^SU(J,I) f2.0^-(F0R(l)+REV(l))

SD(J,l)=SD(J,l)-^.0-"-F0R(l)

DO: 1=2, NP2

SU(J,I)-SU(J,I)+REV(I)

SD(J,I)=SD(J,I)-F0R(I)

i

DO: 1=2, NP2

SJ(J,l)=SU(j,l)ifRHO(l)

SD(J,I)-SD(J,I)^^RH0(I)

R(J,I)=SU(J,I)+SD(J,I)

SD(J,I)=SD(J,I)/F(J,I)

DO: 1=3, NP1

SU(J,I)=SU(J,I)^^PEI0M2(I)

SD(J,I) = SD(J,I)^^PEI0M2(I).
58

K. Description of Subroutine CALC (J)

This is the subroutine which solves the finite-differoiica oqiwtlons

(A42), Appendix H, p. A19. To see how this program functions, consider

a U point grid system (N = 3), with the grid index running from 2 to 5.

There will be two unknowns, and 4>^, for which Eq. (A42) yields

two equations,

A' B> C>

^3 \ °3

^

°3

a; b; c;

4 /+ 4

A', B' , C , and D are given by Eqs. (MO), p. A1 3 and the values

and <i>^ at the hot and cold boundaries, respectively, are known from the

boundary conditions. By simple elimination, we obtain

A' A' _ , A' C' , B' , C'

3 4 <J)^ + 3 4 + _2^2 ~

^3 A'B'

1 ^

Let Us go through CALC step-by-step and see how it solves these equations.

There are $ arrays to be considered; A(J,I), B(J,I), C(J,I), D(J,I),

F(J,I), and also a variable T. Initially, A(J,l), B(J,l), C(J,l), and

D(J,I) contain the values of A', B' , C , and D which were evaluated by

SPALD with Eqs. {kAO). F(J,l) contains the values of <i>^ from the previous step. The contents of A(J,l), B(J,l), T, and F(J,l) as CALC is

executed change as follows: (Note that NPl = 4, and the number in parentheses is the number of the step executed.)
59

A(J,3)

A(J,^) (o)„,

4

(4)

4

^

-

B(J,3)

h

B(J,4.)

(0)

4

(5) /=3*2 °3\

\ - ^3

°3

T (o)o.o (3), .V?

.4 D3

F(J,3) (O). / . N ^^<I>^(o.ld) (7) (*3*5 2if^^*2 °3y C' 1 4

C.

\ ^3 / P D

°3 "4-!I^i

'

(0)(I> (old)

4

°4 - ^4*3

°3

A little algebra will show that the contents of F(J,3) and F(J,4) after

these 7 steps is the same as that given by Eqs. (15) and (16).

After solving (A42) for all the grid points, CALC tests for negative species concentrations. If the test is positive, it puts a very

small positive number in F(J,l).
60

L. Description of Subroutine SMTP/1

This routine calculates the total enthalpy source term which is

given by the sura of Eqs. (A55) and (A56), Appendix I, p. A24.. To see

how it works consider the enthalpy equation (A23), Appendix F, p. A11 .

The complete source term is

SENTP calculates

S = T]

The transport parameters were calculated in TRANS and represent values

at the control volume boundaries. They are

HDRM(J,I) = j^l^i-Pl

/ <M> K

HDRM(J,I-1) =M'\-ip1

TRCP(J,I)

C Y- P

TRCP(J,l-l) ^S'^j^Pl

The derivatives are approximated by

p(yj(M)

)| ^ (F(j,l+i);fF3(JM,I+l) - F(J,l)-"-FS(JM,l))-"-ROMD(l)

etc.

,
61

where ROMD(l) = l/(aa.^^ - oo^)

It would be desirable to rewrite SENTP so that the part arising

from the diffusion of the mean molecular weight was calculated separ

ately. Thus, it could easily be neglected if desired.
62

M. Description of Subroutine ZCALC/1

This short routine calculates the value of the laboratory spa tin

1

variable (=Y(l)) which corresponds to a particular value of the

spatial variable cju^ in the Spalding coordinate system. It does this by

evaluating the integral

by means of the approximation given by Eq. (A70), Appendix L, p. A34..
63

N. Description of Subroutine Q[JrPUT( L,HEADNG,RUNID,X, INPUT, m/lX)

This subroutine generates and controls the printing of a variety

of output information. At intervals of LPRINT steps during the integration a single page of output is generated. This contains the values

at each grid point of the species concentrations, the enthalpy, and the

temperature. It also prints the velocities calculated from the production rates of each species, the grid width PEI, and the entrainment

rates RI^ and RMI. Wlien L = UikX, the integration is stopped and 11

pages of additional information are printed:

Page 1 ) Values at each grid point of uj, the laboratory spatial variable

y in meters, the density in kg/m , the mean molecular weight (M> in

kg/mole, and the total heat release rate in J/kg-s.

2) The production rates R(J,l) in moles/kg-s for each species at

each grid point.

3) The transport parameters at each grid point; i.e., the array

PREF(J,I), A.p^ and \p/C in kg^m^'-s.

J P

4., 5) Values of each rate constant, both forward and reverse rates,

at each grid point.

6,7) The rates, forward and reverse, of each reaction at each grid

point in moles/kg-s.

8) The maximum rates, forward and reverse, of each reaction and

the grid point index at which they occur.

9) A list of the maximum rates of each reaction ranked in decending values. For each maximum rate, its K value and direction are speci'

fled.

10,11) Heat release rates, both forward and reverse, of each reaction

at each grid point in J/kg-s.

OtfrPUT begins by printing the step number L, the step size DX, and

the integration time X. It then prints the heading for the profile

colums and follows with F(J,l), F(JH,l), FS(JI,l), VEL(J), PEI, RMI,

and RME.
64

It tests for L if .TRUE., it proceeds by printing the

values of the grid parameters, N, Nil, R"iE, and OMR; the index parameter

JH, JRAD, JM, jm, JT, JIR, JBODYA, JBODYB; the temperature on the

cold side of the flame TGOLD, the pressure PRESS, and the value of INPUl

FLAG(2) is tested; if .FALSE., then the buffer concentration is

being kept constant, Tiiis means that its production rate R(JIM1,I) is

zero. In this case, the production rate for the next mcajor species

R(JI[M2,I) was not calculated by SCilM2, OUrPUT calculates it from the

relation^jTM^R(J,l)WT( J) ^ 0,0. This section is skipped if FLAG(2) ^

.TRUE.,

OUrPUT next calculates the quantities PliT(j) = h.M. = h'^", and

J

'r(J,I)^^PUT'( J). This latter quantity is the total hsat

release rate in l/kg-s. It then prints OM(l), Y(l), RHO(l), F3(JM,l),

and FOR(I), followed by the arrays R(J,I) and PREF(J,l).

It then tests FLAS(9); if .TRUE., it prints the rate constants.

This requires two pages. Each column represents a grid point whose

index value appears as the heading. There are 1 6. columns per page

giving a total of 32 grid positions. This may not be sufficient to

cover all the grid points. Since points near the boandaries are not

very interesting, we start printing at I = ISTART, where ISTART =

in the present version of OJTPUT. Its value can, of course, be changed

to suit the User's requirements. The next page starts printing at

ISTART ISTART + 16 which equals 20 in the present case.

FIJlG(lO) is then tested; if .TRUE., the reactions rates are calculated at each point in the flamo. Note that FOR(l) which contained

the total heat release rate is now used for the forward reaction rate.

The rates are calculated from scratch; i.e., from the Amlues of the

rate constants MTb;(K,l), RA.TER(K,l) and the concentrations F(J,l).

Also calculated at this time are the heat release rates for the foi^rard

and reverse reactions. These are stored in the arrays HV(K, I), and

HW(K,I), respectively. Note that we have treated the cases K - $ and 10

in the same manner as discussed in SGHM2. To use OUTPUT for flames

(other then H^-O^-Np) which have more than A reactions, it is necessary
65

to remove the following lines:

01/101680

ouroi69o

0UT01960

O'JI 02100

The reaction rates are then printed on two pages in the same format

as the rate constants were printed. Note that ISTART also appears in

this portion of the routine.

The maximum rate of each reaction is now determined. This is

accomplished by running through the grid index I for each reaction K

and setting the vBriable XM/IX(K) equal to V(K,l) if V(K,l) is greater

than the previous value of HylAX( K) . V(K,l) is the rate of the forward

reaction and is in single precision. The same thing is done for the

reverse reaction and yields YRA.X(K). XM./IX(K) and YMX(K) are then

printed along with the grid index value to which they correspond.

A new single precision variable Z(K) is then defined; Z(K) =

XMAX(K) and Z(kTMaTE) = YKkl{K) for K = 1, KRATE. Thus Z(K) runs from

K = 1 to 2^^ KRATE. The ranking of the maximum values of the rates is

accomplished by running through Z(K) and setting ZTAB = Z(K) if Z(K)

is greater then the previous value of ZTAB. KTAB is set equal to the

value of K for this maximum - taaximtom rate. Also defined is a logical

variable SKIP(K), K = 1, 2'^KRATE, set initially to the .FALSE, position.

At this point OUrPUT sets SKIP(KTAB) = .TRUE.. Then the Z(K) are reexmined for the maximum value, but now Z(KrAB) is omitted. This is

accomplished by testing SKIP(K); if .TRUE., then KCAB is skipped and the

next K value is compared with ZTAB. The list from which the maximum

value is selected thus becomes smaller as each maximum is removed. In

this way the ranking is accomplished. The maximum values of each reaction

are then printed in order oi decreasing magnitude.

Finally, OUTPLTl' prints two pages of heat release rates for each

reaction. Note that ISTAxRT appears here also.
66

SUBROUTINE

OUTPUT{L,HEADNG.RUNID.X.INPUT,LMAX)

WRITE: L,DX,X

HEADNG(M)

F(J,I),F(JH,I)

FS(JT,I)

VEL(J)

PEI,RMI,RME

WRITE: N, Nil, NEE, OMR

JH, JRAD, JT4, JMR, JT , JTR, JBODYA, JBODYB

T COLD, PRESS, INPl/r

J1M3 =JH-3

i

. DO: 1 ^2,NP2

I

in

DO: J :1,J1IM3

SlJM^SLJMt-R(j,l)^^WT( J)
67

©

I

OUTPUT (cont.)

DO: I:=2,NP2

F0R(I)=0.0

T=FS(JT,I)

DO: J^1,JHM1

P Ur ( J ) =WT ( J

) 'K h°+d . T +|-e ^. T +|f^. T^

)

FOR(l)=FOR(l)+R(J,l)^^PUr(j)

WRITE:

I,0M(I),I(I),RH0(I)

FS(JM,I),F0R(I)

VfRITE:

RUMID(M)

HEADNG(M)

R(J,I)

I

WRITE:

HEADNG(M)

PREF(J,I)
68

ICOUNT-O

ISTART=:/^

ISTOP^ISTART+1 5

I

WRITE:

I,RATE(K,I),RATER(K,I)

ICOUNT-ICO

OUTPUT (cont.)

I COUNT-;
69

OUTPUT (cont.)

ICOUNT =0

I START

E

ISTOP^ISTART+15

DO: K=1,KRATE

DO: I=ISTART,ISTOP

F0R(I)=RATE(K,I)^^F(LA(K),I)^^F(LB(K),I)^^FS(JB0DYA,I)*RH2(I)

REV( I ) -RATER( K. I

) ^^F( LRA ( K) , I

) ^^FS ( JBODYA, I ) --RHO ( I

)

SFAG=PUr( LRA( K) ) -PUr( LA( K) ) -PUT( IB( K)

)

HV(K,I)=SFAG->^F0R(I)

HW(K,I)=-SFAC^REV(I)

l-'O lU I ) IvATF, ( K , I ) F ( I A ( K) , I ) F ( M( K) , I ) FS ( JBODZB , I ) "-RHa ( I

)

Rp;V( i ) KATKlU K, J

) i!F( ]RA( K) ,l)^'^FS( J'BODYB, l)^^RHO(l)

SFA C- r U{( IiiA( K) ) -P Ul' ( L/\ ( K) ) -P m( IB ( K)

)

HV(K,l)-SFAi;s>'FOR(l)

PIW(K,iy^-SFAG^4lEV(l)
OUTPUT (cont.

FOR( I ) =RATE ( K, I ) *F( LA( K) , I ) *F( LB ( K) , I ) --RHO ( I

)

REV( I ) =RATER( K, I ) *F( LRA( K) , I )*F( LRB ( K) , I ) --RHO ( I

)

SFAG-PUr( LRA( K) ) +PUr( LRB( K) ) -Pl7r( LA( K) ) -PUI( LB( K)

)

HV(K,I)=SFAC-^F0R(I)

HVJ(K,I)=-SFAC*REV(I)

FOR ( I

) ^PJVTE ( K, I ) ->^F( LA ( K) , I

) -"-F( L3 ( K) , I

) -FS ( Jlffi, I ) -RH2 ( I

)

PJ;V( I ) CRATER ( K, I ) *F ( LRA ( K) , I ) -^FS ( I ) --RHO ( I

)

SFAC=P ITT ( LRA ( K) ) -P Ur ( LA ( K) ) =P Ur ( LB ( K)

)

HV(K,I)-SFAC^-F0R(I)

HW(K,I)=-SFAC*REV(I)

FOR( I ) -RATE ( I ) ^F( LA ( K) , I ) -^FS ( J>1R, 1 ) -RHO ( I

)

REV( I ) =RATER ( K, I

) -"-F( LRA( K) , I

) *F( LRJ^ ( K) , I

) ^FS ( JMR, I

) *RII2 (I

)

SFAC=P Ur ( LRA( K) ) +Pin: ( LRB ( K) ) :^P l/X( LiU K)

)

HV(K,I)=SFAC^-F0R(I)

HVJ( K, I ) --SFAC-"-REV( I

)

/(K,I)^F0R(I)

—

W(K,I)^REV(I) -f

—-

WRITE:

K,F0R(I) ,rSV(l)

I

ICOUNT-ICOUNT+1

®
OUTPUT (cont.)

DO: K=1,KRATE

I

XMAX(K)=0.0

YMAX(K)=0.0

WRITE:

IX ( K) , XMAX( K) , lY ( K) , YMAX ( K)

i

DO: K^l,mTE

Z(K)--^XMAX(K)

Z(K+KRATE)-^YMAX(K)

I

XMAX(K)=V(K,I)

IX(K)-I

YMAX(K)=W(K,I)

IY(K)=I

KK-2^^KRATE

ICOUNT=0
72

OUTPUT (cont.)

TAB=3HF0R

NTAB-KTAB

I

ICOUNT^ICOUNT+1

I

WHITE:

1 COUNT, NT A1>TAB,ZT.-lB

I

SKIP(1CTA13)-.TRUE.

ZTAB-Z(K)

KTAB=K

TAB=3HREV

NTAB=KTAB-KRATE
OUTPUT (cont.)

DO: K^1,KU/VTE

1:^2, NP2

V(K,I)^HV(K,I

W(K,I)-HW(K,I

ICOUNT^O

ISTART=/,

I STOP- 1:

WRITE:

K,IIV(K,I),HW(K,I)

IGOlMT-ICOUNT-f-1
74

IV. TEoT CASE

As an example of the use of this program we have utilized the

H^-Br^ flame. The input data were taken from a paper "by Lovachev

~ 2 3

and Kaganova. Spalding and Stephenson have also used the Spalding

method for this flame. The following input information was used:

Thermal conductivity:

where is the thermal conductivity at the temperature T^ of the unburnt gas. It is assumed to be independent of concentration. For

T^ - 323^K, is taken to have the value 3. 34-9x1 0"^J/m-'^K-s.

Diffusion Coefficients:

A. = A. (T/T y-^^/F

J "ju ' u

'

where P is the pressure in atmospheres and the A • ^re considered to

be independent of concentration and are the values at temperature T^.

For T^ - 323°K, we use

\a = " ^HBru = O-IOlxlQ-V/s

^Hu •^^5x10"V/s

^Bru ~ 0.15^-.x10~''m'Vs

Thermodynamic properties:

The specific heat of the mixture is assumed to be independent of temperature and concentration. The value used was C - 532J/kg-°K = d..

P J In this case, the coefficients e. and f. - 0 and all species have the

same value of d .

.

J The enthalpies of formation as taken from the JANAF tables are,

h3^(2Q8°) = 1./.009xl0"j/kg h|^ (298°) - 1 .93505x1 O^J/kg

h^(208°) = 2.16^1x10^ • h|L (293°) = -^^.^SOVxIO^

h!! (298°) ^ 0.0

^2
75

Note that the value for liBr is nli/'ihtly different from the value of

-4. 5304-xl O^j/kg Used by Lovachev and Kaganova.

For the program we want values for h° - HREF(J) = h. - d.Tj,, where

is the reference temperature for h^. . We have

J

.6. hg^ = 1.2425^0x10 J/kg

h° ^ 2.162510x10^

h. -1 .583691x105

h° = 3.513587x10^

^^2

^HBr

~ -6.064-374x10^

Rate Constants;

The reaction scheme used is

Br^ + M 4* Br + Br + M

"2 ^

H + B + M

Br + ii air + H

H + Br, <i HBr + Br

(1)

(2)

(3)

U)

46.09 kcal/mole

104.20

16.70

-41. /^l

The forward rate constants used were

k^ = 7.03x10^''T'^exp(-2.3655x1oVT) m^mole-s

7.63x1o'^T"^gxp(-5.4356x10'Vt)

k^ = 3./i6x10^Texp(-8.3547x10 Vt)

t' - o k, -~ 6.42x10 T2exp(-',. 5866x1 0'V'^')

Pi-ee energies of reaction in kcal/mole were calculated from the formulas,

P - 46.1756 - 2.49467x10"%' - 6.5671x10 '^T

-7.n2

A.F° 104.426 - 2.3895/ixlO"^' - 1 .4.176Ux10'M

AF° = 16.532 - 1.4961 5x1 0"^T - 1 .91 1 12x10"'^T^

A-F^ = -41.7184 - 2.54-835x10~^T ^ 5.69792x10"V
76

For the reverse reactions Lovachev and Kaganova use

kjj - 3.63x10^ mVmole^-3

= 3.63x103

- 9. 06x1 0 ^T2exp( -8. 606/^x1 oVy) m^/mole-s

k' = 6.$2x10'^Texp(-2.U^OxloVT)

4

These formulas yield values

of the reverse rate constants which are slightly

different from those calculated with the AF° values via the equilibrium

constants. For comparison, Table I gives the forward rates, equilibrium

constants K , reverse rates

P

calculated from K , and

P

the reverse rates of Lovachev and Kaganova.

The other conditions

for this test case were,

pressure = 1 atm, cold

boundary temperature =

323°K, cold gas composition, Br^ = /+Omole^,

balance, H^.

The format for this

input data is shown in

Table II. Cards containing the symbol @

are control statements

for the ms UNI VAC 1108.

3

-1000

y (mm)

Figure i

Concentration Profiles for the Major Specie;
TABIi'^ I. Rate constants for the E^-Bv^ flame as a function of

temperature.

V

1 2918x10"''^

K

P

1 .8963x10

k'

1

1 . 2621 x1

0^

4.00 7.

-20

600 2. 1 668x1

0"'^

5.2174x10'

-12

2.04^7x1 0^

800 3. 5799x1

0"^

8.9451x10"

-8

2.6272x10^

1000 1

.

1850 3.1882x10'

-5

3.0499x10^

1200 5. 5762x10^ 1 .6384x10"

-3

3.3513x10^

1400 8. 6271 xlO'^ 2.7841x10'

-2

3.5598x1 0^

1600 6. 6700x1 0^ 2.3689x10'

-1 3.6967x1 0^

^1

3.6300x103

Tempera tiore •^2 K ^2

400 1 .8373x10

-48 P

1 .9216x10

-52

3. 1383x10^

600 5.7569x10'

-29

2.3267x10"

-33

1

.

2182x10-^

800 2.9600x10'

-19

8.6949x10'

-24

2. 2348x1 0^

1000 1 .8881x10"

-13

5.1000x10'

-18

3. 0379x1 0^

1200 1 .3529x10'

-9

3.7477x10'

-14

3. 5547x1 0^

1400 7.4921x10'

-7

2.2544x10'

-11

3. 8178x10^

1600 8.4011x10"

-5

2.8378x10"

-9

3. 8868x10^

^2 , 3.6300x10-

Tempera tu.re

400

600

800

1000

1200

1400

1600

3

1 .1753x10

1 .8616x10

8.0661x10

-2

1

2

8 . 1 41 0x1 0^

3.9318k1 0''

1 .2402x10'-"

2.988^)x10-'

K

P

2.0411x10

-9

-6

2.1356x10

6.9745x10"^

5.6918x10

2.3220x10

6.3739x10

1 .3658x10"

-4

-3

-3

5.7582x10

3.7170x10

6

D

.7

1 .4303x107

.7

k'" ^3

2.1073x10

5.2875x10

D

6

1.1565x10' 8.7390x10

1 . 21 1 6x1

0

1.6933x10' 1.5319x10

1. 9457x1 0'^ 1. 8332x1

0'

2.1881x10"^ 2.1163x10'

Temperature

400

600

800

4

3.1769x10

6.1978x1

0

9.0324x1

0'

K

P

2.0142x10

4.7888x10

7.1752x10

23

15

11

^4 ^4

1.5773x10"''^ 1.3729x10"''^

.-8 -8

1 . 2942x1 0

^ 1.1 838x1

0

1.2588x10"^ 1.1968x10"^
Temperature

1000

1200

UOO

1600

4

1.161 2x1

0

1 .3962x1

0'

1 .6117x10

1 .8112x10

K

P

3.5582x10

1 .0152x10

7.8715x10

1 . V^02x10

9

6

3.2634-xlO

1 .3752

3.1826x10

1 .3611

-2

2.0^75x10' 2.0386x101

1 .5385x1021 . 5802x1

0'

78

TABLE V. Continued.

Reverse' rate constantc calculated from the expressions given by

Lovachev and Kaganova.
79

O

CV

cv \—

T

\

o ^

S

LO Q o

M

pj EH Q Q

W 3 O O

W CQ M ^-D f-D

O

f3

99

o

W II II g Q i-D I-;) M cni M f-;) hs

»>£e! I—I W S 1-3 l>< X!

fin Q

Ph O

is; O

1>

Qo

o

m

I

^^

£3 CD @ ®

CO

rs

CO

. EH

C/) £h Hft ft (i) FH IX, m

a H

i-i li, H EH M fe

I—I C/5

9

COw

f-l

J--^; E^; M M M pq @

Pi

CO W

Eh

in

O O -t O

00 o

o . .

CM O Q Q O

O OO

CV [> M V- O S O X- cr> r- Q P c.-i fr; --d OA O CV O O l>qpQMc^il^^0Cx^f\'00 Q CQ V- V- o vD

T- o in o ^.0 • •

OA

I I I ! I Q Q Q Q Q

vD 00 vO OA rv v- O v- 00 O o o o

o

o o

CV UA 00

????? Q n n Q Q

o o o o o

o o o o o

o o o o o

o o o o o

o o o o o

o o o o o

o o d d d

??Q Q

O O O

O

o o

o

9 9? Q Q Q

O O OO

O O

O O O

o o

o o

o

o

o o o o o

OA OA OA

Q Q i

-4 -J- cv cv cv J> r- o

C\i OA OA OA'

I LTA IXA UA Q . . . C\i O O O

O « O Q

O O O V-

• -vt ^ O ^— UA UA vO

8 0i C\i 4vX) 00

• Cvl V- UA

. . . ^ cv Y- 0^ I

cv cv I> i>

OA HA

UA MA O • . Q

O O • o

ifA o

+ + « Q

I> -<f I 00 Q l/A OA • --t UA O V- v,0 UA o T o • • Q O

r<Avp . T i.(A OA

EH EH

Eh eh

UA UA

OA

• '-n o o

Q Q

o o

o o

UA •

EH EH

Eh Fh OA rv Q Q

t> vO

-4

LrA 00

UA

• * 00 'n

ss

o o

o o• UA

^- --J- T- T- O (5 Q Q vO (AJ

• • • • O £> t> OiA v£)
80

« . Q W O EH CP

EH

o —

i4

o

EH w

<5 «aj

II II T- ^

•N ^»\ PJh ^ ^'S p£-| f2£_(

_ p:^ W -—-—^-

E-iPhWI-IEhOC!--^

W Ph S W H

S Q

M M l-H CO

o

13

•H

-Pc

o

o

EH

EH EH &H EH Q Q O VUA O

Q Q

V— O

I cv

I I I I Q n Q ti

L, \ t> CPv V- O O 00 o--4^ in

cv cv

I I

CV CNi V- I

cv in iXN O ~4

-4 OA

(^4 is; O M

il5
81

Since the format of the output makes it unsuitable for reproduction in this report, we shall

illustrate some of it in graphic form. The steady-state

velocity of this flame is 0.24-1

m/s and the hot side flame temperature is 1$18°K for a complete

reaction. This velocity is ^.2%

lower than that obtained by

Spalding and Stephenson; the

difference probably arises

from the slightly different

values uned here for the reverse rate constants. (See Table

I.)

Concentration profiles for

the major species are shown in

Fig. 4-. These are plotted as a function of the laboratory

spatial variable y in mm from the

hot boundary. Also shown here

is the temperature profile.

Profiles for the two trace species, H and Br, are shown in Fig. $

along with the enthalpy profile. This latter profile has units of

J/kg. The rates of the inajot- ro.-jctlonn are ,",hown .in Fig. 6. The +• sign,

in the figure denotes a forwjrd rate. Figure 7 gives the heat release

rates for the major reactions and also the total heat release rate as

a function of distance through the flame.

0.4

y (mm)

1 .8

Figure 5

Concentration Profiles for the Minor Species
82
83

Figure 7

lient Release liates
84

o

o

ooooo

o

ooo

ooo

oo

o

o

o

oo

o

o

o

o

o

oo

o

o

o

o

o

o

n <t 10 « h- CO a>

o ro in 00o CVI ro * 10 GO 0>

o ro <t

oo

oo

o

o

o

o

o i-H r-l >-i CVJ CVJ {\J (VI (VJ (M CM (M ro fO fO ro ro

o

oo

o

o

o

ooo

oo

o

oo

o

o

o

o

o

ooo

o

o

o

o

o

o

oo

o

o

oo

oo

o

o

o

o

o

o

o

o

o

o

o

o

o

oo

o

o

o

o

o

o

oo

o

o

o

oo

ooo

o

oo

AQ

A

Qo

Qo

o

o fto

o

QQ

Q

a

o

Q

o

Q

a

Q

a

Oi Oh Oi IX a. cu 0. 0. IX 0. CU 0* (X ft. CU a« Oh (X CU Oi 0. 0* Oi a< OU CU (0 M 01 0] CO CO 01 01 0) 01 01 0} 01 01 01 01 Ul 01 01 0} 02 01 01 01 01 01 01 01 01 01 01 01

CMn

ro

ro

r")

(Jl

M < I W I CM <

PU

06«

lJu bu

«<bw

ro • O O• •< CM h) V b• CQ ro

•-• • 0* « H a; « u • O < C\i W CU « z

cu « •

AIww%

(Li *wum<

pqM

0]• H

hIH

(9

COw z<

-<

ft. < ^ w M CK

CM in p ^ n w• M ^ u a w

CQ 0} u u M M

QO

PQ* Q

06H

ft. os

«<06

CM «

:«• • w ^ w a z m •

row

CMXoMu

ft.

ro

in

in »i

in ^ W ftl >. O • ^ K

^01 • H

01 • W lO

OS >H ro

ft.

w ^ • ft. • o u

iJ bu o

a

H in ft: • - «

ft* • O

Ui ro • <

OS in 10

ft. «-t »^ z\01<>M\z

«paa • H n z

01 < J O < Ui

Q ^ a

ro ID w U Q ft. a «

f5< ro U w a 0. o a• ^) M • a

10 ro

ro ro

<t

H

w

a

a CQ • ^ ro ro <t

4-

w

w

a

ao

ID ft) \ \ CM ro

-< < • ft. ft. u w

in

QK

t-i • CU w ^ • w in « ro

ft.w

in ft] ft)

fts

a

m 06 •HW WCM \ 01 01 H Z Z M < < \z Nz \z Nz \z Wz z O O C C ID V C a a a a a a a a a a a a a a C 'D C C V V C w u o u o u o

Q

ro •

CM

ft) •

06 ^ • »^

rowD

ftS ^»

ro

pq

n•

ro

ro z ^ <

in

S

w z a o

CM

a

06

a

pc• u « p\wz

•<zaa

IDU

(O>•

ro

<t ro •

in

wp

>^ >^ 01 « • ro

row• •

in

M

ft]

w ft.

ft,

\

\

u

QC -< < y \ \ z z a a a a u u

o

rowHO

ft) O

ro

•<H

ft]

ft)

HO

•< • ft) ^ ftS o •

row

o

O ro

ft. CM a © « ft) < •

ft. H

o

06 ro \ ^ z ^ o p a © a ft)

o

roOo

ro

ft) O

ft)«o

row

-<a

ft)•

ft] H

06M•

ro•o

ro

06

U)H

•0CK

ro

in

ro ^ ^ ro •o •

ro in

ft]

w

H< w

0< 01 \ \ U »H H a < u « 0] \ \ z z a a a a u u

p

0)

ro• ro

in

w

wZ 9 ft) 01 01 \\CM H a z y ft) 0WJ 01 z z V <o a a a a

JO 10 o u

ro

inp

ro• W

10 ft. • «P

ft] a

06m

ro

in

ft)

rf)

in

a

06 ro

in in <\

ft.

ft] \zaa

10u

•4 ftfi > 10 \ 06 06 ft. H 03 U ^ \ \ z z c o a a a a © T> o u

CMz

CM»

•J

CMH

•J*

CMTO*

CM

h)

ft)

<t X

h)

CMo

CMn

h)•

CM10

paH

hJ z n w « ft.

h]Z

•>

ft)

V ^ u P VI ZM <

S3\ \ z z

ID

©

a

a

aa

©©

o

roo<a

ft.

ft$

o

ro«

ft! h)O

ro<

06 • in o

row

ft.

ft.

o

row

fOo• in o

rowo

-<»J

ft.

ft]

o

ro

CMoo

ro

o

row • in

ft)1 o

CVJ

O

w

ro

W

w <

h) ft. >^

ft] l4

06 ft. \ \ z z

C> ID a a a a

c>

©

in

•<

h-)

o

ft.

\ ft) \ u

fti

H

W

M < * J 01

ft.

c

\

\ z z

JO

V

aa

u (Jw

O •-' P

0(

M

J z z M M

0) OT Z Z

ft] ftJ a ^ M M P P
85

ooooooooo

in«ONooo»o-»<viro

ooooooooo

oooooooo

<tin\ONooo»o«H

oooooooo oooooooo

(/laitnaiuimaiai

o o o o o oo o o

CM ro in vo o> o

in in in in tn mm m vo o o o o o oo o o

o o o o o oo o o

0 Q Q Q Q QQ Q Q

O4 (X4 (Xi (UCk (Li Ou

01 (n 01 Oi ui mm m m

ooooooooo

•-«tvjr)4-m\0f^<oo»

ooooooooo ooooooooo

mmmmmmmmm

o

QO

•*m

o

po

«

«*\

m

IQ

m

ft< \

W (VI 0< »

«• O \ • 2a o

o o

• o

* • ro * ^ m \ *\ m

•» u

O H

\

O

pO

Q

O

• m

S m

o

2

HiK

H

a

-<

(0

(M

a

(Vl

iZi

CM»

©

CMH

1-1*

«

CMX

M

U

•-4 X

•

CM©

CM»

1-1

•> U

CM Q

TOH

PQ • •J m W

0 Q

• * * W <

* « M

n Q

H 0 CM CM w 0 0 •< 0 M 0 0 CM 0 CM 0 ro 0 in M

(M (M vO CM ro K « W • • W • • w •< m H m H m H in m H

w < < "< w w

Q n n 9« Q n

•< PS 06 < P$ < •< PS

• M u Ui M w ©

M PS « (b K K u. 0 (M 0 0 0 0 m

0 0 0

CM CM CM

04

« • % m

Q

• O CM O M

m ro 000 Aw

» • • m H

< m o w <

n w n Q a ^ < < < 06 H H H W ©

iJ ^ Pi b

o

o

in

g

n

n

06H

2i

m

06M

H

M

06W

H

M

t

06H

•

PS

91

P

<

06*9

«i m

CM

10 M

O

in ^ • m H"^

55 S5 w

» • p as*^

CM ro "< psa

0. (ti ©m

% 2 PS

roo

m

H

»9

»-« * n *^w

:« • Oi M M

P4 H

P P • < M

P« P^ PS P P

< < « * P4

06 06 m

m

K * w •» 06 • m « W p« m •1

n) m p w

(Li P • (b

P{ CM M

M W 0.M u i

n

w w ^ H w

CM ro

o

-<

Pi n *

• a* CM D P

a << (Q 06

CM ro

I I n

as

CXI

ro

CM CM N OOP

m m m

• WW

m m H

w <

P 9!

•< 06 U U) ©

PS PS (b

P

<

ro ^

O H ^ m

m H

>^ <

p a

<< PS

Ui V

06 (b

2S

m

p

<

Ui

m O

o

m

H m

-< w

3 p

PS <

© Ui

PS (b PS

CM000ro »^ m >^ ^
86

o o o o o o o O o o o o o o o o o o o o o o o o o o o o o o o o o o o

o CM ro * tn CO o ro IT) h- 00 o (VI ro 10 N oo o n

N ^ 00 CD CO oo 00 OO 00 oo 00 00 (^ 0> 0» (^ o o o o o

o O O O O o o O O O o o o o o O o o o o o o o o o o o o o o

o o o o o o o O o o o o o o o o o o o o o o o o o o o o o o o o o o o

a o a o Q a o a o Q o o Q Q o Q o o Q A o Q o

a* OL. 0. 0. a« 0. Oi cu a< Oi Oh cu a. 0* 0^ a« cu ou

CO CO CO CO CO CO CO CO CO CO to CO CO CO CO CO CO w CO CO CO CO CO CO CO CO CO CO CO CO CO CO CO CO CO

H

H

01H

H

(0

(\j p

O 2 O

•

o

in

in (V)

CO

H M

(0 w

2 « ^ H • H M

Q X «0 X <

O Q

fi o

a. H

<^ A A (y| 00 ro o • o M ^ ro

in (vi tn (fa

M

w

O

•<

i-l m

M

w

O

-<

»J

[fa CQ

CM

in

H

M

w

PQ

CQ

M

H

S

M

M

w

CQ

CkS H M (VJ (VJ [14 U4 «l U H H H

•< « H «< >< < w

u W Cki 06 • o • M M •< • • « w

M m • h4 w M R R R • < w • • M M

H (Vl CQ W C\J (VJ • •

(4 • O Q Q w

m M [LI ^ w M M M (a • • • (^ a; » (4 w

00 (\j (VJ (VJ w M H • M Q M ro ro m H p u y <

>^ o o o (0 H PS (i4 -< M

O (\J >< w • • • [14 U4 (fa [fa « • O (M H in in in CU CO [fa 06 •<

[fa [4 • < « X

M % w w w w w iJ w w in H w

• • M M

«^

(VJ • w <<^ (VJ ro in in in

• in <^ «< 03 M CM o o o h) o o

I O

s

o o o o o o o in

in » in ro M [14 in in (VJ (VI (VJ ro (VJ (VJ fO

o

(VJ

(QH

OS04W

U4

CO

[fa M P

•<

[fa

(VJ

H

n

-<

HinHwinHinHinl-' <w<R w'<w<w^ aipscvipaipapa K-^Qi{J-<06<06<06 Cn4V<WO(dV[x4V

[fa«(L,(fa(ki[faH[faO([fa

in in w < w ^ •-• M

R

(14 H Q a p

< -< OS <

W (4 ID (14 M K [fa K

O O

in in

oo o cvj O iH r4

in in in in

Oii

(fa

in

•-• in

in w

© o

^ in M M w

w w P

CQ CJ •< O O U4

Ui Ui ID 06

in

<

a!

(fa

vO

in

in in in H

w w w ^

P P P 2

•<><•<«

U4 (U [14 ^

OS « CK [fa

fOo

o

(VJ

in H

p :b

04 V

in in in H

w w ^

p p a

< < 06 a

<

[14 [i4 04 ^

(K [fa PS « OS (fa

ino

o

(VJ M

• *

in CO w

p ij

-<

UJ <

K U

<

(fa• H

V

2i

ro

RV

H

CO

ino

o

(VJ

in
87

o

o

o

o

o

oOOOO

O

O

O

oO

O

oo

o

o

o

o

O

oo

o

o

ooo

oo

o

O

O

in oo a»o (\J ro in »o h- oo

ow to inN 00o (\j fO in (Q

oo

O

o

o (VJ CM (M CM CM CM tM (VJ (\J ro ro ro ro ro ro ro fO ro ro

o

oo

o

oO

O

o

o

o

oO

o

o

o

O

oO

O

o

O

O

o

o

oo

o

o

o

o

O

ooo

Q

QQo

O

QQQQ

Q

Q

Q

QQQQ

a* Q< 04 (U CU (1« (U Oi (b At a« CU Oi a* Pu PLi (b (b (b (b (b (b lb (b (b 0. (b (b

CO CQ (0 (0 (0 CQ CQ (Q CQ CQ CQ CQ CQ CQ••

CQ CQ CQ CQ CQ CQ CQ CQ CQ CQ CQ CQ CQ CQ CQ CQ CQ CQ••CQ

N.S 1•• ••• M • • < • Ui 1 OH • t • • u

(0 • ro • H I (b • < • • M * * * u H 1 1 CQ 1 > 1 » % • H 1 • ' « 1 M M • 1 u • w w • H « H • (b D t • Ok 1 • TO I M • 1 • • 1 M * 1 «^ 1 1 M M 1 (b 1 w 1 H < 1 ro ro • • H • H H 1 OS

CQ 1 1 M 1 to 1 1 CQ 1 \

m • • • 1 CQ 1 1 M <j M M i • w r-t w 1 M

3B • H

H

H 1 M

CQ • CQ

O

n

n * Q

Mo * • 1 CQ H 1 CQ M M 1 • 1 w w 1 (M to CQ 1 H Q I m U • 1 •< n * H • '0 M w • M CQ H P • M • • • H 1 U Oi M 1 w 1 M w M w K 1 • M H • 1 H O I • M W w in M in 1 M H I • M M • M as • 1 CQ H P< (b I H -<* 1 (b

a. Oi 1 ro

n *^

s 1 « 1 • n H in Q ^ n Q 1 o « H !S o o (M CO in N (VJ (VJ w K < If) « o • O n

in * o • o in • O

J CM t M 1 w (Ni vO K X X X

M

w

o ?cQ

w

O • «-« • i-« « « « • «» f* * H

s

n • M

M (Q << ai CQ vO

ww

W

W

w

w • < . « M U4 w H H w H w H w H w w H w H < D w •< w < u < w <! w < H M «< u < D

uu

M

Do

U

H

H

H as

HH ae

H :»

HH

H asH

w

o • H •J

w

MM

w

M

MM

M

M OSM •J

(k • U • < •<

W « K TP « « -< •J CQ u o M » (b (b (b (b (b (b O o o o t-« (VJ CVJ o in o o in O in o N «0 «o «o yo «0
88

oooooooooooooooooooooooooooooooooo o^wfO'*in«oi^ooo»o«-«(\j(0*«n*oNcoo»o^{\jrO'<tin<or'00t^o»*w»o

oooooooooooooooooooooooooooooooooo

•4H

(M

H

»}

M

H

•J

(J

«<

O

t o

M

•D

IXZ

O

<^

ro ro

o«

• *

• «M

* H M

m • H H

*-i w w

(Q b (fa«

til o

•-• K • C» >-» H • *

tD M M

O • w

rH (fa (fa w w • s

w w

»

• Q

ro ro N

(M o o • o •

vO H X

15 • • • r-» Ul N Ni^ vO w

Ui w w w H w H

H w w w >< (U H H H H 3

w M M M M M

(h <t> K w K (K 1? M » (fa (fa

•* ro «

O in O

\0 «0

04H

(D

&

(4

H

H

H O

(0

(4 H o

ro * O M

(Q O o w

W •* H (M M

M ro w w ID H o n n

P6 H •J TO O (M ;« * W (D H tfa ID TD

a* ID Z «

•-t O w s H (fa M M

• o •H • in (4 z (4 w

•< a W o • (fa (fa (fa O (0 (0 (Q (I) y n * • ro « a s

• (M o • H M <^

• < O M t o Q w w CVJ H H H

H •J * O D H o ro u M p (fa w w w

• < t H Z • N z P Q (M

• < H V < M • z V M :i O ;« 3

m h) H n H M H « • H n V V V

w < w < •J Z w M •J M Z M M M M

W (fa (fa C> H < (fa (fa m < (4 Ui Hi (I] Ui H M M M M O (fa Cfa U (fa BU (fa O (M o f-« i-t o o O O ro ro o

u u u
89

o o o o O o o o o o o o o o o o o o o o o o o o o 0 0 0 0 0 0 0 0 0

00 o ft ro •* (0 oo o ft ro in oo 0» 0 tvi CO * 10

r- oo 00 00 00 00 (0 00 oo 00 00 0> o a> 0> 0

w

0

(M

0

(M

0

w

0

(Vi 0

CM 0

CM 0CM O O o O o o o o o o o o o o o o o o o o o o o o o 0 0 0 0 0 0 0 0 0

o Q o Q Q Q o Q o A o o o » o a 0 Q 0 a Q 0

a, 0. PU Or (X< eu ou 04 cu 04 04 04 04 04 0. Oh 04 04 04 0< 0. 04 04 04 04 Oi 04

01 0) (0 (0 01 (0 01 01 01 01 01 01 0) 01 01 0) 01 m 01 01 CO 01 01 01 01 01 (0 01 01 01 01 01 01 01

roH

01 00 w o

H r^

row o H

01 N

U

M H

w

(l4 V

«-» roro roro rofo ^ ^

'^ CUO4 O4O4 &4O4 roro ^

roro S5!5 5555 izSe; OhOi

O4O4 «»• 5555 roro ^ ^ ^ v4*H r44H •<» O4O4

4t« II 11 II r-t ^

w1 ^ MM MM MM 11 •4>

II •« MM Mr4 MM ^ ^ 0m, ^ ^ ^ aas MM

Wtd HH OStfl WW WW

•^i^, ^^ mm MM

*<» ««» M^ MM MM MM K« HH II II II II •^'-9 ^ ^ H,H, k^N} ^ ^

m * MM cam

'^'^ ^4-. w w b(b MM MM MM MM DO .»»

• «» • « WOi <^<^

•^•^ •^•^ '^•^ •^•^ •* MM WW WW WW WW I *4>

(14U4 040. MM MM «K UU CMCM WW e««'^ «Q-» OiiM O^K ci^'*' WW*-

aia<(n wwtn hh qq wmin uio)« WW* WW* WW WW KW* (l«U4* WW4' Wm WW-^ Wm WW WW ww-^* wwCM

W pq PQm Wm

<t'*'*(MOO»0*tf)mOU>-*«OrOO'*NrOO'*OOCMCMrO<trOCM^CM OOIA m O wOU) *M 4»Otf) «>MOin 41MOIO «m «00 «m «

• «*M «M attM «M 4»«M <»««M «4t«M VM «t<»M «M \0«04wOw^^0w^}w<0Ow\0<0Ow^0^0Ow^0w^v0w\0w wwwHwHwwH w HwwHwwwHwwwHwHwwHwH

mmmKmMmmwmkhmkmmmMmmmKmkmmosmw

»»»(l4»(l4»»U4»b4»»b4p^»U4»»»b4»U«»»h»(h oinoiOMin 00 •*oo tf)o \o N oocM ro^

u) iomom
90

o

o

O

OO

o

O

oo

o

ooooO

oooo

o

oOoo

o

o

o

o

o

o

oO

O

O

oo

o r4 (M fO U)o COo CM ro in h- 00 0>

o ^ CM ro in CO

o

oo •>« »-t •H CM CM CM CM CM CM CM CM CM CM

n to

n ro ro ro ro

m

CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM

o

o

o

o

O

o

O

o

o

OOo

oo

Ooo

ooo

oo oo

o

oooo

o

O

oO

Qe

oA

oooo

Q

QQ

Q

o

o

o

S

o

o

o

Qoa

A

oQ

Q

a. CU GU a. CU flu eu 0* 0, cu Oh 04 0< 0. o« 0* 04 Oi 0< a* Oh 0« 0«

tn

M (0 CQ (0 (0 Vi

n 0) (0 0) CQ (Q tn (0 CO CO CO CO CO CO (0 CO CO CO CO CO CO CO CO CO CO

CM CM m •

CM (M I « H H

"< < M M

CQH

01whH

U

H

H

<< -< U M PU Hi «bQWH<

CO

CO

WU 1 04

M

M

M 0, Oh 1 CO * • * 1 H

£

MM

M

HW 1 w w u 98 a 1 (14 K 06 in 06 06 { M SI H m u •< H H z H M H 1 M CO CM «< •< • 06 CO (fa 06 06 m H u •-• 06 06 1 (fa (4 m • Z H (4 H

06

u

H

HH TO

Z

Z 1 «M •M a

CM

H 06 06 H H 06 RE

w 14

ZZ ( •<

TE

WH

M CH

o« W

U 1 ^ (1 (fa

H

•< -<

A 1 UI z W

w ^ (0

0«

K in 1 K

K •4 hI in 11 Ui

M

s

w

w • t• NT

WM • 11 M1 06n

:j

W

H

s 1 -<

t CM

Q 1 • U

W

t

O

O 1) Z CO

CO CO

t z w Oh ^ o CM

o in (1 H

w

o

O in • * 1 u z \ in

O

o

o • CM * 1 H

<t

XH

t

H

o

w vO «0 ^ X vO

M < H

H • * • ^ 1 H (0 • h" CM « ^ H» *4 1 M

(0

N «0 tow

t < o z • CO

w \0

w (d

uHH

t

u (4

w

H

w

H 1 04

g

ww

ww

P z • H (Pu

w

u -<

w 1 M H H 3 : M •J Hi M H H H H H a H 38 1 P w M M H M 06 1 91 w M M M « M K 1 o

(h 04

K

s

06O 1 Uh 06

OK

O < • M » h • U M H (fa » fa < z

IDMH

•<•JPU

(4>M

CO

(0

U)oup

CO

(4

04©

(fa W

0)

o CM 0> 00 «0

o

o

in ft CM CM

«o vON

u

u

u

u

u

u

u

u
91

oooooooo

(Vi(MCVjCM(Vi(M(\i(Vi OOOOOOOO

cQcacooicocflcav}

ooooooooo

o^(\jro<ttf)«of-oo (ninioiniointninin

(\J(\JfM(M(M(M(M(Vl(\J ooooooooo

(QUCQCQCQCaOlCQCO

oooooooooo a^o^wro^invONoo

oooooooooo

o o o O o o o o

o C\J n in «0 h- CM (VJ CVi fVJ CVI M

o o o o o o o O

o O a Q a Q

eu PU OU 0, 01 CO CQ to w Oi 01 CQ

COM

9

•J

w

fa

H

S

H

n

M in V

I (\j M • M

M • Cu H ^ I P U

O O B<

0] 2

1 M in w t o t-. n o 1 H

w

iz; 1H 04 S n

ro M^

U 04 M <^ w •J M M 04 M P

K w H 1-1 M a 04 04 w W H w S

-< 5 1 U 1 w ;« < •< 04 p V

u M M U M M P V 1 • M ID a M

pa n CO W W a M 04 4» O

u Q Q u a< 04 H H M 04 M Q M M o M M 04 w 04 04 w \ n M w ta H 04 w \ Ul (Q f 1 n 0. «-> M 9« * * :« 9C 0. 4> * ^ tu ^ n n Q« 91 in

<t 04 W t

-< -< >4 iJ ^ 1

1 • a n CO 1 41 Q Q h) 04 p p 04 to. M M M

« * w 1 * w 1 ro tl4 04 I-) •J M M H H • -< • 4. M i * o< w M t • 1 • W 04 04 «» « 4» 04 1 I-) n n

1 M • • 0* H t H 0} ^ *^ 04 w to. ^ *^ u) in 1 • H 1 n P n n 04 1 w w to) 04 H < •

• «»• M M H • in H M < < »J U« 04 -< a • 1 H to< to4 • • • o w M 1 o :« m p o 04 04 TD 0] w

< t o •< 1 o OS • • M O « OS 04 M M M K K

^* Q 1 ro ;s D • n • (l4 (VJ 04 04 H H 4» * 4t • 04 04 Q H o I O 1 04 •J • • I <-i w 1 • n nJ 04 to. 04 04 V 04 S 04 W W to. 04 S

0. U4 04 < •< H H p H e< H H •< m P H

u u u u

o

o

ro

o

CM u u u u
92

o o o O o o O o o o o o o o o o o o o o o o o O O O O o O o O O O O o

00 o CM to •«t in <0 N CO 0\ o CM ro <t in h- OO o r-4 CM to in OO 0» O •-4 s s oo CO OO 00 OO CO CO CO 00 CO a> 0> 0> O 0* o> 0* (^ 0> o o O o o o o O o O t-4 CM <VI CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM CM ro to ro ro ro ro to ro ro ro ro to O O O O O O o O o O O o o o o o o O o O O o O o o o o o o o o O o o o

o a a Q A a o o Q o a Q Q o a Q o Q Q O A ft ft ft Qu Oi 04 04 04 0. 0. 04 04 04 04 0« 04 04 04 04 04 04 04 04 04 04 04 0. 04 04 0. (0 tn CQ tn 01 (0 (0 01 (0 (0 cn (0 CA 01 0) Ul (Q CO (0 (0 0} CQ (Q (0 CQ (Q CQ CQ Cfl CO

Oi 04 1 H

NTS « n

•

ro

% 1 « 1 to 1

U4 W II • H 1 M M M M 1 H

* 4» n 1 oi; M

FI M

<^

CM

1 CM 1 PA w U. ?8 1 :« H Ui H « n 1 H

CQ CO I CQ o

m M CM V> • • 4t 1 « <t 04 H

Q »

»-4 H4-1 1 * 1 Id

04 M M ID -> r> H

w w M H

S

41 • 1 U

04 • A 1 n ro < < 91 n O < M M M I H o

1 1 O 41 * 1 04 m w m "•9 1 (l4 M H 04 At n m W w 1 n

w w 1 1 n Q in 1 H •J

(A CQ 04 • m w >^ w • t w

n pq M o -< w w w 1 ft

•< < 0< 04 • m H o r-4 Q ft H H n • ^ 4t ^ Q 1 01 •

1 « CM CM •d- •* * ^ 1 D n o

04 a( a« O o «0 « CM O CM o (M 1 >J m 0]

H <0 X vO ^0 X 1 O < t in

• • • M H CM 4» m • 41 41 41 m m ft I • H 0) o

< ^ D (0 vO «0 <0 ^ « <0 <0 w 1 M H n * o M M U <^ w w H w w H 1 * •< •

4» • • M ID w w M < P4 M w w W •< 1 • m H V pa n n n H M H H H a H H :« « w i-l H n

• a M M K M M M M M K 1 < w •< Hi V Oi c> OH W V 1 • lb H U4 H

-< n Q U M O » (l4 » 04 1 ft M M •J

U)

o

o

ro

00 (M

(M \0

N «0

(M

onw

H

tt)

D

^>

CQ

01W

H

w

04

01

Ui

CJ

01

•.J

o

H

01

I*9

ft

w

U4

ft

04

0«

roH

O. 01 ^

!z: w M

ro -St • ft w M 04

CM O

0»CM

1^

ro in o o o

* <

* oj

4 4l M •OH

»^ ft 25 O w • ^

ft U CM U

O CM

u u u
93

o o o O o O o o o o o o O o o o o o O O o o O O O o o O o O o o o O o O

(\J ro <t (0 00 o (VJ ro <t ID CO O (VI ro U) vO GO 0^ o •-I (VI ro IT)

1-1 •-4 •-4 «-« «-« (VI <M w (Vl <M (M CM (VJ (VI (VI ro ro ro ro ro ro ro 10 ro ro <t <t *t tn m n ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro to ro ro ro ro ro ro o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o

Q Q Q Q Q o Q a Q o Q o Q Q Q Q Q o o Q Q Q Q Q Q Q o Q o Q Q

(U CU a, CU CU CU 04 04 Oh 04 04 04 Oh 04 04 04 04 04 04 Oh Oh 04 Oh Oh Oh Oh 04 Oh 04 Oh Oh

0) CO CO (Q 01 (0 CO CO (0 (0 0) W CO 01 01 (Q CQ (0 01 01 (Q 0} 01 01 01 01 01 01 01 01 01 0) 01 01

COH

01 ro W (Vi

IDH

GO

01w

04

o

ro

V

o

(M

04

(VI I

(VI

a

M w

»^ y w (0

»H w

O w ^ pq *^

ro CO« (VI

Oh

ro

*

(VI

01 »^ H * » .-• W « U m

w

o o «o vO «0 «0 H

4t • «> f>« «0 ^ 4)

H

<;ai«

tl4

WWW

H H H

W

H

9f

O vO

• •

w w

H H

M M

m

* (VI

•-• •

vO

<; w

?! H

OS M

X) OS Uh ^

ro

(VI ro

«0 «0

0)

M

H

«<»4D

U

u

w

>H

0]

01w

SSI

D n

01

(VI

W M

O

ro

• 01H

M

u

u

Oh01«

<^ H

M M «»

•<t •t 4t *^ o •-J w H

o w w (l4 «» (VJ D n ro

01 01 in « « -< ro H o o u

p p Dii w

Oh © • • H Uh ro (VI (VI 01 H • (fa (VJ 04 04 M M ro

4t w * 41 • P •»

(VI ro >-> w (VI « w w o • W w OS M (Vi • Hi W M o w M (J P Oh H w M o

w u ro u » « b < < o p

o Q w 01 o o <; fi o

o OS 04 V M »-4 M w (VI •J ID •

ro V Ul H H •-I ro * • > •» (fa ro • w w X

CD • Oh V •< ^) w < (fa p

P H o U P u P Ol w M p

woo

o

o

(VI

o

ro O O O O
94

W w w W w w w o A A A A A A A A A A A A w

00 o CM ro «t in 00 Ok o CM ro in kO 00 o> o CM ro in «0 00 o 4^

in in in in in in in in in in vO *0 «o vO N N ^ 00 CO ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro ro 10 ro ro ro 10 o o o o o o o o o o o o o o o o o o o o o o o O o O o o o o o O 0

o o Q Q Q o Q Q Q Q Q o A o Q o s Q Q Q Q o Q 0

CL, a« 0« Ok ou a< 0« a. 0. cu a< CU CU Oi a* 0. 0. ou 04 O4 04 01 (0 W 01 (0 0) <n (0 (0 0} (Q (Q CO CO (0 01 (0 en 01 CO 01 01 (Q (Q 01 01

w

H

ro * 3 M

M • • M

• u. M

<M D

U) 01 f) •p Q 01

(M

n

»^ W H

^

* D ^ Ul M I • o ^ a

m

w

w

(h ^

W

a< • D

D CM Z

w a M

• m H

:* 35 >^ V

vj U4 O

CM

o.z

CM

CMa

M

M

w

X

1 «^

1 (VI 1 Oi

1 z

1 4t

1

. n

1

11

, CM

11

1

1.

J

! <^ ^

• toz z O4 04 O4

1 M 4t

I M ro

(Q U • ro ro ro w M M • 1 * U4 « M M

U4 ^ 1 •

(M n CM 0« CM CM (M

u n

H »^ » H «

at 01 w • n

n W ro 4» 4» * H • 0 rk z * H 0 • -< H W «^ 4t m 4) • M

ID 41 M M H

0 M 41 4» 4» • 0 * w 4. 01 0 n w W

1 < n P

0 • D w w w

Q (Q w w w

• w n

o

a

o

•

atD

(0

4» 1^

4t M W

<«>

O D

in (0

aip

M M M W

• P P

«^ z z

m H H

z z

w c o

(0 U4 u u

COn

(t4

CM

ro

V

H

«?

o o

4» 41 vO \0

W M

M M

K 0«

in

in

4) «i4 w H

W <

H 38 M 06 « V

00 «

O

• 41

w w

M M

A 1^ A A A A

^0 CM <» ^ < ^

(£> « CM O CM O

«0 X ^0 « O <0 • « 4» 4» 4t \0 w <0 «0 ^0 vO w H w w w w

u n M u

H H H H

M W

K 06

*

CM O (M

«0 «0 ^

<» 4« •

«0 «0 w w

M W W

H H H

H H M

K K «

» » »

(M

U) ^ ro in in

ro ro

o 00

in < o

ro ro ro

ro

in ^0

<0
95

o o o o O o o o o o o O o o O o o o o o o o o o o o o o o o o o o o o

ro <t in <0 OO o (VI to IT) vO 00 0^ o ^ ro <t N OO o ro in o

OO CO 00 00 00 OO OO OO a> <7> (J^ (7» 0> o o o o o o o o o o 1-4 ro ro ro fO to to to to 10 to to ro to ro ro ro ro ro <t «t <^ *t •<t <t <f < •* <* «t o o o o o o o o o o O o o o o o o o o o o o o o o o o o o o o o o o o

o o Q Q Q o Q o Q A Q Q o Q Q Q o Q a Q o o Q a o o Q Q Q Q

IX CU (X a< fX (X 0< (X (X cu cx IX (X (X IX ou (X (X (X PU (X (X IX (X (X cu

01 (0 (0 <n «1 tn (0 CO (0 0} (0 (Q CO CO CO CO CO CO CO CO CO CO CO CO (0 CO CO CO CO CO CO CO CO

H

H

<

D

• a

CO w

(X

•J<

n

H

w w

U) tx

(a

H

0.

n

0* o

M

(X

CO

o

Q

u

N4

H .J •J R

• U

»-) n

w <

O

u

w

A

o

V

n

i4

(X V

IXH

w

(0

i

D

2

M

H

U

fM

ro N \0 00 o o o

^ <t <

41

(X

ID«

w

(X

n

(X

roH ^ m

IX w

35

ro Q w

• IX Z

M R M ^ CO O M « o • o

•-4 PQ a

V w (VJ P u

o

o

H

COH

•QM

H

U

iJ

<!U

w

m

Oh<

H a; H

H >

OS m ^ 1 in

(X *** 1 >* ro

ro IX « ar» o

ro (0 % 1 o

b

* 1 < OO ro 1 » H

w 1 1 H

H M M 1 s

« 1 u

1? M A 1 M 1 w

* 1 n H

m m t H

•-4 *^ 1 M

• w 1 o; OS o U (0 1 ^? CM u w • 1 »J w »H

• n I M

<-» rl 1 (0 n as o

O (M to n 1 D W «) M

o o tO « (\j • 1 • »-« Oh 2: w to to to to X 1 o 1 CO K

H ro • m * * 1 M i>i to • H

CO to to to w tO w 1 u -< « P to H

m w w H w H I M u H ^> O W

M Ui W w •< W << • IX 2; u V o M l-t H H H 98 « (X i4 * H OO • hJ w M M M M OS • u »J ^ w

(X K K V 1 •< • U (X <D b <

M (X » (X 1 ^ M n M

ro tO <* in

to (M ro ro P to to N 0^ u u u u
96

o

o

o

o

Oo

o

CO

o ro

CVi (VI < <t

Oo

o

o

o

o

Q

o

a 8 o

tL, 0* 0* Oi4 CO m 0) CO 0) CO

o

o

o

o

oo * U> <0 N OO (7> <M <M tM (\J (M (VJ <t ^ ^ o o o o o o a a Q a Q Q

O4 CU CU (V«

h Ou

tn (a ui 01 (li (n

MD

asMHZ

<DUOo

00

N

DMMmH

»jD

01m

06

H

Ml«

ro« M

LMAX

n•t

HD

Q«Z

H

M « •• u* Q w H w z ^ 00 * •tt ci w Z W a H H H M 0( ID

pq

a • *

•Jw

X

H •< -< < D Pu 0 • z • 0 Ul • 0 M 0 s p H W 3« M u z Z • P • V M D P H

»J • w

Z

-<

b " • (fa

V MOP M 0 QZ

«0

«0

ooo
o o o o o o o o o o o o

(M ro in 00 o (M

o o o o o o o o o t-i o o o o o o o o o o o O

o o o o o o o o o o o o

n m n n n n n n n n n n

»

tn CO CO CO (0 (0 01 CO CO CO CO CO

ooooooooooooo

ro<(n\osaoo>o^(vir0'(»'in

•^>-lr4v-lt-(>-f«-l(V|(M(\J(\|(\|<\l OOOOOOOOOOOOO ooooooooooooo

cQcococococococococncococo

CMn

CM

ro

(0

to ^

H H

e

Q

TO

PQ

S

Q

H

CMW

-<

GbW

Q

n

n

H

H

w

a;M

H

P

X)

w

n

p

(0

CMn

91

ro •

CM \

^ <

• N

CO

•4

-<b

*o

in o * -4^ ro o \

w (O ro S w CI ro *

9 »

M (0 U

COH

M

(i

%

^ w

o in CO ro ^ J

<< O (h

« ro •

• CM O

O 10 O « It fO \ w in CM n ^ o

•J • « ^« \ fO • O w w

ro «^ CO w O

•< \

S u H H

U M

X « C5

91 a •< a ai H

ID << U U Q

*o

in

S

Hi

ro

10

H

8

w

m

K

o

M

• ^ St »^ •H

W

CM H

m «

Q

«< H H

*4 •

• O

M W

in in

M **" ^ M M

< pq o a

CM

V o m

o Q M

(l4 (l4 V M M (9

t • w w

p p

« K

H H

• •

• M

«^ <^ 53 ^ 25 • •» M

ua M H

w w 25

ro •-« ©

O (9 u

O *4 CM

o

t

IX, b

O 10 •

W M M

P P P

« » «

H H H

• • •

* M ^ 13

<^ » Z

• « • • M W

*^ H M 60 M H P

w www!? tLiCrO'^cMiowK

n

o

w

H Q

10 <t 10
98

oooooooooooooooooooooooo

0OO0OO000«if-i>-t<-4«i«^.-i<-4«-i<-i(M(\i(\j(Vjc\J(\J ooooooooooooooooooooooooo ooooooooooooooooooooooooo

o ooooooooo

ooooooooo ooooooooo
o o o o o o o o o o o o o o o o o o o o o o O o o o o o o o

ro in « CO c o C\J in 00 o ro * in CO o

o o o o o o o o o t-< •-« tvj CM (VI w <M (M (V) (Vi ro o o o o o o o o o o o o o o o o o o o o o o o o o o O o o o

o o o o o o o o o o o o o o o o o o o o o o o o o O o o o

o o o O o 0 o o o o O 8 o o o o O o o o o o o c O o O o o

^ as at :« a as as ^> 10 10 V 10

ro

M

W M ^ PC ro * w ro § w

ro

ro

ro

ro

O

91

-< 1?

•J

CQ

W

o

ro

• •

n z, < ^ w 0*

CO •

< Ui U I-) >

M \

o z

M W

h) 91 0. 91

at O

>

« ro « ro st W w P

•

M » a

© ro

ro ro 4 H

« w i:^ SO©

o n • • © • w ro

J5 ro w w ai H at o

a; \D n v\s «^ (M ro

<< < ^ U O V UNO)

S 3 S

N \ \

:k z z, O O ID X ai 91

91 a ai o V V

u u u

row

* M

0)

u

<

M \ SK

Oi

aiIDN

8 w

•H N ^ (\|

o

o

gM 0) 9« a

H Ul ID P

U Z Z H (O ID

ai ^

H V M

M * >^ Z « Q

« at a

ro ©

• OS M • 91D

O M 01

w

Z

• M

4^ M

04 Z W

M O

f 8

M •

0^ Q

z 3

m «

M •

P

91

asD

U) (VI M

91D

0}

a Oi H

D Z w

Ul 4) P

Phi

• ^

^ ro H

o

P

o

p

at

Oi ^

Z I m M

ro

• 98 H V

I

ro

p P a

© 91 P

P V wi

P

as p o

at pa p o a p

(/J P V

CM

a p w w ^

p © a a a © a

CO P ©ID © P ID

ro

p a

• ©

^ *

« o

p

ro 10

04 f z

w H

a a © ©
100

o o o

(vj ro n n ro o o o

o o o

O o o

a« s s

ID C V

o o o

4^ U) «0 ro ro ro o o o

o o o

o o o

» :m ID ^ ID

O

^

ro ro fo o o

o

o o

oo a>o

o

o

o

<to

o

O O O O

a 34 ID ID ^

a.

o

Io

w ^ WW

w M 3« :« C O w C ID

I a X in <^ ^ ^ ^ \ rt ft O

pu o A« n

w •w tw vwfO

Q

10

• W

M w Z OMwMOMfOK •^wUwtnwHP

Q S S SB H P

o o

in
101

o

o

o

•JH

O O

W fO O O

o o

o o

H H

O

o

o

o

H

o o

10 ^o

o

o

H

o o

CO 0^ o

o

o

o

o

•JH

o o

<vi ro o o

10 o

o o

CO o> o

o

o o o

(VI ro •* o o o

to o s

o o o o o o

00 0> o ^ (vj ro o o

o o

o

o

_)

H

o

o

H

o

o

H

o

o

H

o

o

•JH

o

o

H

o

o

•JH

o o

o o

•J ^

H H

O O

O O

•J

O O

O O

H H

O O

O O

•J

O O

O O

•J

O

(0 ^.H^^^^«-i^_i^(M(vj(M(vi(\j(M(\j(\j(vj(\jrorororororo

o

o

•JH

HH

06W

H

ro

to

roo

o

in

o

m

(M ^ w H

W H

«

n

•<m

H

(Q•H

U

«

M

U 10

Oi ^

K w

• H

fO • w 10 > »^

-< ^

o

ro • w to

n

Q

X>» tx)

H

• n

CQ

Q

a,m

to

ro

ro ^

0< w

N 2 a

I • at V M ^

m Vk m B » ^

* ro

s n

» \

w <

< M

MVS

u z :c M © ©

91 31 X 91

an c (D M y o

t

00*<

m

a*

ro 10

U CO < <

w

z z

•a

0. * Q

•< (O

w

* I-) Q (fa •< K

a (0

« w

tf) * •-• n w

ta> M \ H W \

J H

«< M

D

«

r) O

« Q ^ w • ro <M «o

•«* z

in • *o

w 10 Q

(Q ^ •

(fa w in

to o

ro Q

• «t • in • o

in

I

IDZ

H

(fa ID

(0

SBM

H

(fa

(4

«>

w

z

«

(D

•JM

(fa

<!O

(D

<

m

• Ul M

• (fa O V

«w Z

ID

(fa M

04 H

O) << z P$

w o

pq (fa *^

(fa w

^>

> V

z m

V z

a at 8

UM

n

<

Iz

IDZ

(fa

CQ

(D (0 ^

(4 I

(fa W U

V X U)

(X«

•J

(fa 01

u

w H M M

M Z Z M

M X «

(fa

(US

(fa

z

(1) n

•<M

Q

IZ

z

«<

(fa

(0

(4

h)M

(fa

(fa Z

©

M

H

U

s

(fa

(4

0)

(4

H

•<

h]

D

O

«<u

z

ID

(D

z

H

(fa (4

u

X

m w

M

CO u

M (X] n

O

CM

Z

TOM

H

O

«<

(fa

(4

w

D D • (0 8 w w (4

w

^ ID % z (fa (fa (fa n * « (D P

(fa m

f-i (fa (fa H

• H 0 M •<

(fa o w (M Q w (M P «» n w H (fa w H

< t Z X z Z z

1 1 m M »9 H

w w w M

w M o S (4

(fa w M X o w M e

(fa (d R H 1 • P Vt O (» X o M X o o Z M tn o Oi • ro w M > p H M cc M ifa • (fa • H H

w (fa [fa *-i z M

ID (fa H M > H n

n Q uu Q p p (fa u

o

ro o

(M

o o

<t o

u u u o
102

oooooooo

oooooooo oooooooo .J i!) iJ h)

oooooooo <ir)«0SaD0«O^

oooooooo .wOOOOOOOO

o o o o o o O O O

CVJ fO in CD O

tn in in in in in in in <o o o o o o o o o o

o o o o o o o o o

•J >J _) »J h) >J H H H H H H H H H

M

Q

H

o

•J ro

ro U ro M ro td ro H ro U

•« H CM -< H

•J

• W 2

M •< 0. o o U

u

u

o

u • U)

o • < o * H H

H * A n

2

o

u p m

w

-< m

H

H

o Cb •< ID (JQ

OS »

• D

H

< «< w

M U

• cu (b M ai

[1. w u n «

h1 H ^ m

m

EU w • GU n 1 CM

I-) Z H Q ou -< M »4 H

M Z < «? « o • H f4 n H *^ H hJ > H H I 2 •< M (VJ © H M u w H H w

o

CM

u U o U

o o o ^ CVJ ro o o

o o

•J I-)

o

o

o

o

o

o o

in o

o o

o

o o o o

oD a> o

vO »0 «0 o o o o

o o o o

nJ »J I-)

• COu

•JM

\

©

H

Q

U

H

W

z

u

w

<<

COu

»-) M

h

a*z

v4 H r-t U

n •< m W n

M H >

• Cb • • <

C\J M U CVJ a« o M •J ou O > o

z •< 1 z -< < < M

« M O a • K •

(M • CVJ • • • I

• M o "9 s; Uk Ul

(0 o I a • I D

o o o *^ o o O > ^ o Z M in • M • 1 oo •< (;> M

tj o * H

w 1 CVJ « > »^ Z

cu z > > w

CO » a iX X n p Q » Q Oi 8

O o o o O

>0 in CO (;> r-

>-< o4

O Q

• •<W

CVJ *

ft. z

• o

m ro

b Q

CVI z

•I w

H

U

«

CVI z

a CVJ z

* «

CVJ Cbz

m n

b

o

ro

o

H

<nQ

I-) U

u

n

H

Z

»

H

O

w

H

n

(0

CO

u

QM

COH

n

M

CO

»

n

H

TO

OSb

QM

H

-<

U

>-) «<o

CO

«

D

t-t

U

b

w

H

Z

M

H

U

•<«b

COCO

u

H

H

CO

w

n

H W

CVJ Q

Z b M

© Z W

b

n

H

Z

M

CVJ

w •-> »

09 w ©

H b z

<

n

H

W

H

©

Z

z

©

(/J ©

b

a • © w

M

P w

z

©

b

h) w o

•< n p

m H o

H • z z

o u woo
O

O

•JH

O O

tvi ro o

o

H

O O

in «

N

o

o

H

o

o

o

H

o o

00 0> o

o

o o

(M ro o o o o

OO 0« o

o

o o

(vj ro o o

in \0

o o

OO (7« Nooooooaoooooooooooobo«o>(7«o>0(>(^o<(^(^ ooooooooooooooooooooo OOOOOOQOOOOOOOOOOOOOO

HHHHHHHHHHHHHHHHHHHHH

O

O

O

o

•JH

w

41 H

9 ^

n ro

9

o o *^ O O Q

Q O O O

O • • ^

• • • w ro c

CP u Q o

»9w

H w

^ n. w w

m ^

-< w • N <

u

m

*U

»U

«

• ro * o cvj

roa

o

PQ O \

roro » ro

ro w

P) • o

ro W »

• H

« w »

ro U4 p

n w ©

« R u OHM

*H

w

w «

I H

(\J PQ

PQ N

U

w <

PL. b

U

roCM

P

M

w • < H

P

H ID

* U

H M

fO

(M

V

H

M

U

H

H

iJ

t H

D

u

M

m

H

o

SB

-<

u

•Ju

(Q

H

M

U410K

W

«

P

H

U • m

W M

H

H

I M H

w <C b n

M H

W

a.<

n

IU)

0«

Z H

CM

o

too

Io

p

•

Ip

•

i-i M

a«o

M

H

M

«<H

<

<

p

Q M

ai C P

O

H 0)M

P

o

-<PQ

H

P w

H

o

H

I H

H CM » CM • n Pt H u

w • H ro

• • M pq M

• H M

o • M 1^ w H

Ul p (b

o O

10

p

H P

H is

u
o o o o o o o ^ C<i <t U) O N

o o o o o o o

o o o o o o o

o o o o o o o

Ai Ou (X4 Q4 CU 0« OU Z iz; » S!; SS as 2

0000000 00 0> O ^ (M fO <» O O ^ «4 0000000

Q4 CU fti &« (I4 2; Se; Z Z Z Z 2;

n

m

o

o

1*) N

roo

• z

• • n t4 b4 •< I M « V

• z • »

H (h i • w

w ro

CMw » 00 « 0

M * in in Z ft M 0 I M "< •* in n

U4 z •* 1 00 • 00 • U) M 0 « * z z in • • H H • I • • w • H M M <M M N It p *^ (M

D OT *^ <^ z z

n <C> M Z Z 0 0 M H M CVJ z K

"<i OH •H <VJ • CM • D

n a, ;« z H Q

n D M c> w W Z

M VI M Q z n p a U* (J Q K w

VI• 0 t-< 0 CM

CVJ (Vi CM h4
105

o O O O O O o O o o O O O O o O O o O o O O O o O O o O O O O O o O

CM fO in 00 o «H CM ro in vO N oo 0> o CVi ro -* in <0 N 00 0» o rH CM ro «t o O o o o o o o o «-« •H t-« •-4 «H CVi CVi CVi CM CM CM CM CM CM CM n ro fO o o o o o o o o o o O O o o o O o o o o o O o O o O O O O o o O o o

o o o o o o o o o o o o o o o o o o o o O o O o O o o o o o

2! SB as z; z z i z z g z z z z z z z z z z z z z z z z z z z z z z §

M 06 « 06 06 K « K 06 » K w w 06 « 06 06 w 06 K K 06 Oii 06 « » 06 H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H

o o »^ o p a

Q o o -<t O t • • a •

• CM ro ©

m u Q a

M ^

*^ *

f> w

^ PU ^ w

• « [fci H W

» OS w • m

H >

» •< <<

I H »-<

CM ro u p

N R CM ro u p

>

•<

a

o

p

ro

roro o ro p r) M

in ro *

• • (36 « * 9« CM ro V a

a a U P
106

oO

O

O

OoOO

oO

o

O

O

O

o

oO

O

o

O

O

O

o

O

OO

o

O

o

O

OO

O

o

O

in CO a>

o ro IT) 00

o •-« (M ro in <0 oo 0>

O •-4 (VJ ro <t in ooO

ro ro ro * ^ < <f ^ in in in in in in in in in in vO <0

o <0 <0 <o <o

oo

O

o

oo

Oo

o

o

Ooo

o

p

o

o

o

o

o

o

o

o

o

oo

o

o

o

o

o

o

o

oo

oo

oo

O

ooo

o

o

o

ooo

o

oo

o

o

o

o

ooo

o

oo

o

o

o

o § z s z z z z z z z z z z z z z z z z z z z z z z z z z z K K K 06 06 « K M w w w K 06 K K « OS w « 06 K 06 w K 06 » QC « K 06 H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H * •

(VI OJ Z Z w 'w b (fa I*

(M (Vi w

(fa < w

(fa «n

H

(fa

M

(fa

roM

10

(fa

(QHZ

ID

a*

(VI (VI

H

M < 06

M • O M

w

(fa (faH

H

H * m <5

(VJ (VJww

H ID ID (Q (fa

H • * ro < n (fa (fa o a

M

m (VJ * « M 06

H

w

H «^ (fa (fa OS ID Z ID M M M (fa (VJ o n H H * « o • • w (0 ID * (VJ (VJ •-» (/J ID ID (jg Q « H ^DO

LH

n

ID

NT

m

H «

ro

MM

Q (Q

w

(fa

H

M

a

M z ^D

t • > M w Q (fa ID (fa w ID p • O H Q n < « « « Of; * (fa H M u * >0 w H (VI (VI w B3 CO (VJ w U H V « W Q Q 06 D < • Q o (fa (fa ?. H > M o < n H r*- m X

(fa » « (VJ (D (4 :«

O < w * • Q • * O 06 OS Ifa >^ (VJ • P4 OS

m

o << 1 z • • QP oo

MM z O

U (VI (fa

o 0« •< • W « D o w n ID (U •< w • • o « (fa 06 M n V • H • • • (D H M m * -«t CO • • M • H • < n • << • ro n eg (fa H M M (^ M e • o n \ M M (H >- U •-4 w 06 • • *^ M M n • *^ o > P M * Q a • o • M w z H m n (U • in • 06 H w "< H O • ID < H H CO R ^ m »J (VI at m Z (fa (fa (fa H n CQ M M M O (0 Z ;>< w M r-« eu (fa P • • « p M -i M WJ (fa Ui CO M M (fa (fa • (VJ (fa w ai w w w w Z w (VI (fa M Bt (U (U (4 u (U

(A

w • 5 < • u (fa

W Vi (fa (U

m (fa oa (fa

n Ui (fa Oh OS OS «< OS

(fa

H

M (fa Ifa

M

M (fa (fa (fa (fa

M

o

M (U

H

M

u (fa Oi

X (fa

o 0.

ro

(M

(UpzMHz

CJ

(VJ iH

ft

(M (M

O

O

O {J

o
107

o

O

O

o

OO

O

O

oOO

o

ooo

O

o

o

o

o

o

OO

o

oO

o

o

O

o

o

OOo

oO

o (VI ro (0 «0 CO 0>

O >>• (Vi ro in 00

o •-I (VI ro <^ in vO 00

o (VI ro in

N ^ N

N OO CO CO CO CO OO CO OO CO 00 (^ (^ 0» 0* 0>o

oo

o

o

o

O

O

o

o

oo

Oo

o

Oo

o

o

ooo

ooo

o

o

oo

oooo

o

o

o «^

Ooo

oo

o

o

o

oo

o

oo

o

o

oo

o

o

o

o

oo

o

o

o

o

o

oo

o

oooo « z z z z z z z z z z z z z z z z z z z z z z z z Z z z z z z z

P6

K (US fkS

K » « K 06 06 06 « K p< ptf

K Ptf p( OS P4

K Oi OS OS 06 OS Oi OS OS OS

HH

H

HHH

H

H

HHH

H

H

HH

H

HHHHH

H

HH

H

HH

H

HH

HH

H

H

H • M nHzH • P

(b

M

MK • P A • P • 5 H 1 O W P o <

CM

W xn h) oi

M (0

p •<

w

p lb

& « w (0 (VI u pi! (I) ro

(K -< z H

•<M <

(0

H

H

H CO

u

H CO

w

M z ID

M

M

H

M

P

O

ID

H OS 9 W K U M H W Cb H M <<

Q Ou lb i-l •< H ID • P« H TD p o w Oi * M Cb 01 06 o •< PU M 0* H

;«

H "< •< • H H P (b W H h * CJ H

•<

K )_)

o ^ •< 06 X> \ Cb P (0 z

(b Cb

o (Q (b

H • M (0 H M Z < z K iz; <^ \ ^ w • (b -< H K ID 06 Z M p

06m

H

M P< z

CO

H (M 06

M -<

M

m « * H 06 lb • -< H u O (M Z <, H • lb •< < •< •< • H -< W O Z (0 U ro <D tD lb (b (b • H Cb • w > ^ 9 SB 0 «0 * * M n Pei p <D to U t <^ *9 z 36 z iJ H «i wt 1 I 1 ro >^ H w »J z 1 M Q n p lb lb CM <-« CO in in lb ro w M H (b M Oi O w X o in o w • (VI K piS lb P CO P O o Q in • in K ro fj (b lb * H p • • >< M z < i4 • • • • (b • w a • <^ A z •< lb m w P • H» . •< ro • • 1 4> lb M (b CO CO H c

GU

V

H

M »-•

w

MM

M • • • • p > • -} M M M M • w w (VI H P < •< H M ro • OS

w

w

s * • « •1

nn PQ in

D

o

w

n

M

Z lb lb

M •J ro

H

(0

n

H

n (Vi in I Z «<

w < 06 CO >* • • • • • u w w w w <D H p I © M • H V II UI • • M MU lb <M M lb H H lb lb a, -< n ro ro M w tr* H >• w W u w u :m U Z «< w p 38 u z >^ w • e>

PQ (b «<

M « K ID « « tu PC (D lb >• lb (VI ID

w

M < (b IF

V -<

M (b

Q PU (b

Q lbu 0.

O

u

M

w

M 06

U (0

nH

H

Pp

H

H

H

ro in (VI ro

(M CM (Vj

OU

u

u

u

O

U

o

V

u
108

o

o

oo

oo

o

oO

o

o

ooo

oo

oo

o

o

o

o

oo

O

o

o

o

o

ooo

in en r»\ V CVl rn1 # 10 00

V

w (VI >* U) if) 1 OO

V (NJ ro in CD

V

Q

O

ww CM (VI (VI f\l (VI (VI CM (Vi (VI ro rift

mD rr\ 1 f ro CO 1 9 ro ro

w

W

noo

AW

w

ooo

Q rto

w

O

o

o

w

as 2 a; 25 X* 2; 25 25 25 z 25 55 25 55 25 55

Z z 2; 55 2; 25 OA 06 OS 0$ Qi W(K OS W K 06 K Qi K Q6 06 K IK w 06 K

r'

h ' MLrf IT*

L • H

H Ih ft

H 1 * • n K

f/l r . fti (b z w TO FN M riw H M rr1 U « 111 Tn

(h

#

MW r«H< n [Li w r . H fii W tb

ITS V

m

U ft rnwi

tnp*<

K frl [H

H Q. ft

n

4i

u f \

w

M P*4 (ht ' ft. in w [_j rn 7h

r«M rn Lii

Q rvi

p*i

ft.

W

M PC

M p\\J

H tlx NTM

r/}

w

n * z U

<D • w M M z n

w (VI [Li . 1 M w

m

W PP

m

06 fr, fVI VM Ed

ffi < M ft 1 Bfl C_i n m o l-M r . r* n V M n (VJ p a ITS V pq

ff1 HjWfti (n

W 1 m lb (VI

(K4

p in ttl . 1 H 1*1 M M

ID »«k M p*i vJ

Z -<

H

w

n

rW1 <<

o ft. H%

w •< > o?Vi

H ft

M1 Q lbw

[L|

M

M

OVV ( (b « M rj

»-)0 rn

w p#4 t—« M fti rj

H >u %^ fm 1 Ok (^ ft. I 1 P

M •—* • rn % HJ M z V]

f)

H

W

M -< (0

VM

m 2 W [Li 1*

r\ W [^ in

H 1 Pi'

Z rn

Q UJw

to

oQ [^ (VI 2;

o

m ft, V3 P (b rv <• ID H M CO R R 06 -< p

H

P [tl

W

t U)

#

# 06 PM

M (V,w

M • • M lb

u • ap*i

H ro ro IB

H

M

D

D

H • O

H

%

H

w

M 00

u •< • O * • • Z Z n W cn w z z CO P o n •J

H • M V <! • M M M n oo • t w M • M in (b H P « M H W M H H z H ro ro •0 lb (Q V y

>^

W • O

w

Z

Z

wWu

ww

M z w

Z

M

o

UW •J •

tbQ (b

n •<Q x>

©

n pq

M z K

V (b U4 06 T> « (b ID

n

y • n

m -< (VJ

M

Q CP

M

H

H

Q

H

H (Vi

Q

M

M Pb

u lb

M

yH •< a.H t-t

y

H

U) vOo 00

n to ro ro ro ro ro ro y y y y y y yyyyyyyy
ooooooooooooooooooooooooooooooooooo

^(\jro<i'»o»ONooo>o»-«(\ir)*in»or^oo(7>o.-«(\jro<*^»f)«oscoo>o-<(MrO'*tf)

o o o

06 K OA »

H h

O o

7i o o

H K W K K

H H H H H

fOH

* ro ro H

1 H

P

(M 1 r4 1 p

p

in oo

• in

ro

1 •

(Vi • (VI H

* (M

(VJ in (VI H

H 1 H *

« p *

in (

1 «0 t p

in P N

in in o

• (0

h- • CO

• • ( * ro H

I • * H 1 fVJ * P « 1 ro O P

• < P ro p ft CO 0»

«0

• in

vO ro • fO n w

O

• P • CM

1 in P pa * o s ft ro p ro «0 * 1 in in ^ vO o

• f-i in 00

00 00 P «

• 1 • • 00

• NO

« • • • V

<M (M <M •

<D ca h) »J OS »^ W w w in 8

X

IDU 8 p u

H

WM

u

H

P

w

H • -< W

»-) u

D M

ss

"3

H D

(I) (Q (0

(fa (h

«

• H p4

H

P

• P

H

M

m

o O O o o O o o o o o O O O o o o o O o O O o o o

z z z z z z z z Z z z z z z z z z z z z

06 (K 06 « K « M (K K K M (K (K 0« 06 (K

H H H H H H H H H H H H H H H H H H H H H H H H

• ro ro U H H

ro H H H CO i H « « M M (d

•-4 » in o r-i Ul m » w

1 1 (b H D

1 P P Z « tD -< p 06 M (fa OS (»

ro ro W h) (4 [fa (0 H -< H

CM U H H

In p H < •J n

ro • • M Z h] <

• » (4 p

-< O (0 t Ui * CM CM OS M z 06 CM H H w < n (4 m

H 4> * H U H o

* N 06 H U) o

1 I (4 M

I P P Ph H • H Q« • P CO < CO H CO

0^ ro * z (4 (3t; [fa H

CM «o < >> M (4 M CO Z

«0 00 O 06 H o H M M

<t r-l P O H M (4

O • • P o (fa (4 Z (fa P • ro in <t hI < CO z ID * • 1 in -* (fa OS M P

in 1 H CM ai -< m D H M

H « * • N w o W ID « Oi 06 « <* ro •-« • u H CO 06 TO O

• 10 t • H « H m

1 P P H « H CO x > CO H

P CO N « M M ui M (fa M -< M ^ in in • 1 «t Ui n h) n

« 1 06 P K m • < H

in in CO p H in H H P M n H

CM o 00 •-J o M H M

• • • cn p H Z m • >

w in • • M H •< U] w P M

• X z (U H

o o « * P M D w o H • y

s p p o o o o u U4 U m p <i D

o CM ro N p p p O p M H h) H H P

P CM ro o in P oo O n < •<! M D M <^ Z

O O ro in o ro CM fv 0^ P » » (Q O U H

t P P (M o o ro • 00 CM • ^0 D h) < « o

in ro ro CO 00 o> 00 ro • U Ifa h) < (fa ro ro o CM CM o o vO ro hI • CO h) u < P h) O

in • • • • • ro • o> in M P M u •< P

• • • • CM • • • • (4 W CM o

> P (D U 06 • CM (M CM ^ VJ (li • z w m •

(NJ CM tD m Z CM CM CM (M CM CM M P ID u * n

»J i4 V n Z ID n n (fa (4 Z H H n H M

n4 m ^ »4 Hi Hi • n H

« • • V V * * <t • • « w U) -< CO (U n (U n

(M (M CM CM CM CM CM CM CM CM CM CM 06 « h) (Q h] M P4 • 06 O V V n n m m n n Z » H D D D O

n) »-) 1^ t.) h) C/J h4 U H M t-" (fa X QiS 4

hI •J X (U

P P P p p P p p p P P p H H M < Ui M (fa M

P P P P p P P p p P P p as (fa u ;« (fa

0«
no

oOOoO

O

O

O

O

o

O

O

Oo

o

oo

oO

O

oO

O

OOooO

oOO

OO

o

CO

o r-« C\J ro •«t in COo CM ro U) «o 00 0>

o CM ro If) 00

rN 00 OO 00 OO 00 CO CO 00 00 00 0> o» 0» o>

o Ov

o

CM

o

CM

O

CM

o

CM

o

CM

O

CM

o

CM

o

CM

o

CM

o

CM

o

O

Ooo

o

ooooo

o

oo

o

o

oo

o

o

o

o

oO

Oooo

O

O

O

OO

o z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z K K w % « 04 K OS K K K « w w K K K PS « K OS K « H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H H w*

CM

M ^ in « H ID O H « • « Mp

a ^ M M M M

04o

*^

o • • M • M ^M O M

10

w

Q

W

0<

O

M

n

M

Z

w

yM

M • M ID w

;«

w« O Q«

n

06• in

pq

•-} 1© w H

0«

«H

MH to

CM

in

MwM

MD

MV

P 06 • w

O

n

0<

M • HH oQ M

10

u

ro

W

ro

P

ro «>

ro CQ ro• >* «

CM

a

H -< * HH o U

eu in• OQ « • in H o (0 o • u Q Z

in

M in

M • • M

M < * • •-4 w • H P n

CO

Z

M z «-4

H

P 06

g

b

M « H m M w 06 4t

M

0« • K U4 < M

M

0* EC

oo

ME 04 «

PE

H

JH

n ro

P 04 Q4

Z

ww a;

H Hi

w

M

M 1X4 [I4 OS

Oi • 00 « w

W * o w H > JH M S3 in M PR PR u w M * «

•<

Qw

wH

M*

EF «^

b

IT

H

w PY PY

*^w

>* « H 06w

H

H hJ (b

Oi * O z 04 04 04 HI 06ZZ E3

H

w

h)

n • z z O4 •< UlZ

Z 0«

<< •< • w • * w

P

w

u

w 04

B • h CM CO CM • * • • * H w i • H 1 (b M • • in 0 in <^ n U4 H • M M w t-l • M M M

M

w *^w

M

p • PQ • CM «»

P

w «>

Z -< U2

ooo

m < JH

Z

w

o

n

M <J

Wo

M •-4 w-t z •J O4

w

M < h- • It" M P CM H4 CM CM (I4 lb (b H H Z w h4 m U u u w H u CM Z Z (fa •< at K 0< n bu Z Od n w u M 06 04 Q Q 04 H M 0 w m z

04PH

Q

W

Z

CM

in

o

in in

o

o

oo

CM

00

in

o

o ^ o

CM ro
ni

o o o O o o O o o o o o o O o o o o o O O o o o o O o o o o O o o O

(M ro in o oo 0> o t-* (VI ro •* in 00 O •-« (VJ ro in «0 00 0> o (VJ ro o o o o o o o o o rH t-« •-• •-4 (VI (VI (VJ (VI (VI (VI (VI (VI (VI (M ro ro ro ro ro o o o O o o o o o o o o o o o o o o o o o O o o o o o o O o o O o O

o o o o o o o o o o o o o Q o o o o o o o O o o o o o o o o o o o o

z z z z z z z z z Z z z z z z z z z z z z z z z z z z z z z

{J u u u u u u u u u o o u o u u u u o u o o

« w « 06 « K W K 06 oc 06 Oi Of$ 06 » 0« 0« K M « 06 OS 06 «

m o

>* ro Q w

0

o

•< o

ro ID w

n (QO

M

• o

ro O

« ro o o

H w <

•< <

tJ * o b

d. w 06 H

an • fa w

• H o fa (M •< -< ro X

W fa in (I. * o M o p

ro r) m « n • w M

Q H ro fa w

o M M

• • << w

Q o »J < •< • • ro o fa H H

ro • w ro w W W

•-• -< ro o 06 • PQ m

ro ^ w H 4i 41 • w O H •< ID (> « ro • in •5 ro •< tJ o Z ro

ro (VI H fa ro M in • M •-4 M M \ « w U H • fa <r • w *

S5 1 H • (Q n y • w 06 • P H H

U U • < ro H H H p R H •-J H 0* W <f O V "9 z fa w

< n W o w * fO « 'P -<? tD M (0

M 1 • M ro ro « o w « Ui 5 • M fa M fa

• -< • m w ^ ro fa • p w • w

w (VI « O O w "< ro * K 4r 06 4i (M iz; a « in ro w M H «^ U oo • n e w H fa w p fa M w M M M M M

H « W \ (VI •< n z H w U w W w w

s >< « \ ID CM \0 o (VI o (VI (4 O (VI Q (VI w o

< » © \ K U fa oo fa 06 < fa fa « fa K fa fa p "< n u U W Z M fa z fa z 4i fa z fa z fa U K < « n H "< * (VI • fa pg « • fa • N fa « • • PQ > -•J H •< Z (VJ • fa • (VI n (VI (4 (M (VI w M h» M w O « fa <U • • • (VI • M o • in • w M o • M o • w M

fi H M •J \ N Wfo \ VMM 1 M H M H M H H M M H • o D u •< z z w Z z (a z M y M M

n M U V v> V z «u M O z • •-• w V z w V z o 00 •<! w

04 M ai ai a« at P s as w w • H W H t H fa fa • M fa fa W

m ai at at ai H H 1 H 1 H

n o a« n V V M ID fa fa •< V fa V « < fa O •< O • << M w (J u * u U P P 06 p M M U p M O M O p w 06 O M o p K p fa «

• CO Ot-«(virO'«riniocoNO

W
oooooooooo

oooooooooo oooooooooo

oooooooooo

OOOO0OOOO (n\ONoo9>o«^(Mr)

ooooooooo ooooooooo

ooooooooo

oooooooo

<tinvOh-oooNO^

OOOOOOOO

OOOOOOOO

o

H

Oi O

z • •

* M W

M « p

M Z

N w M

« UJ H

H Z

« O

<0

IDH

O

W M

o

w •<

m b

• • V

O ^

<0 •

Q M

w

(A*M

w

OO

u

IM

n

o

w

I

H 00 *^ ^ (M

CO Z ID - O MHO wr w O

0*

(VJ 0, •< z O X

• WW

(M I «

• I

<M

O

z

«>o

w

M

w

w

H

CM <

Z •

z

o

w

H

•M

w

W

H

CVI •< z •

(Vi w

\0 z

O N o

in vo

w W M

M M (M M (MM (MwCMO aO* • O .^(M'4MMvOMM\0

HmH»OWH>OWH

H

Z

o

w

\

M

w

w

H

CM -< Z I

CM • O

M «0

O

Q

CM O

• z

(M

« • W

in K

vO W

w Z

H

•< ©

w o

\0 w

H

Z

K

O

H O

© W Z

O W W

z

COCM CM <0 CM it CM <0 « in CM vO

-I «0 O

«0 vO «0
Aw

o

oo

oo

fVI rn ^ to

V

o

o

o

o

o

o

Q

o

o

o

o

o

Q

Q

\^

oQOQo

W

HHH

HH

HHH

25 35 S5 35 35 35 J5 55

m rtim

m

m

m

m 1*4

a

hH

b)

1*4 ro 1 p^ w

(M

Q

a

CDM

04

ro

n Stw n ^ m s

rv ff) ^ S

H

PE CM •

ITSm

ro (M • m m1 f ^4 ^ w ro If

rn

M Ll] (VI «»

n

ro < • ID * rn

t-l

w

ii

V 9 ro

<f fL

H in ^ tD a

(0 ai 1f • H

a • n

w *^ • ID ^ pM

PU • • to 41 as M fid ro 10 w • s ro w -* ro • ^ w « ^ b w M w w •

(V)

M

a OS

m to • n Z © w n 1 ' ?i :r> \ wV \ W \ \ (VJ ^ "/ ro H(

UJ en u] < w CI) •< U 35 V Ol. Hi* m < < M W < w * n a « a < w M «p V H H V >- Q H M J \ \ \ N \ \ \ \ U P u < Z 35 35 35 35

pq

M

V <o

<!

M ;« ;«

a

a

aa

a

n

00 Ol

e

aa

a

aa

P

W

M

M

Mu 85

u

w

CQ

41 HH

(I

o

oo

oo

o

o

o

o

o

o

<M ro <t to oo

o •H CVI fH tH l-t •>« fH *H fH (M W o o o O o o O o O O o o o o o o o o o o o o H H H H H H H H H H H

35

Z

ZZ

Z z z z Z

Z z

www ffi tii 111WW fllww ffi w

W 111 w

«H191fM

hTm

IV. s

^1 ' hri <^ M1^ m^ Q^ 1 / r a ^ ^ Hj

fv rn ro

m

1^ in U/ (F^ u/ rnw ^ % fH < w mP*4 ^ m rv h4 H4 w m>*4 [LWW V^

0^ ^ V © LlJ N. O OS n.••4 « z H < CO 4t o

w 2^ Uw <!^ ro (VI w >i a o H 10 z 04 \SN M a CM CM 41 z 4) z z z o p «-4 * V c • o • a a a • ro • • CU (l4 a a a a a a I • <D V O D p p • n • u u u ta p 01 w n m o

ro

113

o

O

O

o

O

oOo

o

O

o

o

n < to 00 0>

o C\J ro «»

(VJ <M CVI (\J <\J <\J CM ro (0 ro ro ro o O o o O o O o o o O o o O o o O o o o o o o o H H H H H H H H H H H H Z Z z z Z z Z z z z Z Z M M w u w u w u u u o

If)

o

oQ

10 • • H O o

u 04

o

o 1 z ID •H fH n 5? 10 •< w n «^ «< n •<

pq •H 04 U4

w

•<

\ • * a

a

M

a

(l4

a

a Oi < 04

a

\ « p

n

w

f

w

o

a • pq tno * 1 C/J

p

p • 1 cq

M

p

au pq • CQ pq a t p u •< fH t U4 w p K (M w w P R M w a hI w M M M M M (U w

«<

p 04

a

a

a

a

a

a

(0 > M o«
o o o

N 00 0» ro n D ro fo o o o o o

o

H

Z

o o

in «o

o o

H H

o

H

O

H

H (I) W W W

;3

» D

» *

p

IH

p MHO

M w z

a« H w

o

o
ooooooooooooooooooooo

r^<Mro<t^n^o^-cDo\o»^(^Jro^ln«or^ooc^o^ 000000000*^>-<«-i«H*ii«>-(>-«>>4<-4<-t ooooooooooooooooooo ooooooooooooooooooooo

o o o o

CM r) in

CM (M (\J (\J t\J (M

o o o o o o

o o o o

Q

n

ro

(M

n

ro

ro

\

n

I

I

CM *N

<

W H

Q H

U P

m V

< pen

n p M Ul

to

H

(I

00

UiK

t-t M

uH

at

Q

M

H

H

06

• Q

•<

069

D •

a« CM % :e • n

CM

• :«

« n

•< •<

u u

ro

in

ro ro

fi in

a« ^

m 06 W Q

a H

ID • * «^ H ro

• in

ro w

w m

Q 06 s a.

row W

cWa;

CM CM

VI V a; M << 8P\ \

V o

ai at

ai 2

u u

ro

•* w

CMn

ro

p

06

ro

•*w

n

OS\ro ^ to ro z

5» <

« «

in H

^ \

as V

CM a

* U

*

ro

>

-> «

to* ro • <t in w

^ Q

to «

(fa• ro ^ <t ro w

in H

u

06 <

< w

^ <

\ \

o

in

ro

•*•

--^ in ro ^ It w «

•I in o ^

WJ w \b

Ui V Nm

» to V

ai as © M

(J n

CM M

CM ai «

04 m 1^ » »^ w

• • b

CM r4 Si^

11 *

06*

M «

at ^w

w CO CO b

fa w

I *

w

w fa

to CQ fa « w ^

I M H

* w

*^ m

w as fa 06 n «

n

* I

I n

7

w

fa

• fa ^ p M Z

« M

»^ »^ H

w Z

<£> V fa V

Q p n u

o o

fa n

ro

H •

o o

CM CM V V

Q »

fa u

«• w <

9 K fa

fa I * fa

fa

06fa

•<fa

I<

fa fa

fa • M P

at « z z

M K

w H P

CO ^ Z H P

C^ O fa Z

• (0 U K W

o

CM
116

o o o

•-• CO 10 o o o

o o

o o

<M CVI <\J o o o

ui 0} tn

o

o

o o o o o

in vo N eo o o o o o

Q O O Q O

0 o o o o

(M N CJ (M CVJ u u u o u

01 (Q 0} 01 0]

o

o

o

o

o o o o o o

<t in o N

o

o

o o

o o

o o

o o

<M

10

in o

* \

n

o

CO

fVI tvj

u

n

n

0)

m H

o

o

o

o

uiaioivivivitnmm

o

ro

ooooooooooooooooo

ooo>o»-»(vjrO'*io«ot^ooa>o^(\jro^

c\i(vjwM<\JW(M<vj<vj(Mforororoio ooooooooooooooo OOQOOOOO o

o

o

o oooooo

C\J(MeM(VlM(M(\iCVICVI(\j(V)CViCM(VJCViCVJ(\J(VJ(VJ(ViC\J(\i(M(M(M

(L, w

Q

u

n

<

M

01

u

CMQ

I CQ

» Q

< n

fO

eo O

ID M

W h)

=> s

Ul M

M

-<

*^

fO «

Oi (VJ

• n

a< • ^ n

z * « n

Ul

-< <

\ \ z z

at

ro

-*w

(VJ n

• n

ro •

<^ in w

D w

M 01 ^ •

ro ro V • n in « 1^

ro uu « \

z

< w

P5\ \

z z

5 ^ o

ro U

•* *

<M O

a n

M CP

Ui o

« •

10 o

sr ro Q •< a oID U

M « u u

CU H

• •<

ro M

<f •

© •

• O ^ ro

ro w

w W

ro H

H <

98 K

V • n

• ro

<t ro •

<t o

ro

10 ^ < \

a •< u w \ \

z z

© ID u u

ro

in

Q

« to fO ^

• in

in

•-• w

w W

P \ VJ ©

N «

(\J GU a W

s \

z z

ai ai u u

in

*-t ^

• o

o ro

ro w

w n

ro D5

in o

.H ro o <

ro «

w

CM * U

<^

• O

<^ ro

in w

^ n

o •

row o

ro O w \ <

u S H H

M U

* <

w w

o w \ N

z z

ID ©

o

roH

U

-<W

ro o

w ro H H

u

n <; * w

tt o «

row o O ro h CM

< • H O

» ro \ w

z »^ o a

a m

V • o

M M

• M

M PQ M »

M ^

n w

« M

M M

OO • M

ro »-» M M

M •< M P6 >

00 w CM

XT w (fa <t • (fa

ro A M H

iU M * H • M

Pki CM < • CM w CM o O M w (fa ^ PS (fa Z Q w Z (U W Z

• o O • < m <^ O • H H • z CM • • in M CM •< -< CM

1 * I • < n • • w • 0! O * M M «^ M H I • a o r-l M CM M

M M • M M w P P

z o • • 00 z M M M M • • o M M

w t-« • M M M M m in w w ft m in a w M M «< n w w w >

p K? (fa -< m « w (fa (fa (fa w (fa Q Q (Q Ul P H M H M p (fa 06 p

o O fi

in St
11

ooooooooooooooooo

rorororon4^-*<t*4 ^ <t < tn u) ooooooooooooooooo

(M(Vi(VJ(M(\Jca(M(M(M<M(M(M(VI(\JCViN<M

MCQCQcacQCQoieacqcQcqcQcacacQcooi

oooooooooooooooooo

c\jro<ttf)<ONooavo^(\jro<»o«oh«ooo»

oooooooooooooooooo

(M<MC\J(VI(\J(MC\i(M(Ni(M(Vi(M(VI(MC\i(M(M(M

<^M H H

w w w

ID M

w

«

M

*H

S

H

S

*H

• « * n

H at >* •M w ID x>

(Q PQ m

fa* w

(Q to w CO M

fa H fa * fa * « •1 « -< « m

M M

H • M O H Q

* W M M • «> <^ M 9« M n M n *^ M b M M H

M w M w w

n cQ w n n n (0 > >

•-I fa w fa 3 fa fa w M

tn 41 w *

fa fa fa fa * M « M « M « H

M

M M M M H M M M w H M

* M • M *M « M (K w w

M M M M M M M M O PS >

M << M 't) M •< M fa fa u

M M 53 M W M Pi (•> w fa

H> M < ^ M M M M * •d «

•J w h) w h) w O w w w w O o o o

w fa w fa fa w fa > w w > CM Q Q n

fa« * * *

fa

J* fa 4i • • • • * c fa fa w (M < CVJ H * 1 • 1 1 M M M • M ID M M M M H H M M

M w M w CM M w CM CM • CM * • CM m • CM • « CM w K fa w M fa w K fa fa fa *» fa fa •-J fa W W fa M fa fa Z fa W <^ w w z w Z Z

H H • H • H H * H H • • P P * & P * P P •

CVJ -< CM CM (M CO (0 CM (0 CO CM CO (0 CM CO CO CM o • o R o R « M • H R R CO R R R 00 • R R R 00 R R R fa R • • M • • H • M R H D M M p M

<^ w M M M M w M M M M z M M V ro M M V M M in M M M ro CM • • V * • CM 00 * • M o

w w to w H in w H in H O O to H in in H in ro

OS M > P6 > > Z

V fa V V O w O fa fa D P V D P fa V D p 03 P P c>

fa Q fa » O fa 06 O n fa K M M Q CO (0 p CO CO o M p (0 CO O p (0 CO o p

<M ro oo in O o 0^ 00 r-l 0* 00 10 in in in in CM in in CM in t-«
ooooooooo

ooooooooo

CVICM(VJ(M(\i<MCM(M(M OOOOOOOOO (naiaitnuimtnaif/i

© to

BQ n w w

« « n b

<\J CM O ID H M

s

• *

^ w w w

w Q • D P

(fl n CO (Q

• • •

« ^ M «^ M M M 2 M « • • W

• ro p

*^ >^ w <^ H P

w P V P P W 2

« W P (Q m K M

o

ro ro
o o o o o o o o o o o o o o o o o o o o

T-i ro <t in CO o w ro <* in CO o

o o o o o o o o o .-1 (NJ o o o o o o o o o o o o o o o o o o o o

o o o o o o o o o o o o o o o o o o o o

o u u u u o o u o u o u o o o u o u u o

•J »-) •4 *-) •4 o U o u u o o y u o u u u

n

(M

ro

in

fO

m

ttm

<

Sis

2"

S

a

n

H

*>

a,

•<

Q

-<K

"J

ro

«

in

ro

in

w

m

ro ^ ft,

c (M b

n SIS .H

• « aiM W fti

« cq

><>

CO

ro^ ro

<t in •

t-» in w

0) w

ro ro

• *

in in

<

M M

H

P

«

n

p

01

M

(i

N \

L> 9S z

M ©

^ a a

0. a ;« :« V ^

M U y

\ \ 1-4 (fcl <

> ^J

qp o

row

ro«

ro •

\^^ ro

ro •

H

^

M I « M

»^ •

w

m

o

ro

IQ

<

M \ •

f-i • M ro M • N w CQ ^ I w PQ ^ '< AM ^ ^ ^ I roroo'^MMo-*

• •(Mw • wro*

O »^ 35

o

(Vi

O U4 M ro w H

A (5

O

n

o

H

• w "J

<! • M U4

• O « P

-« W 1^ Z Z

H • w M «

• B4 H P

*^ w w 25 M Q

w b b o w z

(L, M M U W U4 o

fO
120

o

o

o

o

o

ooo

o

o

o

o

oo

o

o

o

ooo

o

o

o

o

o

ro <* in

N 00

o •-4 (Vi rn in vO r- COo C\I ro in

o

ooo

o

o

o

o

o f-« 1-4 f-t •-t C\J

w CM CVi

o

o

o

o

o

ooo

o

o

o

o

o

oo

Ooo

o

o

o

o

ooo

o

ooo

o

o

oo

o

o

o

ooo

o

o

o

o

o

o

o

o

o

o

cu p. ou a« cu (t OU 0. BU Ou Pu ou QU z z z z z z z z z z z z z z z z z z z z z z z

CO (0 CO CO

w CO (0 CO CO CO CO CO CO CO CO CO CO CO CO CO CO CO CO CO CO

msnan

(M

roCM

10r4

U)N

roH

CO

CVi <

(It

U « z wQUm<nM

COHW

ro

ts)

Z I *

10 CM • a* n z I *

•<

Z

CD « * z

hJ

\w

H

M

M

\

o z H ©

»-)

M

a,

X

S3• ft

:« ai10 • K a c

-< •

ro

CVJw« • a ro m ^

t» ^ \ CM

CO < < (4 z z

ro

inKH

ro

in

ro

in

U4

ro

ou

OS

ro

<*• w><

ro*Q

ron«NCM

COz ro

-< ^

?5 in \ r« Z w a (MQ

(Mn• ro^

ro •» ^ in w •-• D w W CO

ro

fO

w <f n in «

ro tb

CO

\ z

•< « K < M >• \ \ z z a a a a

(0

uu

row ro

{H *• w M Z H U)

(Vi CO

a ^

10 •a M w• CO a b

CO ^ b • * w • w w I UU a M w• CO a Pu

CO i-t * ^ * w M « u

•J< z

t>

w ro

H vj WWz z V o o a a CM a a c o c UUP

EUZ

O• QO•Izw

CO

a *

pa

M • a

<^ (36 p o pq ©P b

ap

—

H ^ • i U H < w b P • a

•< cm

w

M

10

s

K

"> «

r11 M h-l • w(X4 (fa 1 1M• •

-Jw

*^

•H

M 1 u w

MM• 06

*^w

cu

MM

w

w z Z

H U) 1 CO

w (0

Uo • R •< -< z

lb

M

M • I w u z z z

-< (U

M

Ut CO

u

DZMHZ

o

o

CM

ZKOH

PZ

06

W
o o o o o o o o O

ro tn CO o o o o o o o o o

o o o o o o o o o

o o o o o o o o o

•J •J •J I-) •J »J o o o o o o o

SI N N M

M

H

04

*

ro

ro >

<<

(M a«

ro ro 06 (M •

1*) n

;i ro o • M

M

W

•-«* o

ro ro

ro '<t in w • M P

1 «

(Vl 10 • ro

{J ro M * ro <t M

<

w h w m

CM pq M n

00 • ^ K (D

M ^ \ 4

Hi \.(M \ •>« N < W < Oh M

U •J U Z w

U 0( •< W •< < V

* at; K m

w M H >> « s

Q H M w w o

P u 7C 1 • m O M V V O

« 10 M

cq w

m p »l ID p H M o o u U a W >*

* o

M in M

(I

o o O O o O o o o O

o (\J ro in oo

»-• •-I t-» •-t o o o o o o o O o o

o o o o o o o o o o

h) h4 o o u o o o O o o

N N N N N

'/

N

1M

p o

n \ in

•

1 X

M

•-4 Q 1 31 M

ID >^

4i M >4 M • M

C\i * <H

di <M 04 04 Z z z H

• z • « w

o ro t (M M ro >4 S • M • w N W

o H w • M H * • • z

o • ro o O M

o« ro M p

H z w > H P

w >^ o © M z

Q Q P W w

o o o

CM ro <t
122

oooooooooooooo

000000000«-ir4^ oooooooooooo oooooooooooo HHHHHHHHHHHH

O O

P P

oooooooooooooooo

in«osooo>o^cMrO'4-invOh-ooo>o

oooooooooooooooo

HHHHHHHHHHHHHHHH PPPPPPPPPPPPPppp

o o o ^ c\i n

fO ro ro m

o o o o

o

« n p

o o

o

o

roo

(\i

o

row

X

<

O

ro ro

H

P

A*P

CM

i

«

%

o

"<w

n

• »

w <

H w

P

(t

P

O

ro

M

CO

w

zM M

H H

P V

M

n

p

U2

> ^ m

ro w

5t in •< •O w (b

ro H • w P

> cu o

• « \0

n w

a* H M

« M

H

N

ro

O

row H J

N

u

H

ro

•

CMz

n

>*

n

n

M

K

m H

a

«»

a,n

<:

PS

n

z «

CO (/]

•< •< W

n

H

QH

«

ca <^ w

ro «

w • Q Q

;«O ID » u

* H

M • o in H w

;«ID

« K ^ 0* ro M

* • « w

a Q in « ai ^4 « ID W • Ul W M

z ro « M ^ U)

M SB •-• Z w \\M

W X

< N

O O H

W W M

a ;« z

V V M

z z z

V V C7

ai a

a X a

V ID ID u y w

K

(i] H 10 M <t

in

in ^w

• H «^ •

in

ro w «*• 0, • u in

in « ^ o

w M

w ro

in in

o

roO

n

o

row

S

n

o

ro (VI

0* u« u w

•

in ^

^ ro w< <

(t< w

w «

» m

n 06 ^ (\i z z ro < < KM*

H H U)

Z Z w

<D ID a

31 98 CM

ai a K

<D V a

u u «*

ro

row > a

CM <

ta tH• X

• ro Q

ro m

ro • -St W

< U)

w Q

p W • m H

ro

ro *

ro

in n in M

o

roH

O

<

a

u

«w

H

ro

<t

<t o

ro

Z

< «

K -<>

N \

z z

a a

a a

O ID

a. w

ro^ o

• CM inw <

© \ PS Wi

z z

V o

a a

o

CO -< \ \

z z

a a

a©

ro

M CM « P

• m

<-» it ro ^

^ o

« n

-<w

O

ro o

w fO < w

H H

W U

n <

• w ro

pa o « « ro o

w o ro C ro

pi w

(fa CM

PS m N -< • u

a. H

H O

PS ro w

S w \

z ^ z

(D P O

a V a

a m a O • ID U U

H

•<PS

-g

o »

ro • w a

m H

PS p

iJ ID • «

<-> H

O H

ro P

w ID

PS n

• p >-» o • ^

ro ro w 4 <^ w •-• ^

H< w ^

• H W ^ P P5 O •

ro

ro

vo <

CM w

- 0< H w l> U U (fa < z

U) Q

PSz m CO IDa ^

a •< ^ w

U PS

oH

P

cu

(D

H

(Q

H

(a

<

PS

(fa

S)M

CQ

(fa (q

H

01

*X

ro*

ro

Z

V

z

(fa ;a

Ul

l-l <

(fa *o

\

p*M

M

m

N

o

•

\ m

P w

o

H

P

'D

O

H

P

ID

a

PS

(fa

«0

M

<!H

CO

Oi

(4

H

Ul

(fa

X PS * w

X m

o a

* P

hI Z

w

m O "H ^ n

* r1H

-<a

PS

(fa

in

CMo

II (D

a

(u

CO

(fa <

X

ro (7>

in «

CM >^

^4 (I] Q H

m M

i* PS

» »
123

o

o

o 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

in (Ni fO in vO 00 0 «-« (VI ro ro ro ro <! •«t «t <t <t in in in o o o 0 0 0 0 0 0 0 0 0 0 0 o o o 0 0 0 0 0 0 0 0 0 0 0 0 H H H H H H H H H H H H H H H H H H P D 0 D D D P P P P P P P p P © TD T> « © Ci ID M• nHD OS «o ro (M * H tn • U

(M in ^

ft • H td a H l-J riw • •

in

s

H • ttJ 0

tu

a *

•<!

P in

W • M

M

# PQ •< w r> .

NN

w «0

M 1*1 (VI <VI • M * X* • • • H in in 0

r«

M

n 9 w

X rvj if) II # • to • • w W CM M • * > w M w 01 (0 P • vO H (L w H w w H

H

w (M ro 1 • w < M

W

w "<

W •<

to

H

H

H (VI >^

H

HH

H

M •A

M

M

H

M

Pi

DS

K CI

P » a

10 in (M in

(VJ (Vi 0 0 0 0 0 0 0 0 0 0 0 0 000 0 0

ro «!• in 00 •-• (VI ro in \o

N CO 0* in in in in in in in ^0 « <o <o «0 ^0 VO <0 0 0 0 0 0 0 0 0 0 000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ^» %^ 0 0 H H H H H H H H H H H H H H H H p p D p p p p p P P P P P ^ © © © ©a

It

ffl

(VI ID

III

M CD UlZ 1 «

at SB

X

09 (VI

(t) (VI

fVI

X )-)

M (VI « H It 

(VI

K 03 H ^ M

M CDQ *^

as (Vt 4! *^

X

(VI

w

M

M rM

CO CMM

M 0) 10 M 1 H (\|

n

w

to• in

(VI •

«-» (M ro

M til

a

Q (b

P

Q ^ - 111

(/] •t z « M :s X • H M « W W ^ ro cvi

W

W

H

P «l « H (0

m1*1

H

w Ul ^ ft a» 9 w n w P M n «» K n w H Oi H H ^ U •=11 • P (VI On M (Q (Vi Pt! M p H •» ID M ^# Q^1 a^ M BP in z P en • H M (VI Q (VI z K M ><! in am (VJ PC p < » ^9 M

CO (K CO 0^

H II (VI f* n »-« X «i X <^ (VI • • • Q> rH •» V It *4 • «t vO M yO w *0 vO w W 0 w w H w w H w w H w H 35 w w U t V W «<; W M <; W •< w H H H H H H (VI H S - H H M w • M M M * M P< • M w Oi X (b ID W 5? K (VJ W © X (VJ U4 M » M » (b (VJ H * 41 «

CO CO

(VI f-4 •-»
124

o

oO

O

o

o

O

O

OOOOOoo

o

ooo

o

o

oOOO

oOo

o

o

o

O

o

o

o

o cvj ro 1/) vO 00 0>

o •-• fO in ^ oo

o (Vi <t in <0 oo

o

n <t h- N N N 00 GO CO CO CO 00 00 CO oo 0» 0* c^ o o o o o o o O o O o O O O O O o o o o o o o o o o o o o o o o o o o o o O o o o O O O o o o o o o o o o o o o O o o o o o o o o o o o o H H H H H H H H H H H H H H H H H H H H H H H 8 H H H H H H H H H H H D D P D D D D D D P p P P P P p p p P P P P P P P P P P P P P P

10

©©© ID

© t> TO

wKP

CO

nowHP

-<Om

06Ui

X (fl ro W

in

M• b< oM

CO

p

Q

h)O

Mw

H - < • P ^ H M < • K => W * 04 Hp M a*

I*

M

H2HZupHP

cuHP©WH

inswHWWMnHH

HiH<wH

•<U)HWnH

B

-<Po

•<O • QU

COPOHw

CO

CO CO

CO S5 w © « M

ro

in iJ . • H - « H ^ z • Q M vO * * * H©

©oH

COP

CO

«0

UW

Pu

CO

(MnunH

U4©wH <MH

0.

-<

cvi HH

CO©o

©w

HP©

M

(I) (0

Pd

HHm

H

ill

H

K

U4

H

O©H©O CM•

CM fO (M I • pa

»^H .-I 0>

w at

(1.

n

M

Q »^

(Vi

HH

ro

a

a

M £3

m .COI * w • ^ « M

:«p

CM WJR © p

CVJ

CVJ z o

*Q

(VJ

O

M

R • • M R H O M ^ w tfJ © © R P (b H

O

(VJ

(Vi

o

c^ o» • wu CVJ

(fa

o •J Oi

(fa

© z p fO n fO (VJ fO >

<! ro < M

ro GS 4» CO ro M • Mw H (VJ w M H ©

(fa

(fa

H

M

M

CO

CO (fa z

(I] d) M (fa u

(4

o CO

eu

o

X (fa CO in (^ 4k • •

5iM

H >^

m * e

©

w CO «» H W X

w » (^ 41

n (fan 41

H

U < M

(U

H

ifa

K v3 in

© * UJ

% •

(fa > m> (VJ M © M fH o Cfa • (4 w Q

C^ CJ • in Q "<! X © <•

M Pti • 4k

X • « CO • M (VJ ai H H IB 4» w pq H 4» M w (0 « •-J w X B> H

ro 41 (VI 41 C^

H

W

© (VJ (VJ

X < (VJX

"<

R (fa (fa U4 • «» • (VJ

M 41 4t «t

n (U

R

R

P »0 EC] vO «0 \0

H Pi <^

Z

ww

HH

w

H >^w

•-•

m *^

M

M UJ (U -<

w •< « (4 (JU (U

Ui * w

w

M

H

XH

H

H

u

H

K

Z

M

M oi:

M St

H

MM

©

V •<P

©© (K

© (K

©© (K

K

Z

a (fa (fa (fau (fa (fa (fa

r4o ro t-4 (VJ

(VJ (VJ tvi

u

o

u
125

oooooooooooooooo

lO<o^*ooa»o•-•<M^o<ttn^o^-oo(^o

00000»^»^»-«»-i«-»«H«-4»-«<^»HtVI

o o o o o o o

(vj ro ^ If) \0 N

o o o o o o o o o o o o

o rj oo n ro ro ro fO ro ro ro w w o o

H H H H H H H H H H

P D P P P P P

"»»' TD ID 85 83 TO TO

«^ u

o

ft. >**

»

H

M

H

M

ft. »

U

M

<

H

5

STi

©

toH

w w

Si/ H H

«)

II H

W H \u <

HmWJ to M

W *

»—

P H

l/> «« V' M

> w

w

•< • ^

ro < o 0»

V m

it M • o UJ r-4 vO R • M w vO <o w Z H H » H w

W M EU w -< w w

H H • P <^ V ?9 H H

M H H M M M

W u Ul w

U M M M M

in

ro o O

o

H

P

O

H

P

ID

(\J

ro

(M H ^ <

H

N U,

SI M

* H

;« w

w p

O M

H

z z

p p

W 06

o

H

P

^)

O

u

oM

:«»

U

COp

w

«

H

(bt>

01M

H

-<K

Z

M

H

P

Q

«

W

0.W

o

ooooooooooooooooooooo HHHHHHHHHHHHHHHHHHHHH PPPPPPPPPPPPPPPPPPPPP

K •

»H

\0 w \0 w

M

-

U< M

*

(M

CO

R

m

o

z

o

<

n *H w CMz

\0 •-« (M • • R R \0 «0 M

w w

W M O >^ H H ro ro K W V O

» » P Q

w o

« • • R R ^ ^

(VJ w ro w

H H POP <D P <D

fO ro

a

*•»m

*i ro w H

H « P ^ C • w Q\

ro «» ro H

• ro w

vO w P

w M Z

W ^ M Han

M W Z

Oti ID tD ^ (b U

ro o

ro ro

01K

W

H

W

A*w

z

H

00

R

z

p

<

w

«

ro m

w H

S3

W ID

ro

(7« ^ ft r-l • m *

«0 ^0 w w *9 W W R H H

M M PL4 w » m

^ »

O4 R <w Z *^ « • U

CM P<

R R (Li M R

Www

(b

n

>^

R

<->

in \o >^ ro ro w

H

ID <>:> P

P P ID

vOro

ro

IDP

O

•

R<^

»^ W M

H

P

ID

ro
126

ooooooooooooooooooooooooooooooooooo ©•-•<viro^i/)*ONoo(7>o^(\iro^invONoo(7>o^(\jro^in«ONooo>o.^c\jr)-*

ooooooooooooooooooooooooooooooooooo HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH

H

(fa

•JH

u

COp

H

.-^

(0

!;M

M

>

H

Qt

K

<

m

H

M

H

H H

(0 CO M H

W

H H

(fa TOH

(0

M

« "<

<D (^

(fa M

• I M M

o o

W W

H H

< <

PC «

P4 >

(fa W

Ifa K

• • w ^ >0 vO P

M M ^ w Z

(X >

V w

[fa K

Ifa 0«

roo

M M

u

o

H

2SP

M

IH

Z

P

o

o

H

a

H

P

©

10*

H

IDO

^ p

H H

o u

o

o

(fa

tn

en<

H

M

tD

H

U

(4

H

•<

H

H

5j

(4

K

n

(fa

w

H

< M ro K WW

(fa H

CQ

5S

H ro M CV) R M

M •

ro

w n

M a w

w « H

ft

<-» *

00 •-' ro n ^ w

w H

w

» (fa

GOro

^ ^

«0 vO w

w w H

(U W <

M H a

M M OS

06 K V

» » (fa

(U

H

• M

IM

w «

(fa «>H

01

<^ ^ P5 M M <

w w M

00

<»

H

ro M

06

*M

M

m

w

o

o o

• • w w

(fa (fa

(K

M

M

M

m

w

(fa

*M

M

M ^ .J

(fa

M

C • U M

w ^ W — H

M <

w PC a

n w

w «

(fa V

H (fa

H

* p

cu

M M

(4 M

w

H

P

(faI

• M

M M

M M

M PQ 2 PC

h) w

w H

(fa P

« (fa M

• M

-< w

p

RU

<!

[fa

ID

(fa

PS

(fa PC (Q (fa
127

oooooooooooo

ln^o^-ooo>o^<vJn«^•u)^o Nf^h>KNQO0OCO0OOOCO0O

oooooooooooo HHHHHHHHHHHH DPPPDDPDDDPC}

ooooooooooo

f^oo(^o.-i(Mn'*tf)«OhOOOOOOOOOOO HHHHHHHHHHH PPPPPPPPPPP

OOOOOOOOOOOO

(oo>o.-«wro<i-in\Of^ooo> 0«Q^OOOOOOOOOO

oooooooooooo

PPPPPPPPPPPP

M M M

w w

(VI (Vi (M

°H H »

w K K

(XJ * « * H

K H M M

* « • <^

PS •< M n H

M < 5h w >* w

m O T) Q ID K M M «5 H

w n Pi (X) OS © « *^ « w w * w w 1^V (0 M

* lb • (b • * Mm • « < * M • M M

M M M M Q M o M

w M M M •© M M

w

LBK rv

LRB (LA MM

M (

JB

LBK M

M

M (

JB

LBK n (0 m 0} w n w

Hi (b H lb lb H (b H

w P Wl * (b w « P w « P

(b ftlb • (b (b (b (b

• M * M 1 M f • • M «

«^ ^

H M M M M M M H M M

• M M « M • M M « M M

M M M M a m M M M M M M

M •< M K M -< M •< M •4 M w »4 M QH M PS

•< W •0 i-i w *^ -< •-I w •< w

H (M »-l M w H

<^ M H M h) H

M w « w P H w w (b P H w lb P M w w lb P M

w > * w >• (b (b l> Cb w (b « (b w

w w X> 1 K w M 1 PS Ul 1 w H M K H M M ^ €> w M (L M <» lb » M • M lb M m M lb * • M M (b * * M M lb M 14 lb « • M M lb H * O M w M M M M w M M * w M M M « u M M M M 41 w u "< K <5 w K w < W K w <! OS << w U

> •< »< M W « > < (b M W K W > < lb M W > •< lb w w OS OS > "<

(1) H H •J w OT H H h) M lb H H © lb 01 H h) lb

f M -< i

U) 1 •< w lb W (Q 1 «< w lb K Ul 1 '< w (b OS CQ

* • 1 o « f I 1 1 o » OS • R N • o » n ft H o PS • II R • p <^ • • p « P

#^ -* N • p

M M M n Oi M M M PU M M M M <^ (b M M M M (b M M

M M H w • • M M w w • • M M • w • M M « w • > M M H n u K >• M M w K H w w > H W U OS >

w > <; w w « > < ^ > m > •< w

M (b w u. (b PS > V Ul lb lb PC > O lb lb > © U lb lb OS >

n M lb n lb « C/i lb m m lb lb w n n O lb PS (n lb OS n

(VI 00
128

ooooooooooooooooo

o^^l^o•<*ln»t)soo(^o»-•(\J^o*ln^o

W(M(M<>i<MCMtM(\JCVJtVJW(\J(\J<M(\JC\JtM OOOOOOOOOOOOOOOOO HHHHHHHHHHHHHHHHH

oooooooooooooooooo

r^ooo>o»^<Mro<t'«*Oh-ooo>0'<<\JrO'*

oooooooooooooooooo HHHHHHHHHHHHHHHHHH

7iM

H

U

«<U

«

U

"<(4

«

10

0)

M

H

«<U

K

U

K

>

m

06

<^ Q

(0 -< M

H H H Q

CO CO •< W

M M

« • H w 06 w

< •< <

H H w

01 W

M M u

N « H H

H M tn • in ID Ui o < H

M H (VJ u •< w • n «

H K CO > H

w w M M H

> o cm X» w w p P

» X 1? « w < < M

OS * M • M

« M M H -< <

M M w • "< (M • • M o

< K tb « • K > M P m M Z M M in •

ID • • O o « • w t w w o in

1 U. o H o • • (NJ > t (VI n <o n

w U< cx W • * w M z (VJ W a N I M ro M a M w M <t o W «» i

1 P w p 1 p I P M M M « oo M P p Z w (VJ p ^0 w

M z w H H H Z M M M M M M 1 Z z M w H z H

* M M M M w M u M in w <0 w w W w <^ M M u* w < V M W •< H • • H t-« p M < H oo ?<! oo >• H H M S

M M M M Z w Z •< < < w < z z w M 06 z M «

w K ft o u VJ V tK ^) V © 06 n > U U M M u n X >* Q M M M u u M b O u

O O ro in (VJ ro <0 in in in

in ro «* ^ 00 00 o o

(\l CVI
129

oooooooooooo

<\l(M(VJ(\|(M04(Vi(VJ(>i(M<MCVI oooooooooooo HHHHHHHHHHHH

OOOOOOOOOOOO r^ooo^o^(vjro^in<ONoo

W(\JtM(M(VJ<\J<M(VI(\J(\Jt\i<VI OOOOOOOOOOOO HHHHHHHHHHHH

OOOOOOOOOOO o»o^wrO'*in»Oh-ooo»

OOOOOOOOOOO

DPDPPDDPPDP

P

M

o

w

asP

PS

s s

^R

Cfl ^ (VI OS • > r

VH

• ¥m I

CVJ M ^ M

~ W M

M M >^

1

-

3 SI

in M • w

w 00 « H X K

H

• H H

2; • • (M ©

qj ^ CO <* M M w M ©

H w • M w H

M H M m

V] M (f) >< H •< * 06 •<} X 8 H

iW H M N H • M N M • • <» w w * H

X QO n w H O

• •- • • t-« •< f o a • so w w M < o M 0. H w H H X H • M

w <; w w •< OO M « « y» M

M P W N

M M W M • V < w w

«:> M « W M H

» bu N M M N Q M M

«» CO n

CO oo

H

MHO

n

-< • H

w

• m X

<

CO H « W

•t ro 00 m <

H -< • ©

X H

©

* • ©

H (VJ Z M • M P • Ui H © X p M

H w M H M « H • CM • O

• M P M R w

1 © «^ • > CQ 00 X n H

M U4 • w OO V < w M • (M <! Z ro u

P m H CO b H p R <£) w H P CO p

N Z <! m M n M H w H M © z

• • M H fO R © r) R M W «< u V M

n n H M 1 n N n H P H :m M H

•< •< ^ n < CQ © H 06 M w Z

H H © < H © <! H © W © M fti © ©

N M M H H z; U M M

O 00 r4 CM

00 CO OO 00 00
o o o O O o o o o O o o O O O o o o o o o O o

Q CM ro tn to oo o (VI ro «t in oo (^ o (M

00 oo 00 00 OO oo oo 00 00 oo 0> 0» (7» (^ Ov Ov o o o

(M (M CVJ w w CM CM CM CM CVI (VI (VI (M CVI CM CM CM (M r> ro ro o O O o o o o O o o o o o o O o o o O o o o o

H H H H H H H H H H H H H H H H H H H H H H H

Q CJ D D D

ID tD o 10 TO l> ID

n

1^V

s

CL. ft. r \ MSs mi*< rn

tm

fi.

rAi/i

fii a< rt) IDH

V tm (Q

1 1\ M it Ul

"V cn H

4k m1*1 5S ID ID M H W w CO H

Vi/ mMH Ik m

CM M

• S Ifl H w 41 41 *^ lAUl * H MW ¥S •H /VI *H M w (VI w n

.

I*< MM w 41 • • M o

hm ^ • w D

W < o (0 (0 •

H • (VI pq • (M (VI X o (^

w WJ n • • M 1^ Ul M CM o

n W vO .r\ t5M • & * (M H H 1 w w CM M 35 H V

O o M M Cb < w W O w M M X V ft! • CVJ tM • • P < V M H (VI H D M < H

M M V H H M M M M M ft! ;z; ID w H o

(b <D w u W M K ft! ft! ft: TD U (b (/J V W M Q O ?^ M M M I» Q U M M M O Ui w

O 0> O (VJ 0> o o o in O

(VI (VI (VI (VI (VI
131

1) Patfjnkar, S. V. and Spaldin^^, D. li., Hoa t a_nd J4a, r3G , Tranr^fer in

Bourida rx J^J_QS_s_ (Intertext Books, london, 197C).

2) Lovachev, L. A. and feganova, Z. I., Dokl. Akad. Nauk S3SR 138,

1087 (1969).

3) Spalding, D. B. and Stephenson, P. L. , Proc. Roy. Soc. Lond. A3 24.

315 (1971).

4.) Hirschfelder, J. 0., Curtiss, C. F. , and Bird, R. B., Molecular

Theory of Gase s and Liq uids (Wiley, New York, 195-4).

5) lAndsay, A. L. and Broroley, L. A., Ind. Eng. Chem. ^2, 1508 (1950).

Al

APPENDIX A.

8A

DERIVATION OF TIIE ONE-DIMENSIONAL LAMINAK FLAME EQUATIONS

Consider a control volume 6V

oriented as shown in Fig. Al . Since

the flame is one -dimensional, there

will be gradients only in the y-direction. The gas mixture flows from

right to left through the + face of

the control volume and out through

the - face. The area of each

face is 6 A, and the length of the

volume is 6y. The velocity of the

mixture v is the mass average flow

velocity defined by the equation

by

y axis

Figure Al

.

Orientation of control volume

V = p .CM. v.

'J J J J where v. is the average velocity,

J C, the concentration, and M. the

molecular weight of each chemical species; p is the density of the mixture.

There are three conservation equations to be derived;

a) the conservation of total mass;

b) the conservation of each chemical species;

c) the conservation of energy.

a) Conservation of total mass

Let p^ be the density of the gas mixture entering the + face of

6V, and let p be the density of the gas leaving at the - face. Let

v_|_ and v_ bo the mass average velocities of the gas at the + and -

faces. If p is in units of kg/m and v in m/s, then the number of kg

of gas entering 6V per second is 6Av_|JD^, and the amount leaving is
A2

6Av p . The difference between these two quantities gives the rate at

which the mass of gas in 6V changes. This is 6p6V/6t, where 6p is the

change Ln the average density of the gas in 6V in the time interval 6t.

Thus we get the equation,

6p6V/6t = - 6A(v^^ - v^p_).

Note that v^ and v are negative numbers because the flow is from right

to left. Dividing this equation by 6V, we obtain the expression

6p/6t = - (v_|_p_^ - v_p_)/6y, since 6V - 6A 6y. If 6y is sufficiently

small, this equation can be written in derivative notation;

ap/at = - b{pv)/^y. (A1)

This is the mass conservation equation.

b) Species conservation equations

Let Y. be the mass fraction of the chemical species j. Then I.p

is the density of species j in the mixture. We can derive the same

type of equation for the rate at which the mass of j changes in 6V as

we did for the change of total mass. We get,

6(pY.)6V/6t = - 6A [^(v + V.)^(pY.)^ - (v + V.)_(pY.)_j + R.6V.

Here, V. is the diffusion velocity of species j, and R. is the mass

production rate of species j arising from chemical reactions in 6V.

Dividing this equation by 6V and passing over to derivative notation,

we get the species consei-va tion equations,

a(pY.)/^t = [(v ^ \^j)(pYj)] + ftj. j-1,n-1 . (A2)

There will be one of these equations for each species except one.

The mass fraction of the n'th species is related to the mass fractions

of the remaining species by the requirement that the sum of all the

mass fractions equal unity; i.e. . - 1.
A3

c) Energy conservation equation

Let h. be the enthalpy per unit mass of the pure species j at the

J temperature T. The enthalpy of the mixture is then h = .Y .h .

.

Let A. be the thermal conductivity of the mixture (J m s °K~'').

The rate at which energy enters the + face of the control volume

due to heat conduction arising from a temperature gradient ( ^T'/^y)^

at the + face is (A. ST/ ^y) ^6A. In addition, energy is carried into

6V through the + face by specios j. This occurs at a rate

(v + V. ) pY .h . ) _^6A. Equating the rate of change in the total energy

in 6V to the difference between the amount entering and leaving by

heat conduction and the amounts carried in and out by each species,

we get,

6(ph)6V/6t = 6A [(?.aT/ay)^ - (^^T/^y)_]

Dividing this eqiiiition by 6V and going to derivative notation, we

obtain the energy conservation equation,

a(ph)/ ^t - ^T/ ay) - ^( Ej-lv >- Vj)^jPhj)- (A3)
A4

APPENDIX B.

DETERMINATION OF THE VELOCITY OF A STEADILY PROPAGATING FLAME

The equations derived in Appendix A give the time dependence per

unit volume of the tot'jl mass, the masses of each chemical species,

and the total energy in the system. For a steadily propagating flame,

the time derivatives will be zero in a coordinate system which moves

with the velocity of the flame front. Let us consider the consequences

of allowing the time derivatives in E(|3. A1 , A2, and A3 to vanish.

From A1 we get B(pv)/^y = 0; thus pv - M is constant throughout the

flame. M is the mass flow rate in the steady flame (kg m s ).

From A2, we get,

A (v + V.)(pY.) - R..

'

Integrating this equation from the hot to the cold side of the flame

gives,

a, (v + V.)(pY.) dy =X = "j^C^'^Pc - '-^ * ^j'h'p^o'hIf the integration extends sufficiently far intothe hot and cold regions,

the diffusion velocities (V.)„ and (V.) become negligible compared

to Vj-, and v^^ since the concentration gradients become small. Thus, we

can write this equation as,

^ R.ay = (pvy.)g - (ptY.)„ ^ M(Y.(, - Y ,

since (pv)q ~ ^P'^^H ~ ^ from the steady-state mass conservation equation

A1 . Any of the integrals

J'^^^^

^® used to calculate M or v^-,, which

is the steady flame velocity referred to the cold side of the flame.

Thus, the e(iuation for v^ is.

The solution of the flame ecjuations in the steady-state limit yields

the mass fractions Y. of each chemical species as a function of y, the

distance through the flame. From these values the mass production rates
A5

R can be calculated as a function of y. If the flame equations have

j r

been solved correctly, each of the integrals / j-^^" should yield the

same value for v^-,.

Consider now an adiabatic steadily propagating flame. In this

situation there is no heat conducted through the cold or hot boundaries.

Thus the temperature gradients vanish at the boundaries; i.e.,

( ^T/By)^ =^

( ^T/ 3y)jj = 0. From this, together with the condition

= (^•)u - O5 g^t from the steady-state form of the energy

equation A3, the relation.

Since pv is a constant this becomes,

= or simply, h^. - hjj.

Therefore, in an adiabatic flame, the enthalpy of the gas mixture is

the same at the hot and cold boundaries. Inside the flame, however,

the enthalpy will generally be different from h^ because ^T/ ^y and

the diffusion velocities V. are not zero there. There exists, however,

J a condition in which h remains constant throughout the flame. If the

energy flux from heat conductic^n \ BT/ ^y equals the energy flux from

diffusion p/^.V.Y.h.. at each point in the flame, then we see from A3

that the quantity pvh is constant throughout the flame. Since we know

from Al that pv is invarient, then h must also be invarient. This

condition is equivalent to assuming that the Lewis number DpC^/X for

the mixture is equal to unity at all points in the flame. The reason

for this equivalence will be discussed in Appendix G. Although this

condition is only approximately satisfied by real flames, it is often

envoked to avoid having to solve the energy equation.
APPENDIX C.

CALCULATION OF DIFFUSION VEIjOCITIES

To solve the species conservation equations it is necessary to

assume that Fick' s law adequately represents the mass diffusion processes in the flame. This law assumes that a particular diffusion

velocity V. is proportional to the concentration gradient of species

j only. If we consider each species in the mixture in turn to be one

component of a binary mixture, with all the other species lumped together as the other component, then it is possible to derive an expression for the proportionality factor in Fick' s law from the kinetic

theory of gases. This factor will be a function of all the possible

binary diffusion coefficients in the mixture and the concentrations of

all the species.

According to Fick' s law, the diffusion velocity of species j is

given by the eqiiation,

V = - A.x ^x /^y, (A5)

tJ J J (J where A . is an effective diffusion coefficient for species j in the

J mixture, and x^. is its mole fraction. Our task is to find an expression

for A. in terms of the mole fractions x. of all the species and the

J J binary diffusion coefficients T)^^ for all possible pairs of species in

the mixture.

The kinetic theory of gases'* gives the following expression for

c^x./^y in terms of the binary diffusion coefficients and the diffusion

J velocities of the different species y,

^x. Z3 EX. V.

§y" -'j^ - J

These equations are subject to the condition,

y),x.V. = 0.

E. X. - V.

i^j D.J

j^l,n (A6)

(A7)
A7

We now assume that all of the diffusion velocities except V. have the

J same value V^, Then A6 becomes,

^x./^y = x.(V^^ - V.) / ,x./D...

From A7 we get the relation,

-1

Substituting this expression for into A8 gives,

V^^ = -x.V.d - X.)

Solving for V^. , we get

(1 - X ) 1 ax

V. = - — ^ — . J X, x^. dy

J

(AS)

k

The proportionality factor A. in Pick's law is thus,

(1 - X.)

A, =,^^. (A9)
I A8

APPKNDJX D.

CALCULATION OF THE THERMAL CONDUCTIVITY OF THE FLAME

An empirical equation was used to calculate the thermal conductivity of the mixture of gases in the flame. This was developed by Lindsay and Bromley. It reqiires a knowledge of the pure component conductivities, heat capacitie,5, boiling points, and molecular weights.

The expression for A. is,

X = .X.

1 1

X.

1

.A. .X.

1 ij J

-1

(AID)

where X^ is the thermal conductivity of species i, and the quantities

A^^ are given by the formula,

A. . = i 1 +

2

(All)

The viscosity ratio ^^M^ is given by the expression.

= {X^/X.){C^. + 1 . 25R/M^. )

( Cp . + 1.25R/M.)

-1

where C . and C . are the heat capacities, and M. and M. the molecular

pi PJ 1 J weights of species i and j. R is the gas constant. The quantities

and S. are the Sutherland constants of i and j. For hydrogen, the value

J is 79° K. For other pure gases it is taken to be 1'5Tg where Tg is the

boiling point at 1 atmosphere pressure. The constant S^^. is the geometric mean of and except when one gas is strongly polar; then

S. . = 0.733(S.S.)'2,
A9

APPENDIX E.

VON MISES TRANSFORMATION OF THE FLAME EQUATIONS

Preliminaries

Before applying the von Hlses transformation to the flame equations

let Us insert Fick' s diffusion formula into the species conservation

equation and eliminate temperature from the energy equation.

Substitution of A5, the expression for V., into the species conJ servation equations A2 yields the equation,

bY. ^ bY

P J + pv

6t by hy

Using the identity,

/ Y.6x.\

\^ X by /

+ R,

Y. bx.

x.by

bY.

b^ <M>

b<M>

by"

bY.

+

where (M> = ^ ^ .x .M. is the average molecular weight of the mixture

we obtain.

Consider next the energy equation (A3). The enthalpy h is a function

of the mass fractions Y. and the temperature T; i.e., h = f(Y.,T).

The partial derivative of h with respect to y is.

(AI2)

hh ^ bY. ^dli^ hi,,

by L^i bT-

(A13)

where we used the relations bh/ bY . -- h. and bh/bl - C , the heat

J J P'

capacity of the mixture at consi:ant pressure. Solving AI3 for bT/ ^y,

and substituting this as well ati A5 into the energy conservation eqation A3, we get after some algebraic manipulation,

bh , bh _ > / X bh\ ^ i ^ ^"dT - bT^-^p^Tj "b

hA . Wp^ A^p 6(Y.<M» _ ^ 6Y.

by ^dy

(AU)
AlO

We now change from mass fractions to a different concentration

variable cp. defined as

cp. = Y./M.. (A15)

This has Units of moles/kg. The species and energy conservation

equations then take the form,

P^t " pv hy \ 0 oy / dy y

r / m.

{A17)

In A17 the quantities h"^ = h.M. are the enthalpies of the pure species

in J/mole.

Von Mises Transformation

We introduce a new spatial variable ^ in place of y. It is defined by the equations,

Wby ^ P, and bi/'/dt = -pv. (A18)

Applying this transformation to Al6 and A17 eliminates the convection terms pvbcp./by and pv6h/by. The mass conservation equation Al is

satisfied automatically because of the relations,

dp/dt h^mtby b^4^/hydt -b(pv)/dy.

We are left vrLth only the species and energy equations in the form,

dh

bt

\

J b^y b\^\j <M> b^

y
All

APPENDIX F.

SPALDING' S TRANSFORMATION OF THE FLAME EQUATIONS

The technique developed by Patankar and Spalding for solving the

two-dimensional boundary-layer equations can be applied to the solution

of the flame equations. To apply this method, a new spatial variable ou

Is Introduced. This is defined by the equation,

The quantities ^q^^h a^© "the values of ^at the cold and hot boundaries

of the flame. They are functions of time only. The advantage of using

0) is that if the calculation procedure is limited to the range 0 ^ (jo ^ 1

,

then it is automatically limited to that part of space where the important changes in the dependent variables cp. and h occur.

J After the transformation, the species and energy equations A19 and

A20 become.

^t

/a + buX p2 cp

\ (M) / dM,

(A22)

^h /a + baXbh ^ 1 _b /A^ ^\

L <M> S

bo) bu J)

(A23)

where ^1 ~ ~ ^j^* -^^ making the transformation we have used the relation for ^o/^t in the ^,t coordinate system, ^o/^t = (a + bcjo)/T].

The quantities a and b are defined by,

a = - d</^/dt, b = - di^Q/dt + d\l^/dt.

The derivatives d^^/dt and di/^/dt are the mass flow rates across the cold

and hot boundaries. In the steady-state flame they have the same constant
A12

value. When this limit is rea ched ^ cp ./^t -» 0, b 0, and A22, inteJ grated from o) = 0 to 1 , becomes,

This relation arises because a = PyV^j PqVq in the steady-state limit

(see A18); also the gradients of cp. and (M> are assumed to vanish at

the hot and cold boundaries. This gives an equation for v^ in terms

of R. expressed in the u), t coordinate system which is analogous to

t] M derived for R. in the y, t system. For v^, we have,

(A2^)

(A25)
A13

APPENDIX G.

UNITY LEWIS NUMBER APPROXIMATION

In Appendix B we demonstrated that in the adiabatic flame, if the

energy flux from heat conduction equals that from species diffusion at

each point in the flame, then the enthalpy remains constant throughout

the flame. We stated that this assumption was equivalent to the assumption that the Lewis number was unity at all points. We shall now demonstrate how this equivalence arises.

The Lewis number is defined for all species pairs as,

Le. . H pC D. .A, ij P iJ

where D. . is the binary diffusion coefficient for species pair i,j and

is the average specific heat at constant pressure for the mixture.

If we assume that Le = 1 for all species pairs, then we must assume

that all the binary diffusion coefficients D. . are equal. Suppose

that . = D. When this is substituted into A9, the expression for A-?

we get,

A. = 7-^=^-*^ = D, since 1 - x = •

Thus, the assumption Le 1 giver.,

pCpDA pCpA.A = 1; or Dp2 = A.p^ = Xp/C^.

Consider the energy equiition A23 in the o), t coordinate system under

steady-state conditions. Here, 6h/^t 0, b = 0, and h is a function

of uo only. We mugt make one further assumption, namely that the change

in the average molecular weight <M) is negligible throughout the flame.

Then the last term in A23 vanishes when A.p - A.p/C and in the steadyJ ^ p

state A23 becomes,
A14

Next, insert into A26 the variable U = dh/do). We get,

71 da)^Cpj Tl Cp doj

•

Solving for dU/doo we obtain.

Integrating this expression starting from the cold side of the flame

we get.

ln(U/l^) = (A27) J

f(u))du) , or U = Ibexpjy' f(u))dcu

where is the value of dh/dcu at the cold boundary. For the adiabatic flame tfe = 0 at the cold side; from A27 we see that if lb = 0,

then U = dh/dcD is zero for all values of u). Thus h remains unchanged

throughout the flame.

Therefore, if the enthalpy is to remain constant in the adiabatic

flame we must assume that A.p2 = Xp/C and that d(M) = Oj it is not

necessary to assume that Xp/C remains constant.

P
A15

APPENDIX H.

FINITE-DIFFERENCE FORMS OF THE FLAME EQUATIONS

The species and energy conservation equations A22 and A23 have

the same general form;

b4) + (a + 3(1))^^ ^ ^-efm ^ ^^^^^ ^t bu) \ "b^ /

where a = a/T|, 3 = b/T], y is the transport coefficient,

source term arising from the diffusion of average molecular weight,

and g^ls a soiirce term arising from chemical reaction in the species

equation, and from the difference between conductive and diffusional

transport of enthalpy in the energy equation. y> and

shown explicitly below;

are

m'

Energy Species

Tl^ Cp Tp /

J h hn A .p2 ^ b(M>\ 5r' A -p^-^ "^M^

1 „ b , /. 2 . W.. 1 R.. (

J

To set up a finite-difference equation for A28, we divide the uu axis

between 0 and 1 into N strips which need not have the same vridths.

The profile of the dependent variable <^ Is assumed to be linear between grid points. We consider control volumes whose boundaries lie

midway between the grid points. A portion of the grid is shovm in

Fig, A2. This shows the control volume and illustrates the linearprofile assumption. Equation A28 is now integrated over the control

volume shown iii Fig. A2. This leads to the following equation;
A16

The subscript P on<I>is used to denote its value before the time step 6t

is taken. An unscripted <I)refers to its value after the time step.

It can be shown after some surprisingly lengthly algebra that the

integral over the control volume of any function like which is

assumed to be linear between grid points has the value,

where Q = u>.i+1 - U)i-v

1-1 1+1

n+ - - to.

Note that a)_^ = '^'^^±+^

and U) = i-CuD + oj. ^

) . We introduce

1 1-1

some additional symbols

P = Tl/6t

G = -TIP = b

s + 3uu^) = a - ao^.G L_ = a - u)_G (A31

)

The derivatives of<t)at the control volume boundaries are approximated by

- , and (b^ = ^ (A32)

Figure A2.

+

Next, we define T_j_ and T_ by the equations

ll^Y = <I),.^ -d)), and Tl ^ = T_( cj) -<!>__). (A33)
A17

Thus vre have

=

TlY/(^i+-l -

'"i)

and T_ = TlY/(cOj_ - u)._^). (A3^)

We further define

= 11 jjd^"^ and S = Tl jQdto. (A35)

For the species chemical source term ;eJ^ only, we take S to have the form

S(species) = Sp + Sp , where ^ is the value after the time step fit

is taken. S ^ and Sj, will be functions of all the <J> for the different

species, but the values will be taken as the prestep ones. For all of

the other source terms, only the prestep values of the variables are

Used. Thus for these, S^ = 0.

Let Us then generalize Sp somewhat, and consider it to be compounded

as follows;

Sp = S^( species) + S^(energy) + S(energy) + Sp( species)

= Sp( species) (A36)

Thus,

S^ + S = Sp + Sp(|> , (A37)

where Sp and Sp are defined by A36.

We now multiply A29 by T] and insert into it the quantities defined

by A30, A31 , A34., A35, and A37. The resulting equation is linear in

the poststep values of cj) ; it can be solved for <I) at a particular grid

point in terms of the values of $ at the grid points on each side of it.

The equation one obtains has the form

* = A <I)^^ + B ct>__ + C, (A38)

The coefficients A, B, and C contain only prestep values of cj)

, and

are given by
A18

A - A'/D

B ~ B'/D

C = C'/D

A. H 2T^ - - i(P + G)J2^ (^^^

B' = 2T_ + L_ - i(P + Cr)fi_

C = iP(3*pfi +<J^p4.A ^ ^p_.J^_) ^ 2Sp

D = 2(T^ + T_) + - L_ + f(P + G)fi - 23^

= A' + B' + Pfl - 2Sp

According to Spalding, these finite difference coefficients have

a defect. This can appear if -g" L I , (the convective flux terms), become

greater than (the diffusive flux terms). It is not clear to me at

the present time why this problem exists. Spalding does discuss the problem at length, so the reader can look there for an explanation. The

defect is remedied by replacing T_^ in the formulas MO with the quantity

TJ defined by,

TJ = ^[Ti + I^LJ + - I^I^ill]. (Ml)

This has the following consequences;

i) If I^L^I € T^, then 2TJ + - 2T^ + L^.

Thus T^ = T, and the replacement of T by T« has no effect.

ii) If ^L^ > T_^, then 2T^| - o and 2Tf; f L_ 2L_.

iii) If iL < -T, , then 2T^|: - L -21^ and 2T^^ + L =0.

This modification of the formulas MO has the effect of neglecting the

diffusion term in the flame equations whenever it becomes somewhat

smaller than the convection term.

Let Us now see how the difference equations A38 are solved.
A19

Consider the equation for grid point 1;

^1 1 1+1 1 1-1 1 \ /

We consider N+1 grid points. (Remember that N is the number of strips

into which the au axis is divided. ) The hot boundary is arbitrarily

set a grid point i -- 2; thus the cold boundary will fall at i = N+2.

The grid points 1=1 and N+3 will not be used in the present calculations?

Equations A^2 can be written in the form

\ = A| (D.^^ + B|, i 3, N+1 (A^3)

where

A|.A3 .

.

Af = A. /(I - B.At^ J 1 r 1 1-1

B^ = (BiB|_^ + C.)/(1 - B.A^._^).

After calculating A| and B| for i = 3, N+1, one can easily calculate the

<t>^ from M3 starting from '^^1+^ ~ -^1+1 ^N+2 ^

^N+1 working down

to Note that the values of <t>.^ and ^-^^2 determined

by the boundary conditions as discussed in Appendix K.

@

Spalding unes the grid points i=1 and i=N+3 in the boundary-layer

problem in a special way to cope with large gradients in (J) at the

boundaries. Because such gradients do not arise in the flame equations,

this particular treatment is unnecessary. Although we could have

started the grid at i=1 , thi;; fact was not realized until later and

so the computer program was written with the hot boundary starting at

1=2.
A20

APPENDIX I.

CALCULATION OF THE SOURCE TERMS

There are four source terms to be calculated; for energy and

species arises from the diffusion of average molecular weight; (these

are normally quite small and can usually be neglected; ) and -ejffor

energy and species, which arises from the difference between diffusion

and conduction in the enert:^'" case, and from chemical reaction in the

species case. By far the most important source term is that for the

species. As described in Appendix H, this terra is given special treatment. Because this term can undergo very large changes during a time

step it is necessary to estimate at the beginning of a step its value

at the end of the time step. The chemical source term for a particular

species j is a function of the concentrations of the various species and

the temperature; = J^^{^>^>^>^, jT). Thus we could approximate a change in by a linear expression,

6jzd. <^ ^jda Sep. + 6cp, + + ^jda 6T.

Thus {jdPf ^ (^j-)p +

Instead of using (^j)p, the prestep value of y in the difference

equation, we want to use (>ei£j)pj 'the poststep value. We should, in

principle, use the above equation to calculate 6>€jf. • This, however,

involves the poststep concentrations of the other species as well as

j, and would miike the difference 0(iu/itionr; unpleasantly complex.

What we do then is ne^iiJect all of tho terms but the first, and write,

Spalding has found this approximation to be satisfactory. In actiaality,

we need only a rough estimate, because when the steady-state condition

is reached ^g/. ceases to change with time. In fact, it would probably

be possible to solve the difference ecjuations using only the prestep

value for j^^., but then much smaller time steps would be needed.
A21

This treatment is unnecessary for the other source terms because they

change little with time.

We shall first show how the chemical source term is calciilated.

A typical two reaction mechanism will be used to illustrate the method.

Next, the source term for the energy flux will be calculated, and

finally, the source terms arising from the diffusion of average molecular weight.

a) Calculation of the chemical source term

To show how the species source term is calculated, consider a

reaction mechanism consisting of two reactions

A+BiiC+D (l)

A+A+M^tE+M ,(2) (A4.2)

Reaction (1) is a second order reaction in both the forward and reverse

directions, while reaction (2) is third order in the forward direction

and second order in the reverse. These represent the two types of

reactions which we shall consider in the flame chemistry.

Consider the rate at which species A is formed. This is

[A] = -k^[A][B] + kJ][C][D] - 2k2[A]2[M] + 2k^[E][M]. (A^3)

The bracketed quantities signify concentrations in moles/m . We wish to

express ^(A) , the source term for species A, in terms of the concentration variables cp (moles/kg). From the table on page A15, we see

that

-^A ~ ^pj^^^k} kg"''s"''.

—3 —1

(Remember that R. is the mass production rate of A in kg m s .)
A22

From a dimensional analysis, we liave

= [A]/p, and [A] = pcp^. (AU)

Thus we need only divide A4.3 by p and substitute for each of the

concentrations their equivalent values in terms of cp; Ai43 becomes,

Let Us now calci:ilate S(species) = T|J^diu. To do this we will assume

that j^is constant throughout the control volume. Thus we can write,

S( species) = T]^(uj^ - au_) = Ueift"(^++ " ^__) = '^2^ (M6)

Note that we could have assumed that ^was linear between the grid

points as we did for the concentration variables cp (see Fig. A2). This

would have given (see A30),

S( species) - T] i(3f2^+ + fi_^__) (M?)

Thus A4.7 Uses the values of

^

at three grid points to approximate the

integral J^pjfduu, whereas A4.6 uses only the value at one point. We first

work out the implications of the one point approximation in terms of

our two reaction model A/^2. This is the approximation used in the computer program. Afterward, we shall show how the three point approximation

could be incorporated into the calculation.

1 ) One point approximation to the chemical source term integral.

As discussed earlier in this section, we want to express the

chemical source term as the sum of two terms,

S( species) = Sp + S^cp --^ {TTI^fi, (A^S)

here cp denotes the poststep value of the concentration of the particular

species, e.g., species A in A/+5; and we use A46 to approximate . S( species)

.

Let Us expand

^

in a Tayloi- series in cp about the point cpp, the prestep value.

/ e\ r P P

= Op + Opcp, (A^9)
A23

where

All these terms with the P subscript have prestep values.

As an example, let us calculate and a-p for the two reaction

mode] Ay, 2. Taking the derivative of A^5 with respect to cp^ we get,

and

2 2

cJp(A) = [kjjpcp^cpp + 2k^P cp^cpj^ + 2k^pcpgCf^]p. (A52)

The subscript P means that all of the quantities have prestep values.

For this case, we thus have,

S(A) = ^TlJ20p(A) + -^Tlfiap(A)(cp^)p = Sp(A) + Sj,( A) (cp^) , (A53)

where (^y^^j' is the poststep value of cp^.

2) Three point approximation to the chemical source term integral.

If we were to approximate S( species) by A47, the finite difference

coefficients of A4.0 for the species equations become,

A' = 2T^ - - i(P + G)Q^ + il1^F++"f

B' = 2T_ + L_ - i(P + G)fi_ + illap__fi_

C = iP(3cppn + 9p^fi+ + cpp__12j

+ 2S^( species) + fTla^,Q + iflap^fi,, ^ i'^o^_Jl_

D = 2(T^ + T_) t- - L_ f- f(P + G)fi - fTlap . (A54.)

These equations are not significantly more complicated than A40, and

should be more accurate. They have not, however, been incorporated

into the present computer prograjn.
A24

b) Calculation of the energy source term

The energy source term is easily calculated. From A35 and the

table on page AI5, we have

,2 _

J

(A55)

As Usual the ± subscripts refer to values at the control volume boundaries. For the first term we have

and for the derivative

- 9.

The - boundary is treated the same way.

c) Calculation of source terms arising from diffusion of molecular

weight

From A3$ and the table on page AI5, we have for the source terms

arising from the diffusion of the average molecular weight

MA.p2 ^1*^.

m

11 Z*^ dc = 1 lYh-A o^ll - ^-\y h-A p2jPi

J

- " <M>

for the energy equation, and

m

2 cp. ^ (M)'

<M> J+ ^L^ <M> J-

(A56)

(A57)

for the species equation.

It is worth noting at this point that we have been somewhat

careless about how we calculate quantities at the control voliome

boundaries. Throughout Appendix H and I, we have frequently been

confronted with expressions involving products of various quantities
A25

which are to be evaluated at the control volurae boundaries. For

example, let (<^J-|'l2'^3 ''^^^^^ ^^'^ \ various quantities

like p, h^, A-, X, etc. All of these quantities can be expressed in

J J terms of the values of cp. and h at the grid points; indeed, this is how

J they are evaluated at the end of every time step. This, of course,

yields their values at the grid points. The values at the control

volume boundaries are always taken to be the average of the values at

the two grid points on each side of the boi:indary in question. Thus, for

the + boundary we take

^k+ "

'^^'^k^ ''k^

To calculate we often take products of the q^^^.

Q+ = qT+q2+^3+ %+

Actually we should evaluate from the expression

= i-[(q^q2q3 %}++ ^ ^^1^2% ^^-^

If the q-^ do not change much from grid point to grid point, the two

ways of calculating Q_j_ would give about the same answer.

Let Us consider a case where Q is a product of two terms and compare

the two methods of evali-teting Q^.

Suppose that

q-| = q^ + Aq.^ and q^^^ - + ^<\2 •

The correct way to calculate gives

' 2(q^q2)-H- ^ hi(i^^i2^ - q-,q2 ^a(q-]Aq2 q2^^ii^ ^ -2^'i^^^2

The approximate method gives

Q^(appr) = hin^+++ q-,)i(q2+i- q2) qiq2

^ 'Hq/q^ + iAq^Aq^

Thus Q^(appr) - iAq^Aq2

For Aqj^ values encountered in the present calculations, the difference

between the two methods of calculating products at the boundaries is

negligible.
A26

APPENDIX J

CALCULATION OF THE ENTRAINMENT RATES

The so-called entrainment rates are the mass flow rates across

the hot and cold boundaries. In the steady-state flame, these will

naturally be equal to each other. During the calculation, however,

they are automatically adjusted to keep the grid centered on the flame

front. As discussed in Appendix F, the mass flow rate across the hot

boundary is

a = -d^-^/dt = m^^

Across the cold boundary it is

a + b = -d'/'o/dt = m^

Consider a j)Osition somewhere near the hot boundary; let cppjj^ be the

value of a concentration of one of the species at this point. Let cpQQ

be its value at a point near the cold boundary. The entrainment rates

are calculated from the following formiilas

mp - m

C m

^cc " ^C

xlOO

mu = m 2 - < H mL \

L % - ^c

" ^C

V

X100>

m = - Tl m '

' -^0 pM. /

- ^C

(A58)

(A59)

(A60)

where ^C and

911

are the values of cp at the cold and hot boundaries.

Note that cp^, cp^^, cp^^, and cpju^ have an implied subscript j which

denot(\s the species used to calculate the entrainment rates. In

principle, any species could be used because the integral in m^ will

be the same for all the species. However, it is best to choose one

of the major species for this purpose since the numerical evaluation

of this integral will not be as accurate if the species is a trace one.

To see how A$8 and A59 automatically control the entrainment rates.
A27

and thus the grid width, consider the situation in which the concentrations cpQQ and cpj^jj differed from (p^ and cp^^ by ^% of the total difference in cp across the grid. Then (cpqq - cpq)/(cPjj - cpq) and (cpj^ - cpj^j^)

/(cpj^j - cpq) would equal 0.01 and from A58 and A59 we have m^^ = m^^ = m^.

This would be the situation in the steady-state. Suppose this condition

arose by accident at some early time in the calculation. Suppose also

that as the calculation proceeded (cPqq - cPq)/(cPjj - cp^) increased. Then

would become somewhat greater than m^ and cold gas would be pulled

into the flame front at an increased rate. If m^ had not changed,

this would mean that the grid would start increasing in width. Increasing the grid width would tend to lower cp^^ for future times thus bringing mp back to m_ thereby slowing down the rate of increase of grid

yj m

size. Through the use of A$8 and A59, the grid is automatically adjusted after every time step in such a way as to make (cp^Q - cpq)/

(cpy - cPq) and (cpj^ - cpj^jj)/(cpg - 9q) tend to 0.01. This condition will

finally be reached in the steady-state when all the cp cease to change

with time.

The parameter v in A58 and A59 is used to control the effect that

departures in the 1% condition have on and m^^. For a stable calculation V must be of the order of 0.1. Note that v is piirely a numerical

device. We desire only to have some relation between "Iq^'j^ and the cp values such that ni^^jj -» m^ when the cp reach certain specified values.

Other relationships which served a similar purpose could no doubt have

been devised to perform the functions of A58 and A59.

It remains to show how the integral /(R./pMjdcju is calculated.

Recalling that R./pM, = (see page A21 )

, we used the three point

J J J approximation A<i7 to give

The jd. values Used in A6I are calculated from the values of cp at the

end of each time step.
A28

APPENDIX K

BOUNDARY CONDITIONS

Consider the finite-difference forms of the flame eqiiations in

final form. This is Eq. A43 derived in Appendix H.

^, = Af <^.,. + (A4.3)

1 1 1+1 1

^-^^

The subscripts refer to the grid points. There is one of these equations for each species and one for the enthalpy. To solve these

equations we must assign a value to ^^+2' grid point on the cold

side of the flame. Furthermore, to calculate B!; = B (J) + C we need

to specify a value for (I)^. This section is about how and 1^-^-2 are to be chosen. The values will depend on the boundary conditions.

We consider two types of flames; a) the freely propagating flame, and

b) a flame stabilized on a burner.

a) Bo'undary conditions for a freely propagating flame

On the cold side of the species equations we set 9jj^2 QQ^l "to

their values in the cold incoming gas. On the hot side, we assume

that the concentration gradient vanishes; thus we set cp^ = cp^. For

the enthalpy, we calculate ^-^^2 temperature and composition

of the cold gas mixture. Since the freely propagating flame is considered to be adiabatic, the value on the hot side is the same as that

on the cold side. Thus we set h^ - ^^+2*

b) Boundary conditions for a flame stabilized on a burner

The situation here is considerably more complicated than that

occuring in the freely propagating flame. Consider the cold side of

the grid. Locate the grid point N+2 at the surface of the burner.

Figure A3 shows this portion of the grid along with the burner and

the last control volume. The burner is considered to be made of some

porous material which allows free passage of stable species, but not

radicals. Cold reactant gases at a temperature T enter the porous
A29

material at a rate (kg/m s).

The composition and enthalpy of

the mixture at this point is cp-^^^

and h^. Since the porous plug

is absorbing heat from the flame,

there will be a temperature

gradient within it. The temperature of the gas will thus increase to some extent as it passes through the plug. It is also

conceivable that the concentrations

N

Figure A3.

cp. on the warm side of the plug may be somewhat different from cp. J J

°° because of diffusion into and from the flame zone. Thus on the warm

side of the plug which is located at the last grid point the gas properties have the values h^, and T^. These values are not known

beforehand, but must be determined by solving the flame equations.

Consider first the boundary conditions for the species equations.

Suppose Y.Q is the mass fraction of species j at the flame side of the

plug. The HB ss flux of j at this point will be (v + V.)„(pl.)„. This

J J must be equal to the rate at which j emerges from the plug. This will

be equal to the rate at which j enters the cold side of the plug, i.e.,

(pvY.) . As an additional source or sink for j we might also consider

radicals diffusing out of the flame zone and combining on the burner

surface. This effect is likely to bo of importance, however, only for

H atoms because of their large diffu^ilon coefficient. To include this

process in the model we define a surface distraction or production

rate for the species j. Call this If ®. is in units of moles/m s,

J p J then M.0. will be the mass rate in kg/m's. For the particular case of

J J H atoms recombining to give H^, we can let @jj =

~^j^(^f|)Qj where (C|^)q

is the concentration of H at the burner surface, and B^^ is a rate constant for surface recombination. We will let negative values of @
A30

correspond to destruction of the specie;;. For the formed as a result of this, we would have ®^ - 't^-^^^n^ q' general then, there

will be an additional source of j equal to -M.@,, and so we must have

J J (v + V ).(plj)p = (pvY )^ - M 0 (A62)

Since v is negative and denotes a flux from right to loft along the

y axis, -M.0. for 6. > 0 also gives a negative flux which corresponds

J J J to production of j

.

Since pv = m^ throughout the flame (note that will be specified

as part of the experimental conditions), A62 becomes

m I.p + (pV.Y.)^ = m Y. -M.0. (A63)

Inserting the expression for V. (Eq. A5) into A63, we get

Y,

„ = m (Y.„ - Y. ) + M.0.

Using the identity shown on page A9 and dividing the resulting equation by M. gives

J

If the gradient of (M> is neglected, this becomes

Transforming to the u),t system with the relation ^/^y = (p/Tl)VbuJ gives

Since the cp profile is assumed to be linear between grid points

(see Fig. A3), the concentration gradient at point C, i.e., i-N+2, will

be the same as that at the mid-point between i=N+1 and N+2. Thus

/^cp.\ /bcpA .

- cp. . ^^^2 - ^jN^I (A65)

Inserting this into A6/4, and solving for cp.-vr^p ~ ^-c ^®
A31

L^^wvi) ^"^

^^J

- -7--i- N ''co (A66)

At the hot boundary we assume ^cp/bou = 0, and thus we set cp^ ~ in

the freely propagating flame.

Next we consider the boundary conditions on the energy equation.

The net flxix of energy at the burner surface will be the sum of that

caused by the mass flow of the reactant mixture, the flow caused by

diffusion of species to and from the surface, and the flow arising

from a temperature gradient in the gas at the surface. Let (F']^)q be

this net energy flux at the burner.

= -(^a|)^*-„\+Z/Wj'o (167)

Introducing the expression for V. and the identity on page A9, neglecting the gradient of (M), converting to cp. concentration units, and

going to the ou,t system, we get

Since this is the net flux of energy into the cold side of the flame

front, it must be equal to the net flux out of the flame on the hot

side. Because we assume that "bT/by and"bcp./by = 0 there, this flux

out is ^J^^ where is "the enthalpy of the hot gas mixture. Assuming

a linear variation of h between grid points gives

/bh\ ^h\ _ ^

N+2 ~ _
A32

Inserting this into A68, and setting (F^)q ~

^-qo^h'

Si^es

This equation can be solved to give ^-^^2 ^~^C^ terms of ^^^^ and

hj^, and other known quantities. The enthalpy at the hot side can be

calculated from the measured final flame temperature and its composition. From h^ and the calculated composition (cpj)Q at the burner surface, we can also evaluate Tp.
A33

APPENDIX L.

TRANSFORMATION FROM THE uo,t COORDINATE SYSTEM BACK TO THE y,t SYSTEM

The steady-state solution of the species and energy equations A22

and A23 in the u), t coordinate system yields values of cp. and h as functions of lu, where 0 ^ co ^ 1 . If we want to compare these profiles with

experimental ones we must express cp. and h in terms of y, the laboratory

spatial variable.

The variable cp is a function of 'A and t, and ^is a f unction of y

and t; thus,

u) = (d( \^,t) = u)( ^{y,t),t)

A small change in ^ can be expressed as

d^ = ^

dy +

^

dt = pdy - pvdt

Here we have used A18, the expressions defining ^. A small change in

u) can likewise be written

dcu = ^

d^ +

^

dt = ^j£(pdy - pvdt) +• ^uu dt

^t ^t

~ ^

pdy + (-pv 3a) + ^

)dt

(^4/ 3\A ^t

Using the facts that

^ ~ r and

^ ^ — 1

^"^

, we get

^ 3t 'I

da)- ^dy M-£V ^^L^ ) dt

In the steady state b ~» 0 and a m^ = pv. Therefore duj -» ^

dy.

Solving this for dy and integrating from a) = 0 (the hot side) to some

arbitrary o) value, we get

In the present calculation, we arbitrarily set y^j = 0. Thus we can

calculate the value of y which corresponds to a particular value of ^

by evaluating the integral
A34

This integral is approximated in the computer program by the summation

(A70)

where the subscripts refer to the grid points.
A35

APPENDIX M.

CALCULATION OF THERMAL AND TRANSPORT PROPERTIES FOR PURE SPECIES

An expression for the effective diffusion coefficient A . of a

J particular species in the flame has been derived in Appendix C (eq.A9).

This gives A- in terms of the concentrations of the various species

and the binary diffusion coefficients for all possible pairs of species.

In practice, the summation in the denominator of A9 need not be

evaluated over all the species since x^, the mole fraction, is large

only for the major species. Thug, we must calculate only the binary

diffusion coefficients between major species, and between major and

minor species. Those coefficients involving pairs of minor species

need not be evaluated. In Appendix D we showed how the thermal

conductivity of the flame was calculated. The expression used contains

the thermal conductivities of the pure species, and also their heat

capacitiej . The heat capacities of the pure species are also needed

to calculate C^, the average heat capacity of the flame. This quantity

occurs in the energy equation. For the energy equation, we also require

values for the enthalpies of the pure species. Thus, we must evaluate

the following properties involving the pure species (or pairs of species

in the diffusion case):

a. Binary diffusion coefficients D. ..

b. Thermal conductivities \..

J c. Heat capacities C .. PJ

d. Enthalpies h. (or M^^h.M.).

a) Binary diffusion coefficients

For the diffusion coefficients, we used Lennard-Jones potentials.

This potential function has two parameters a. and e ./k. a. gives the

a

J J J distance in a where the potential goes from repulsive the attractive,

and e ./k in °K gives the depth of the potential well. a. and e . refer

il J J the the interaction between two molecules of the same species. For

the diffusion coefficients, we require analogous parameters for the

interaction between molecules of different species. If a., o. and

e., e. are the parameters for the species i and j, we take as the
A36

parameters for the interaction between i and j to be.

CT. . = -g-Ca. + a.) and e . . = -Je .e ..

The formula for the binary diffusion coefficient is

D.. = 1-86x10-^ [(M. .M.)/M.M,]^T^ -'-^ ^^^^j

where P is the pressure in atmospheres, M. and M. are the molecular

O ( 1 1

)

weights, and . ' is a fi-mction of a reduced temperature T^'-^TkA.

To a good approximation this function is proportional to T where a

takes on two constant values, one in the range T ^ 3^ . ./k and another for

T< 3e . ./k. We have calculated D, . from either of the following formulas

..„.a.. » ». ....i. „ ....

„ I.MO-' [(«, .»,|/.,«jtT'-'' I, 5.,.A (.7!!

Pa?, (e. ./k)°''^'^

^ , 1 .26x10-^ [(M. + M.)/M.M.]i t1 ^ ^ ^^^^^

^ J - .- I »i I

I I ! .1 I I n il II 11 Pa?. (e../k)°*^

These formulas will give D. . in units of m /s.

b) Thermal conductivities

Lennard-Jones potentials can also be^sed to calculate the thermal

conductivities of the pure species. The formula for the thermal

conductivity is similar to that for the diffusion coefficients. In a

like manner we have found that it is possible to express it approximately

in terms of a constant power of T, where the appropriate value of the

exponent depends on whether the temperature is above or below 3^ ./k.

J The formulas obtained were,

, , i.6o^x10-Vt° -67 ^^^^^^^ ^^^^^ a2(e,A)°-"
A37

' a' (e./k)°-"

'

J J

These two formulas give the thermal conductivities in units of cal/cm-deg-s.

To be Usable in the program they would have to be converted to J/m-deg-s

units.

In the present version of the program, experimental thermal conductivity data were Used for the major species in the ^^-O^ flame.

The data were fit by a least squares calculation to a power series in

T and the least squares coefficients were used in the program to calciilate

the X. values.

J

c) Heat Capacities

Heat capacity data from the JANAF Tables were fit by least squares

calculations to a power series in T. Terms up to the second power were

Used. Thus, heat capacities were calculated from the formula

C . = d. + e.T + f .T^ (A76)

PJ J J J where d., e., and f . have values such that the C . are in units of

j' J J PJ

J/kg-deg.

d) Enthalpies

The enthalpy of a species in J/kg is given by the formula

/•T

h. - h^ + /n, C ,dT (A77)

J J '•j^ PJ

R

where h. is the enthalpy at the reference temperature Tp.

Substituting A76 , into A77, we get

h. = (h^? -d.Tj^ - ^e.T^ - If .T^) + d.T + ie.T^ + if .T^ (A78a)

= h^ +d T + ie T^ f if T^ (A78b)

This formula was used in the program to calculate h. values. Values of

o R

h. were calculated from h. values taken from the JANAP Tables, the T^^ values,

and the heat capacity coefficients d., e., and f..
A38

APPENDIX N.

INITIAL PROFILES

The mixture of cold gases through which the flame will propagate

is naturally therraodynamically unstable. Nevertheless, in most systems

of practical interest the rate of approach to equilibrium at low temperatures is negligibly slow. This is because the free radical reactions

which destroy the reactants normally have large activation energies.

There are two ways to speed up these reactions. One can raise the

temperature of the system or introduce a sufficiently high concentration

of the necessary radicals. For a flame steadily propagating into a cold

unstable gas mixture, both of these processes occur. The cold gases

receive heat by thermal conduction and radicals by molecular diffusion

from the flame region. To solve the time dependent flame equations,

one must begin with the system in such a state that the radical reactions

are proceeding at an appreciable rate. Furthermore, it is desirable

to start the calculation with concentration and energy profiles which

are similar to those expected in the final steady state to minimize the

number of steps required for the integration.

These goals were satisfied by starting the calculations with

S-shaped profiles for temperature and for the concentration of the major

species. These profiles were obtained by using the followirg function of

cu, the spatial variable in the Spalding coordinate system:

P^((u) = lOu)^ - 15/ + 6isP (A79)

This has the boundary values, = 0> and ) - 1* Since uu = 0

corresponds to the hot side of the flame, A79 represents a decay profile.

Formation profiles were obtained from the function

P^(aj) = 10(1 -(ju)^ - 15(1-0))'^ + 6(1 -u))'^ (A80)

At the start of the calculation, we specify the mole fractions of a major

species on tho hot and cold sides of the flame. The former is calculated

by nssujriing tliat tho overall reaction has gone to completion. Also

specified is whether the species decays or grows as we go from the cold

to the hot side of the flame.
A39

For example, suppose the molo frnction of species j is x.p and

x.tT for 0) = 1 and 0, and also suppose that x.^ > x.^, i.e., we have

a decay profile; then

x.(u)) = x.„ + (x.^ - x.JdOuP - 15(1)^ + 6cu^)

J J*^

To calc'ilate an initial temperature profile, we first calculate the

enthalpy on the cold side from the cold side tempera tiJre with the

expression,

is the concentration of j on the cold side in units of mole/kg.

We then set the enthalpy on the hot side equal to that on the cold

side. (Note that jj^ is the index of the last radical species; i.e.,

we sum over the major species only.) The expression for the hot side

enthalpy is

r = . a.. + d.T„ + ^e.T^ + %f -^2) (A82)

= A + BTjj + JCT^ + §DT^

where

M.

J

M.

J V'.^d.

M.

J

M.

J

We want to calculate T^^ in terms of h^^. This was done by finding the

appropriate root of the equation

F(Tjj) = A +BTj.j + ^-CT^ + iBT^ ' " °

A good first approximation to the desired root can be gotten by neglecting

( 1

)

the non-linear portion of A82. This gives T^'^ - (h^^ - A)/B.
A40

The Newton-Raphyon method was used to got u better approximation.

The second approximation is

m(2) _ A^) _

) w„, )

N

H

^

H ^/ ^ H

^

"4^^ + |h^ - (A + BT^''^ + ^CT^'^^^+ iDT^^^^)|/(B + CT^^ ^ + DT^^

(3) In practice, T-^ was a sufficiently acciirate value for T-^,

The Newton-Raphs on method was also used in the main part of the

program to calculate the temperature from the enthalpy value at a particular grid point.

An S-shaped starting temperature profile was calculated from the

formula

T(ci)) = T^ + {T^ - T(.)(10{1-U))^ - 15(1-0))'^ + 6(l-uj)5)

For the enthalpy, a constant starting profile h(uo) = h^ was used.

It should be noted that this profile is not consistent with the starting

temperature profile. The enthalpy profile should be calculated from

the temperature and concentration starting profiles. It would be

worthwhile to examine what effect, if any, the starting enthalpy profile

has on the course of the calculation.

The radical concentrations were set to zero on the cold boundary

and remained so throughout the calculation, since the cold boundary

concentrations of all species are never modified. Throughout the rest

of the flame they initially were set to some constant value comparable

to their average expected steady-state concentrations. For flames like

0^ decomposition and - Br^ the initial radical concentrations can be

set to zero throughout the flame zone. This is possible because the

radical generating reactions like Br^ + M ^ 2Br + M proceed rapidly

on the hot side of the flame. For the - 0^ flame however, such

reactions are too slow for cool flames to initially generate enough

radicals. The flame dies out before the radical concentrations become

high enough to sustain the flame. Thus, one must begin the calculation

with non-zero radical concentrations.
NBS-114A (REV. 7-73)

U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET

1. IHIHI.K A riON OR RI'I'ORT NO.

NBSIR 76-994

No.

3. Rcc ipicnt '.s Accession No.

4. in i-H AND SlIIVl lTI.l'.

A Numerical Solution of the Time Dependent Partial

n-i T" on t"1 fl 1 Fmiat'lons wTilrh D©j3f"rll>es a one—dimensional.

Laminar, Premioced Flame

5. Publication Date

January 1976

6. Performing Organization Code

7. AJJTHI)R(SJ

R. L. Brown

s.V . r^prfrirmino 1 1 1 1 1 1 1 1 1 ()rpan 1 fL nil. Reriorr l^v^l/ \J 1 L No1 . yj . NBSIR 76-994

9. PERF^RMWA' ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

10. Project/Task/Work Unit No.

3162675

11. Contract/Grant No.

12. Sponsoring Organization Name and Complete Address (Street, City, State, ZIP) 13. Type of Report & Period

Covered

Final

14. Sponsoring Agency Code

15. SUPPLEMENTARY NOTES

16. ABSTRACT (A 200-word or less factual summary of most si^ificant information. If document includes a significant

bibliography or literature survey, mention it here.)

The set of time dependent, parabolic differential equations, which describe the

physical and chemical processes in a one-dimensional, laminar, premixed flame is

solved by adapting a solution procedure originally developed to solve the twodimensional steady state boundary layer equations. The flame equations are

integrated by an implicit method until the steady state is reached. This corresponds

to a flame propagating steadily through a mixture of combustible gases. By a

suitable choice of boundary conditions, it is also possible to model a flame

which is stabilized on a burner. Solution of the flame equations yields the

concentration profiles of the different chemical species as well as the temperature

profile. From these one can also calculate the production rates of each species,

the rate of each chemical reaction, and the heat release rate at each point in

the flame. The velocity of the freely propagating flame can be calculated from

the integrals over the whole flame zone of any of the species production

rates. The model incorporates realistic thermodynamic data and transport property

data that are functions of both temperature and concentration. A complete

documentation of the computer program which accomplishes the integration is

presented.

17. KEY WORDS (six to twelve entries; alphabetical order; capitalize only the first letter of the first key word unless a proper name; separated I'y semicolons) Combustion; differential equations; flame chemistry;

flame modeling; laminar flames; numerical solution of flame equations

18. AVAILABILITY [X Unlimited 19. SECURITY CLASS

(THIS REPORT)

21. NO. OF PAGES

For Official Distribution. Do Not Release to NTIS

UNCLASSIFIED 175

Order From Sup. of Doc, U.S. Government Printing Office

Wa.shinRton. D.C:. 20102. SD Cat. No. CH

20. SECURITY CLASS

(THIS PAGE)

22. Price

[IXl Order From National Technical Information Service (NTIS)

Springfield, Virginia 221^1 UNCLASSIFIED $6.75

USCOMM-DC 29042-P74