NASA

Technical

Paper

3315

1993

National Aeronautics and

Space Administration

Office of Management

Scientific and Technical

Information Program

A Critical Analysis of

the Accuracy of Several

Numerical Techniques

for Combustion Kinetic

Rate Equations

Krishnan Radhadrishnan

Lewis Research Center

Cleveland, Ohio

Summary

A detailed analysis of the accuracy of several techniques

recently developed for integrating stiff ordinary differential

equations is presented. The techniques include two generalpurpose codes, EPISODE and LSODE, which were developed

for an arbitrary system of ordinary differential equations, and

three specialized codes, CHEMEQ, CREK 1D, and GCKP84,

which were developed specifically to solve chemical kinetic

rate equations. The accuracy study is made by applying these

codes to two practical combustion kinetics problems. Each

problem describes adiabatic, homogeneous, gas-phase chemical

reactions at constant pressure and includes all three combustion

regimes: induction, heat release, and equilibration. To

illustrate the error variation in the different combustion

regimes, the species are divided into three types: reactants,

intermediates, and products, Error versus time plots are

presented for each species type and the temperature. These

plots show that CHEMEQ is the most accurate code during

induction and early heat release. During late heat release and

equilibration, however, the other codes are more accurate. A

single global quantity, a mean integrated root-mean-square

error, that measures the average error incurred in solving the

complete problem is used to compare the accuracy of the

codes. Among the codes examined, LSODE is the most

accurate for solving chemical kinetics problems. It is also the

most efficient code. in the sense that it requires the least

computational work to attain a specified accuracy. An

important finding is that use of the algebraic enthalpy

conservation equation to compute the temperature can be

more accurate and efficient than integrating the temperature

differential equation.

Introduction

Many practical problems arising in chemically reacting flows

require the simultaneous solution of large sets of coupled

ordinary differential equations (ODE's) which describe the

time rate of change of chemical species concentrations and

temperature. Examples of such problems include the development and validation of reaction mechanisms, combustion of

fuel-air mixtures, and pollutant formation and destruction.

The main difficulty in using classical methods, such as the

popular explicit Runge-Kutta method (e.g., ref. l), to solve

large sets of chemical kinetic rate equations is that of

"stiffness." The property of stiffness arises in chemical

kinetics because of the widely varying time constants fi_r

different species. For free radicals the relaxation time is on

the order of microseconds, whereas the nitric oxide R)rmation

time is on the order of milliseconds. To satisfy the stability

requirements that errors in the numerical solution remain

bounded as the calculation proceeds in time, classical methods

must use extremely small step sizes, as illustrated in

references 2 and 3 for the explicit Runge-Kutta method in

solving combustion kinetics problems. Consequently, these

methCxls require prohibitive amounts of computer time to solve

a practical chemical kinetics problem.

Numerous approaches have been proposed for stiff ODE's

to remove the stability restriction on the step size. In Part I

of this effort (ref, 2) and other recent publications (refs. 3 Io 5),

several techniques were examined, and detailed comparisons

of their computational work requirements fi)r solving combustion kinetic rate equations were made. The methods

examined in these studies include the general-purpose packages

EPISODE and LSODE (refs. 6 to 9), which were developed

for an arbitrary system of ODE's, and the specialized codes

CHEMEQ (ref. 10), CREKID (refs. l I to 14), and GCKP84

(ref. 15), which have all been developed specifically to

integrate chemical kinetic rate equations. In the present work

the accuracy of these techniques in solving combustion kinetic

rate equations is examined.

In general, numerical methods generate approximate

solutions to the governing ODE's at discrete points in time.

To maintain accuracy of the numerical solution, they require

that the estimated error incurred on each time step be less than

a user-specified local error tolerance. This result is usually

achieved by restricting the size of the time step. Some solvers,

in addition, adjust the order of the numerical approximation

when appropriate. In either case, only the estimated local error,

that is, the estimate of the error incurred in advancing the

numerical solution by one time step, is controlled. However,

the quantity that is of interest to the user is the global error,

which is the deviation of the numerical approximation from

the exact solution and which generally accumulates in a

nontrivial manner from the local errors.

In the present paper, a detailed study of the estimated global

error incurred by the above techniques in solving combustion

kinetic rate equations is presented. Also presented is a study

of the variation of the global error with the user-specified local

tolerance and an examination of the computational cost,

measured by the required CPU execution time, associated with
attainingdesiredaccuracy.Thepaperconcludeswithtwo

appendixes: Appendix A describes the methods examined in

this study, and appendix B describes the procedure used to

solve the algebraic enthalpy conservation equation for the

temperature.

Symbols

Aj,A_) pre-exponential constants in tbrward and reverse

rate coefficients for reactionj (eqs. (6) and (7)),

units depend on reaction type

local absolute error tolerance for ith component,

required by LSODE (eq. (20)).

local absolute error tolerance used with LSODE

for all species mole numbers

exponent-on-ten in pre-exponential constant for

forward rate coefficient of reaction j, where Bi =

Iogl0A;, arbitrary units

local error test constant used in GCKP84 (eq. (21))

constant-pressure molar specific heat of species

i, J/kmole K

estimated local truncation error in i _hcomponent

at l,,

cumulative difference between converged and predicted values of (dYi/dt) at t,,, used by GCKP84,

units depend on component i

activation energy in forward and reverse rate

coefficients for reactionj (eqs. (6) and (7)), cal/mole

mean integrated root-mean-square global error

(eq. (34))

for EPISODE and GCKP84: local relative error

tolerance for species with initially nonzero mole

numbers and temperature, and local absolute

error tolerance for species with initially zero mole

numbers; for LSODE: local relative error tolerance

for all components, for CHEMEQ and CREK 1D:

local relative convergence criterion for all

components.

relative error tolerance for Newton-Raphson iteration

for temperature

local error weight for ith component, used by

LSODE (eqs. (19) and (20))

estimated global error in ith component (eq. (24))

estimated global error for ith species of type j

root-mean-square norm of the estimated global

errors for all variables (eq. (27))

root-mean-square norm of the estimated global

errors for species of type j (eq. (26))

estimated global error in temperature (eq. (25))

ATOL,.

ATOLSP

Bi

CG

di. tl

Et,I1

Ei,E- )

EPS

ERMAX

EWT,

ei

ei ,.i

eFIlll_

'Cr uv",,I

_"7'

f

/-/O

h_

h0

IERROR

ITOL

kj.k_j

N

N ,N_j

NR

Ns

PI ,P2

P

R

Ra

Rj,R_)

RTOL

T

AT

TSITINY

t

tend

it1

to

V

Xi

xi

Xi,xt

Xmin

Y).,,

ylml

i.ll

generalized algorithm for Y,, (eq. (14))

time rate of change of i _hcomponent, units depend

on i

initial mixture mass-specific enthalpy, J/kg

molar-specific enthalpy of species i, J/kmole

step size used on the n th step, s

initial step length to be attempted by integrator, s

error control indicator for EPISODE

error control indicator for LSODE

forward and reverse rate coefficients for reaction

j (eqs. (6) and (7)), units depend on reaction type

total number of first-order ordinary differential

equations

temperature exponent in forward and reverse rate

coefficients for reaction j (eqs. (6) and (7))

total number of elementary chemical reactions in

reaction mechanism

total number of chemical species in reacting gas

mixture

test problems I and 2, respectively

pressure, N/m 2

universal gas constant in cal/mole K

universal gas constant in J/kmole K

molar forward and reverse rates per unit volume

for reaction.j (eqs. (4) and (5)), kmole/m 3 s

local relative error tolerance for LSODE

temperature, K

maximum temperature change allowed before

reaction rate coefficients and thermodynamic

properties are updated in CREKID, K

standard solution value for temperature, K

minimum species mole number values allowed in

CHEMEQ and CREK 1D

reaction time, s

final time (_> I ms) at which numerical solution

is generated, s

time reached on the n _hintegration step, s

initial time, s

reacting gas velocity, m/s

chemical symbol for i th species

mole fraction of species i

standard solution value for mole fraction of species i

mole fraction value corresponding to omi. (eq. (33))

numerical solution of the ith component at t,,, units

depend on i

value obtained for Yi.,_ on mth iteration, units

depend on i
yIm}

i.n

yPq

i,Ii

Yil,

VIlla x.l

yi

O_

(-jl

u'_;,, % ,,

p

0 i

_min

value obtained tbr Y,.,, on mth iteration, units

depend on i

predicted value of Y, .... units depend on i

numerical solution of the ith component at 6,

generated by using exact past valucs, units depend

on i

local weight fi)r im component, used by EPISODE

and GCKP84 (eqs. (17), (18), and (21))

cxact solution for the i th component, units depend

on i

constant in gencralized algorithm for Y,, (eq. (14))

global error at t,, (eq. (15))

stoichiometric coefficients of species i in forward

and reverse reaction j (eq. (1)); numbcr of

kilomoles of species i in elementary reactionj as

a reactant and as a product, respectively

mixture mass density, kg/m 3

mole number of species i. kmole species i/kg mixture

mole number value at which local error control in

LSODE is equally relative and absolute (eq. (32))

Governing Differential and

Algebraic Equations

The ordinary differential equations describing homogeneous

gas-phase chemical reactions of the type

Ns Ns

i=l i:I

j = 1..... AiR (i)

are as follows:

do t

-- = f(Ok,T )

dt

i,k = 1..... N, (2)

oi (t0) = given

T(to) = given

where f, the total formation rate of species i, is given by

N S

/

R, = _',H (poa'"

/=1

(4)

and

N5

R i=k j_I(pol)";'

/=1

(5)

where the forward (k/) and reverse (k_/) rate coefficients are

given by the modified Arrhenius expressions:

ki = AJ '% exp (-E/RT)

k w = A_iT '_' _ exp (-E_//RT)

(6)

(7)

In equations (I) to (7), v,'i and e,; are the stoichiometric

coefficients of species i (with chemical symbol X,) in reaction

j as a reactant and as a product, respectively; N s is the total

number of distinct chemical species (reacting and inert) in the

gas mixture: N,,¢ is the total number of independent reactions

in the mechanism: di is the mole number of species i (in

kilomoles of species i per kilogram of mixture): t is the time

(in seconds): p is the gas mixture mass density (in kilograms

per cubic meter): T is the temperature (in kelvins); R is the

universal gas constant (in calories per mole per kelvin); and

Ai, A_i, N_, N_j, Ei, and E_j are constants in the modified

Arrhenius expressions for /,) and k_i. The reverse rate

coefficient parameters are calculated from the forward rate

coefficient parameters and the concentration equilibrium

constants by using the principle of detailed balancing (ref. 16).

In this paper, as in the companion paper (ref. 2), attention

is restricted to adiabatic, constant-pressure chemical reactions.

For such problems, the following enthalpy conservation

equation constitutes an algebraic constraint on equations (2)

to (7):

N A

Z oig_

i=1

= Ho, (8)

where gi is the molar-specific enthalpy of species i (in joules

per kilomole) and /4o is the initial mixture mass-specific

enthalpy (in joules per kilogram). Equation (8) can be

differentiated with respect to time to give the following ODE

tbr the temperature

N8

: - P-' E (v,) - vij)(R j - R_))

j=l

(3)

The molar reaction rates per unit volume. Rj and R_j, are

given by the law of mass action (e.g., ref. 16):

Ns

'_ f_i

dT i= 1

dl N,

E OiCp,i

i=l

(9)
where _),., is the constant-pressure molar-specific heat of

species i (in joules per kilomole per kelvin). Either cquation

(8) or (9) can be included in the equation set. We explore the

use of both these equations and examine their effects on

solution accuracy and computational cost.

The mass density of the mixture is given by the ideal gas

equation of state

p =pI(R.To,,,) (10)

where p is the absolute pressure (in newtons per square meter),

R,, is the universal gas constant (in joules per kilomole per

kelvin), and _r,,,, the reciprocal of the mean molar mass of the

mixture, is given by

0,<#1= _] a, (11)

i=1

For constant pressure problems the following density ODE

can be obtained from equation (10) by differentiating it with

respect to time and then rearranging terms in the resulting

expression:

dt <,,,,,=

12)

Either equation (10) or (12) can be used to compute the density.

The code GCKP84, which allows the pressure to vary, solves

for p by integrating its ODE (eq. (12)). With the other codes,

however, we obtain p by using equation (10). Indeed, ,o is

implicitly replaced by the right-hand side of equation (10) and

does not appear as a variable. We, therefore, exclude density

from our discussion, including statement of the problem, and

restrict attention to solving for the other Ns + I quantities.

Problem Statement

The initial value problem may be stated as follows: Given

(I) at time t -- to, values for the species mole numbers, oi

(i = 1..... Ns), and the temperature, T, (2) the pressure, p,

and (3) the reaction mechanism, find, at the end of a prescribed

time interval, the mixture composition and temperature.

Methods and Codes Examined

The codes examined in this study include the general-purpose

packages EPISODE and LSODE (refs. 6 to 9) and the

specialized techniques CHEMEQ (ref. 10), CREK1D (refs.

I I to 14), and GCKP84 (ref. 15). The methods used in these

codes are summarized below and are discussed in detail in

appendix A.

The packages EPISODE and LSODE consist of a variableorder, variable-step implicit Adams method (suitable for

nonstiff problems) and a variable-order, variable-step backward

differentiation formula method (suitable for stiff problems;

e.g., refs. 1 and 17). Both methods use a standard predictor

and a variety of corrector formulas--from functional iteration

to a modified Newton iteration--is included. The Jacobian

matrix Of/03,, where 3' is the vector of dependent variables

andf = dy/dt, is computed either numerically or with a usersupplied subroutine. In part I of this investigation (ref. 2) all

options relevant to the problem of chemical kinetics were

attempted, and the stiff method with Newton iteration and usersupplied analytical Jacobian matrix was found to be the fastest.

Therefore, only this option is used in examining the accuracy

of EPISODE and LSODE.

The general chemical kinetics program GCKP84 uses the

integration technique developed by Zeleznik and McBride

(ref. 18). The algorithm is essentially a revised version of the

GEAR package (ref. 19), which contains the same two integration

methods as EPISODE and LSODE and several iteration

techniques. For reasons given above we restrict attention to

the stiff method with Newton iteration using an analytical

Jacobian matrix. GCKP84 includes corrective actions if thc

physically impossible situation of negative concentrations,

temperature, density, or velocity arises.

In CHEMEQ, at the start of each time step the ODE's are

separated into two classes: stiff and normal. For equations

classified as normal, a classical second-order predictor-corrector

method, the trapezoidal rule, is used. For the stiff equations

a simple stable asymptotic integration formula is used.

The code CREKID is based on the exponentially fitted

trapezoidal rule developed by Liniger and Willoughby (ref. 20)

and Brandon (refs. 21 and 22). This code includes special

treatment of ill-posed initial conditions and automatic selection

of Jacobi-Newton iteration or Newton iteration.

Error Considerations

In this section the error controls used in the different codes

examined are discussed. In general, numerical methods replace

the differential equations with difference equations and solve

them step by step. Starting with the known initial conditions

y(t o) at t0 numerical approximations Y,, to the exact solution

y(t,,) of the ODE's are generated at discrete points in time

(t,, (n = 1,2 .... )), until the end of the integration interval is

reached. At each 5, the numerical method provides a rule tbr

generating the approximate solution Y,, in terms of computed

quantities at one or more previous times.

For the scalar differential equation,

dv 1

dtt = f(y)

(13)

y(to) = given
thealgorithmsusedforY,, in all of the codes can be written as

Y,, = _Y,,_j + h,, _(h,,,Y,.];,,Y.-I,Z,-_ .... ) (14/

where o_is a constant, Y,, _y is the approximate solution at t,,_y,

h,, (equal to t,,- t,,__) is the step size used on the step

[t,,_ +,t,,] and f,;_.] = f(E;-/). Because equation (14) involves

the unknown quantity Y,; its solution generally requires an

iterative procedure. Starting with the predicted value (an initial

guess), denoted by yI,ol improved estimates El,;'';) (m = 1,2 .... )

are generated until the iteration converges, that is. until the

difference in two successive approximations approaches zero

within a specified accuracy.

During the calculation procedure, errors called discretization

or truncation errors are introduced into the numerical solution

because of the approximation of the ODE's by difference

equations. Two measures can be defined for this error, which

is a property of the numerical method (e.g., ref. 23). The

global discretization error e,, at any t,, is the difference

between the computed approximation y,, and the exact

solution y (t,,) :

_,, = Y,, - y(t,,) (15)

It is the quantity that the user wants to know and control. The

local truncation error d,, at t,, is the error in the numerical

approximation Y,] that is generated on the step [t,_l,t,,] by

using exact past values:

d,, = Y;, - y(t,,) (16)

It is the quantity that ODE solvers generally control. The two

discretization errors are illustrated in figure 1 for a single

ODE.

The codes examined in this study require the user to specify

values for one or more local tolerance parameters, which

control the accuracy of the numerical solution. Now, as

discussed below, the same error control is not used in all codes.

Nevertheless. for convenience, for all codes the same notation,

EPS, is used to denote the local tolerance quantity, or the

primary one if several are required. In EPISODE the local

truncation error vector d,, satisfies the inequality

di tl -

,=,l<)) (17)

where N is the number of ODE's, di.,, is the estimated local

truncation error in the ith component at t,,, and for the error

control used

Yllll_a x i

for i that satisfy y,(t0) # 0,

for i that satisfy yi(tu) = O,

(18)

y,Y

r.

Y(t n)

Yn-1

Y(tn-1 )

Y2

Y1

y<to),v0

Exact solution

0, • Numerical solution

,r-"/,

t_

-, i

) I I

J I _ _m-" 'n-_l .

t o t 1 t 2 tn_ 1 t n t

Figure l.--Nutnerica[ solutions and truncations error types for the single ODE

dy/dt = ./(3'). The exact solution to the ODE is denoted by y(t). The

numerical solutions oblained with the initial condition Yo -- Y(to) are

denoted by solid circles. The solid square denotes the numerical solution

obtained at t,, by using exact past values. The local truncation error is

denoted by d, the global truncation error by c. and the step length by h.

where the vertical bars denote absolute value. The error control

selected to be performed by LSODE is given by

,:l \ewr# /

(19)

where

EWT_ = EPS iYi.,,-tl + ATOL_ (20)

where ATOLi is the user-supplied local absolute error

tolerance for the ith component.

In GCKP84 the local error test satisfies the inequality

,=, _ co EPs (21)

where, El.,, contains the cumulative difference between the

converged and predicted values of the derivative (dYi/dt) at

t,, and where Ca is a constant. The quantity Y,_,,_.i has the

same meaning as in EPISODE (see eq. 18)).

The codes CHEMEQ and CREKID do not control the

estimated local truncation error. The solution is accepted when

the magnitude of the normalized difference in successive
estimates(_,,,+11_ EI,,,,I) is less than a specified amount.

Therefore, these codes control only the error in the solution

to the difference equation of the method. In CHEMEQ each

component i satisfies the inequality

iyz,,,+'J i.II -- _,,,l t,_t i

min tlYi.l;"+_li, IY,,II"II}

< EPS (22)

The convergence criterion used in CREK1D is given by

t) 1 ll<l,,,+_J- _!;;,i 2 -

_i.,, __ < EPS (23)

i-- vitYJ -

It is clear from the above discussion that the user-supplied

local tolerance EPS does not have the same meaning for all

codes. In LSODE it is the local relative error tolerance for

all variables and is a measure of the number of accurate

significant figures in the numerical solution. In EPISODE and

GCKP84, however, as discussed in the section "Computational

Procedure," EPS is the local relative error tolerance for only

variables with nonzero initial values, such as the temperature.

For species with zero initial mole numbers EPS is the local

absolute error tolerance and is a measure of the largest number

that may be neglected. In contrast to these three codes,

CHEMEQ and CREK1D do not control the local truncation

error, and EPS is the local relative convergence criterion, or

error in the solution to the difference equation. However, as

described in appendix A, although CREK1D does not test that

the estimated local truncation error is within a prescribed

bound, the step length calculation procedure attempts to

achieve this result. The step length to be attempted next is

selected such that the current estimate of the local truncation

error normalized by the solution is at most equal to EPS.

Because of these differences in the meanings of EPS it will

be referred to as simply the local tolerance.

Evaluation of Temperature

Of the codes examined in the present study, only GCKP84

and CREK I D were written explicitly for nonisothermal chemical

reactions. These methods, therefore, have built-in procedures

for calculating the temperature. For the other codes, however,

the temperature has to be calculated along with the mixture

composition. In the present study (as in ref. 2), the temperature

was computed using one of the two methods outlined below.

In method A the temperature was calculated from the initial

mixture mass-specific enthalpy /4o and the solution for the

species mole numbers returned by the integrator by using the

algebraic enthalpy conservation equation (8). This equation

was solved for the temperature by using a Newton-Raphson

iterative technique, with a user-supplied local relative error

tolerance, ERMAX (as described in appendix B). In this

method, the temperature is not an explicit dependent variable,

so the number of ODE's is equal to the number (Ns) of

species and the Jacobian matrix is of size Ns x Ns. The

integrator, therefore, tracks only the solution for the species

mole numbers. The temperature was also computed when the

species time derivatives and the Jacobian matrix were evaluated.

In method B the temperature was treated as an additional

dependent variable and evaluated by solving its ODE (eq. (9)).

In this method, the number of ODE's is equal to N s + 1, the

Jacobian matrix is of size (N s + 1) × (Ns + 1), and the integrator tracks the solutions tbr both the species mole numbers

and the temperature.

The following naming convention was adopted. Techniques

using method A were given the suffix A (EPISODE-A, etc.),

and those using method B were given the suffix B (EPISODEB, etc.).

The code GCKP84 allows for heat transfer between the

reacting gas mixture and its surroundings and must therefore

use an ODE to solve for the temperature. It also includes the

density and velocity, V. of the gas mixture as dependent variables

and evaluates them by integrating their ODE's. (For the static

test problems used in this study the velocity ODE is given

trivially by dV/dt = O, V(to) = 0.) Consequently, the number

of ODE's solved by GCKP84 is equal to Ns + 3, and the

Jacobian matrix is of size (Ns + 3) × (Ns + 3).

CREKID computes the temperature by solving the algebraic

enthalpy conservation equation (8). However, the calculation

procedure is different from that used in method A. In CREK 1D

the mixed differential-algebraic system of equations (2) and

(8) is solved simultaneously, whereas method A solves

equation (8) after the species ODE's have been integrated

over a time step. Thus, although the number of ODE's solved

by CREK1D is equal to Ns, the Jacobian matrix is of size

(Ns+ 1) × (Ns+ 1).

Test Problems

The algorithms examined in the present study were applied

to the same two test problems used in our previous work

fief. 2). Both problems describe adiabatic, constant pressure,

transient, batch chemical reactions and include all three combustion regimes: induction, heat release, and equilibration.

Test problem 1describes the ignition and subsequent combustion of a mixture of 33 percent carbon monoxide and 67 percent

hydrogen with 100 percent theoretical air at an initial temperature

of 1000 K and a pressure of 10 atm. It comprises 12 reactions

which describe the temporal evolution of 11 reacting species

(CO, CO2, H, H2, H20, N, NO, N 2, O, OH, and 02). Test

problem 2 describes the ignition and subsequent combustion

of a stoichiometric hydrogen-air mixture at a pressure of 2 atm

and an initial temperature of 1500 K. It involves 30 reactions

among 15 species (Ar, CO 2. H, HO 2, H 2, H20, H202, N,

NO, NO 2, N2, N20, O, OH, and 02), of which two (Ar and

CO2) are inert. The reaction mechanisms and forward rate

coefficient parameters for the two test problems are given in

tables I and II.
TABLE L--REACTION MECHANISM AND FORWARD RATE

COEFFICIENT PARAMETERS USED FOR TEST PROBI.EM 1

IRate coefficient kI = lO&/rv,_ exp(-E/RT).]

Reaction

number,

.i

Reaction

I CO + OH = CO, + H

2 H + O, = O + OH

3 H, + O H + OH

4 H,O + O =OH +OH

5 H + H,O = H_ +OH

6 N + O, NO + O

7 N, + O = N + NO

8 NO + M = N + O + M

9 H + H + M =H, +M

I(1 O + O + M = O_ + M

II H + OH + M - H,O + M

12 H, + O, -OH +OH

Rate coefficient parameters

Bi NJ Et

kca],/Mole

11.49 0 0.596

14.34 16.492

13.48 9.339

13.92 18.121

14.(I " 19.870

9.81 1.0 6,250

13.85 0 75.506

20.60 - 1.5 149.1125

18.00 1.0 0

18, 14 - l.O 0.34(I

23.88 -2.6 (1

13 .(R) 0 43.([011

TABLE II.--REACT1ON MECHANISM AND FORWARD RATE

COEFFICIENT PARAMETERS USED FOR TEST PROBLEM 2

IRate coefficient k_ = 10_,T v, exp( -Et/RT).]

Reaction

ntlnlber,

J

Reaction Rate coefficient parameters

1 H + O, = OH + O

2 O + H_ = OH + H

3 H, + OH = H:O + H

4 OH + OH =O + H,O

5 H + O, + M = HO_ + M

6 O+O+ M =O, +M

7 H+H+M=H,+M

8 H +OH + M = H,O + M

9 H, + HOt = H,O + OH

10 H,O, + M = OH + OH + M

II H_ + O, = OH + OH

12 H + HO, = OH + OH

13 O + HO 2 = OH + O,

14 014 + 140, = H,O + O,

15 HO, + HO, = H,(), + O,

16 OH + H_O, = H_,O + HO,

17 O + H,O, = OH + HO,

18 H + H,O, = H,O + OH

19 HO, + NO = NO, + OH

20 O + NO, = NO + O,

21 NO +O + M = NO, + M

22 NO, + H = NO + OH

23 N + O, = NO + O

24 O + N, = NO + N

25 N + OH =NO+ H

26 N_O + M = N, + O + M

27 O + N.O = N, + (),

28 O + N,O = NO + NO

29 N + NO, = NO + NO

30 OH + N, = N,O + H

Bi Nj Ej,

kcal..'mole

14.342 0 16.790

10.255 1.0 8.9(X)

13.716 0 6.51X)

t2.799 1.093

15. 176 - 1 .O(X)

13.756 - 1.788

17.919 1.0 0

21.924 2.0 0

I 1.857 0 18.700

[ 7.068 45. 500

13.0(0) 43.0(X1

14.398 1.9(X)

13.699 I ,(XX)

13.699 l ./;_)

12.255 0

]3.0011 I. 8(/t)

i

13. 903 I. (XIO

14.505 9.(X)0

13.[)79 2. 380

13.(X10 0. 596

15.750 - [. 160

14.462 i I).795

9.806 1.0 6.250

14.255 0 76.250

13.602 0

14. 152 51.280

13.794 24.520

13.491 21.8(X)

12.556 0

12.505 , 80.280
10 0

10"1

10"2

10-5

10"6

3500

1500

10 -7 1000

10 -6 10 -5 10 -4 10 -3

Reaction time, t, s

Figure 2.--Variation with reaction time of temperature and species mole

fractions for test problem 1.

100

10-1

10-2

_- 10-3

t..-

10-4

10"5

10 -6

4000

35OO

3OOO

25O0

2000

10 -7 1500

10-6 10"5 10"4 10-3

Reaction time, t, s

g

O.

E

Figure 3.--Variation with reaction time of temperature and species mole

fractions for test problem 2.

Figures 2 and 3 present the variations with time of the

chemical species mole fractions and temperature for test

problems I and 2, respectively. These solutions were generated

with LSODE-B using a small value (10 -5) for the local

relative error tolerance. Both test problems were integrated

over a time interval of 1 ms in order to obtain nearequilibration of all chemical species and the temperature.

Computational Procedure

For each method, global errors in solutions generated with

a certain value for the local tolerance EPS were estimated by

comparing them with results obtained with the same method

and a reduced tolerance. The solutions used as a basis for

comparison were the most accurate generated and are referred

to as standard solutions. For example, for CREKID solutions

used as standards were generated with CREKID and

EPS = 10 -5. These standard solutions were used to estimate

the global errors in results produced with CREKID and

EPS = l0 _, 10 -3 and l0 4. The above procedure for

estimating global errors is reliable provided the technique is

effective in the sense that reducing the local tolerance actually

reduces the global error (ref. 24). In any case, in the absence

of exact solutions the only method for assessing the accuracy

of an algorithm is to compare the solutions that it produces

with those obtained with a reduced local tolerance using either

the same algorithm or a different one. The use of solutions

generated by each technique as a standard of comparison only

for itself ensures that the accuracy comparison is not biased

in favor of any one method or code.

A typical computational run was performed by first initializing

the time (t, set equal to zero), species mole numbers, and

temperature. The integrator was then called with values for

the necessary input parameters, including the local tolerance

and the elapsed time (equal to 1 ms for both problems) at which

the integration was to be terminated. After each step successfully

executed by the integrator, the current time and the solutions

for the species mole fractions and temperature were saved.

This procedure was repeated until the time reached by the

integrator was greater than or equal to 1 ms. The saved time
values served as input data for the output stations at which

the standard solution was to be generated. At each of these

discrete time values global errors in the species mole fractions

and temperature were estimated by comparisons with the

standard solutions as follows:

x,(t)

e,(t) - 1 i= I ..... N_ (24)

Xi.sr(t)

T(r)

er(t) -- 1 (25)

Ts7(r)

where e i(t) and eT(t) are, respectively, the estimated global

errors in the mole fraction xi(t) of species i and the temperature Tit) at time t and where XLST(t) and Tsr(t) are,

respectively, the standard solution values for the mole fraction

of species i and the temperature at time t. To prevent the

possibility of requiring accuracy in species with immeasurably

small concentrations, global errors were not measured for

species whose standard solution mole fractions were less than

O. 1 ppm. For such species e_(t) was set equal to zero. In this

way time histories of the global errors in species mole fractions

and temperature were generated.

For each technique, standard solutions were generated with

a small value for EPS. In addition to EPS and the elapsed time

at which the integration was to be terminated, other input

parameters were required by all codes examined. In this paper

only those input parameters that affect the accuracy of each

code are discussed. A more detailed discussion of these

parameters can be found in part I (ref. 2).

The user-supplied parameters relevant to solution accuracy

that are required by LSODE are the error control flag, ITOL,

which indicates the type of local error control to be performed,

and the local relative, RTOL, and absolute, ATOL, error

tolerances. Both RTOL and ATOL can be specified either as

(1) a scalar, so that the same local error tolerance is used for

all variables, or (2) an array, so that different values of the

local error tolerance are used for different variables. In the

present work the error control given by ITOL = 2 (for scalar

RTOL (equal to EPS) and array ATOL, see appendix A) was

used for reasons given below. Since the same number of

accurate significant figures is acceptable for all solution

components, RTOL was specified as a scalar. Now, for the

test problems examined in this study, the species mole fractions

and temperature vary widely (figs. 2 and 3), so relative error

control is appropriate and is the reason for designating the local

relative error tolerance as the primary tolerance (eq. (20)).

Pure relative error control can be achieved by specifying a

value of zero for the local absolute error tolerances. However,

since many of the species had zero initial concentrations, pure

relative error control could not be used. To make the error

control mostly relative, small values were specified for the

absolute error tolerances for the species mole numbers; for

convenience the same value (equal to ATOLSP) was used for

all species. Since the temperature can never be zero, pure

relative error control was used for this variable, that is, the

local absolute error tolerance for temperature was set equal

to zero. Thus, ATOL was specified as an array.

The values used for ATOLSP were those obtained in part

I of this study (ref. 2) for LSODE-B as follows: With

EPS = 10 -5, ATOLSP was progressively decreased until the

temperature-time trace showed essentially no change with a

further decrease. The values obtained for ATOLSP by using

this procedure were 10-14 and 10 -_ _, respectively, for test

problems 1 and 2. For consistency, the same EPS and ATOLSP

were used with LSODE-A, They were, however, checked to

ensure that reductions in ATOLSP resulted in essentially the

same solutions. For reasons given in the next section ERMAX

was set equal to EPS.

To make accuracy comparisons among the codes meaningful,

the same value of EPS (i.e., 10 -5) was used to generate

standard solutions for GCKP84, CHEMEQ-A, CHEMEQB, and CREKID. With EPISODE, however, an EPS value

of 10 -6 was used because larger values produced physically

meaningless results for test problem I--little or no change from

initial values after an elapsed time of 1 ms. The error control

to be performed by this code is selected by means of the flag

1ERROR (appendix A). For the reasons given above, pure

relative error control (option IERROR = 2) could not be used

and the option IERROR = 3 was used, instead. This error

control is semirelative (sec eqs. (17) and (18)). It is relative

for a variable that is initially nonzero. But for a variable that

is initially zero, it is absolute until the variable reaches unity

in magnitude, when it becomes relative. Since none of the mole

numbers attains a value of unity, the error control is always

absolute for species with zero initial mole numbers.

The solution generated with EPISODE depended on the

value specified for the initial step length (h 0) to be attempted

by the integrator. In generating standard solutions with this

code, ho was progressively decreased (with EPS = 10 -6)

until the temperature-time trace showed essentially no change

with a further decrease. The values obtained for ho by using

this procedure were 10-9 and 10 s s, respectively, for test

problems 1 and 2, for both EPISODE-A and EPISODE-B.

However, an h0 value of 10 -9 s was used tbr test problem

2 because it resulted in smaller execution times, as shown in

table III. For EPISODE-A the savings were modest, but for

EPISODE-B they were significant.

GCKP84 uses the same error control as that selected to be

performed by EPISODE. It also requires the user to specify

h0. Since details of the integration technique used in GCKP84

were not known, a default value of ho = 10 `6 s was used in

our previous work (refs. 2 to 5). However, Bittker and Scullin

(ref. 15) have since then set the default value for h0 at

5 × 10-s s. Nevertheless an ho value of 10 6 s was used in

this study to be consistent with part I (ref. 2). In addition, as

shown in the next section, the 10 -6 value generally produced

more accurate results than the new default value, while requiring

comparable execution times for all EPS used in this study.
TABLE III.--EFFECTS OF INITIAL STEP

LENGTH ON EXECUTION TIMES REQUIRED

BY EPISODE-A AND -B (EPS = 10-_')

FOR TEST PROBLEM 2

Method

1SODE-A

'ISODE-B

Initial step

length,

ho,

S

I0 8

10 _

10 -s

10-,J

CPU

execution

time.

S

3.1

3.0

14

7.8

In contrast to EPISODE and GCKP84, the other codes

automatically compute the h0 value to be attempted by the

integrator. In LSODE the calculation procedure for h0 employs

the user-specified values for the first output station and the

local error tolerances. The computed initial step length can

have an adverse effect on both the computational work and

the solution generated by the code (ref. 3). The calculation

procedures used for h0 in CHEMEQ and CREK 1D are based

on the problem physics (see appendix A) and the computed

ho did not cause the above difficulties.

Both CHEMEQ and CREKID use a relative convergence

criterion (eqs. (22) and (23)). The difficulty of applying the

test when the solution vanishes is avoided by setting mole

numbers less than a suitably small value, TINY, to be equal

to TINY. In this study a value for TINY of 10 -2o was used.

The only user-specified parameter required by CREK 1D that

affects its accuracy is AT, which is the maximum temperature

change allowed before the reaction rate coefficients and the

thermodynamic properties ¢,_and c_,._are updated. Use of this

parameter increases the efficiency of numerical techniques in

solving combustion kinetic rate equations (refs. 2 and 3). To

ensure that the most accurate solutions were used as standards,

a value of AT = 0 K was used for both test problems.

Results and Discussion

The procedure described in the previous section was used

to study the global errors incurred by the different techniques

in solving the two test problems. All results presented herein

were generated on the NASA Lewis Research Center's IBM

370/3033 computer using single-precision accuracy, except

GCKP84 which uses double-precision accuracy.

Both temperature calculation methods A and B were

attempted with EPISODE, LSODE, and CHEMEQ. The error

control used in method A is pure relative, and the local relative

error tolerance is equal to ERMAX (see appendix B), In

EP1SODE-B and CHEMEQ-B the error or convergence

control for the temperature is pure relative and the local relative

tolerance is equal to EPS (eqs. (17), (18), and (22)). For

reasons given previously in the section "Computational

Procedure" the above remarks apply to LSODE-B also, To

make accuracy comparisons between the two temperature

calculation methods meaningful, ERMAX was set equal to

EPS, thereby imposing the same local accuracy requirements

on both methods. Thus, both methods A and B used the same

error control (i.e., pure relative) and the same local tolerance.

To facilitate accuracy comparisons among the different

techniques, the species were divided into three types: reactants

(R), intermediates (I), and products (P). At each discrete time

at which global errors had been computed, root-mean-square

(rms) errors, erms,/(t) (j = R.I.P). were computed for all

three species types as follows:

erm_j(t) = e?,j(t) j = R. 1, P (26)

i=1

where ei,j(t) (i = 1..... Ni) is the estimated global error at

time t in the mole fraction of the i_h species of type j and

where Ni is the number of species of type j. The values of

N i (and the species that comprise each subset) are as follows:

For test problem I, Nn = 4 (CO, H2, N:, and O:), Nt = 4

(H, N, O, and OH), and Np = 3 (CO:, H_O and NO). For

test problem 2, NR = 5 (Ar, CO:, H e, N_, and O:), N/ = 6

(H, HO2, H:O:, N, O, and OH), and Np = 4 (H,O, NO,

NO2, and N:O). Although Ar and CO, are inert species, so

that their mole numbers do not change during the course of

the reaction, they are classified as reactant species because

they participate in three-body reactions as catalysts and their

concentrations affect the rates of these reactions.

In addition to the rms error for each species type, a single

rms error for all variables, erm_(t), was computed at each time

t by using

e_,._(t) = INs ll2

e;(t) + e_(t)

i=1

Ns+ l

(27)

Figures 4 to 9 present the variations with time of the percent

rms error in reactants, intermediates, and products and the

percent error in temperature for test problem 1. Similar inlbrmation is presented for test problem 2 in figures 10 to 17. For

brevity, test problems 1 and 2 are hereinafter referred to as

PI and P2, respectively. Note that for clarity the actual errors

have been magnified in some of the figures. The maximum

percent rms errors incurred and the reaction times at which

they occurred are given in tables IV and VI, along with lhe

values used for the input parameters discussed in the previt_us

section. For each code (except GCKP84) and EPS, these input

parameters, obtained in part I of this study (ref. 2) by a trialand-error procedure, minimized the execution time required

to solve the problem. To prevent the possibility of generating

physically meaningless results by using too large a value of

I0 Continued on p. 29
10

0

.

I.

tU

(a)

, I

L

1llllW

tolerance,

EPS

.......... 1x10 -2

................... lx10 _

1_o-'

l

"_..:__-fl'ltllt'Tt_"H_'{'t_

16

12

(b) A_

ll"i Illliti -_-_rlIi[IT;TT_:t't:H_B

14

(c)

12

10

8

4

2

0

lo4

,,

Reactk)n time, t, s

4 m

(d)

m

3

2

!

0

lo4

1IIIII :_,:t-.I-t-tl)ll. ...I...L.I_tlIIII

10...5 10 -4 1@ _

Figure 4.--Variation with time of percent rms error in (a) reactants. (b) intermediates, and (el producls, and of (d) percent error in temperature. Test problem

1 results generated using LSODE A.
10

8

0

G)

e,,

10

R

(,_

m

I

Local

tolerance,

EPS

......... 1x10-2

................... 1_0-3

1YI0-4

n

(c)

m

-

- t

) I I Illlll Ui..._l_II lJ]lll

16 m

(b)

m

m

B

12 m

$°°o.

o -1,I"IIllIIIl

m

(d)

t

o7 0

lO-S lO-a 10--4 10-3 104 104 10-4 10-3

Reaction time, t, s

Figure 5.--Variation with time of percent rms error in (a) reactants, (b) intermediates, and (c) products, and of (d) percent error in temperature. Test problem

I results generated using LSODE B.

12
800

50 Local

F(a) ,_ tolerance,

i-- i_ _P_ _o

,oF-- il ......... _o_

/ _1 ................... 1×1o-3

10_ ill _tt 200

o

II

45

35

25

(b)

-

-- II

- IZ -/ii

- /

-A_I III'

m

(d)

L t 1t 1_11111I I l!11111

15

J

5 m

IIIIIH ___ I

10 -3 0"6

Reaction time, t, s

l

_J_

!1IIIIII

lO--S

I

I ! lilllll

10 -4

I !1IIIIH

1u-3

Figure 6.--Variation with time of percent rms error in (a) reactants, (b) intermediates, and (c) products, and of (d'l percent error in temperature. Test problem

1 results generated using GCKP84.

13
50 ill

(a)

ii

40 -

m

30 -

20 "

i!

10 ,-

iii

0 I

0

_25 ,-

W

2oi

15 -

10--

5 D

10 4

toleranL°C%e, /

(error ×

3°I j

b)

25 .

20

t 10-32_

15 (error ×

10

s I- r'/ IL-lo -_

o _o1

(c)

r- 10 .4 /

, I I III '_L-'-

10_s 10-4 10 3

.5--

(d)

L_J

0 ..... _--Ii

_o5 mm

-1.0 -

-1.5 -

-2.0 - t I IliUm

lo-_ lo-s

Reaction time, t, s

10-3

10-4

(error x 1O)'''_

I I I Jllil! I I I IIIlii

_0-4 10-_

Figure 7.--Variation with time of percent rms error in (a) reactants, (b) intermediates, and (c) products, and of (d) percent error in temperature, Test problem

I results generated using CHEMEQ A.

14
"E

tolerance,

7O a) EPS

60 ---- 10-2 "--'/ 50 -- i10-3 j_ //_ 401" (err°rx2)--Jt"/

40 -- (error x 5)--. -.

30

30 ---

20

20 = L---IO -4

10 -- _" f (error × 251 10 xl0)

-I I lllllll- LgllIIIIL/ IIIIIIII o

0

o

n

3.0- U

25 (d) 1

c) 10-3 --_ .. 2.5 --

20 --

2.0 "

lo-= ----. / ,-lo_

15 -- 1.5 , //_,_# t

1.0 -

10 --

.5 D

5-- 10.4 __ .A.

rror x 10) 0 - V (error x 5)

iI zI Zjl,, IIIIIl[- --,-,-,ll,,--_ t-lVVlllllil jl i •I ,,,,,,, lllllll -o.5 i i Iilllll I I 1111111 I I IIIllll

0 10_ 10 4

10-4 10-3

1 _ 10"8 10-4 10-3

Reaction time, t, $

Figure 8.--Variation with time of percent rms error in (a) reactants, (b) intermediates, and (c) products, and of (d) percent error in temperature. Test problem

I results generated using CHEMEQ-B.

15
30

25

2O

15

10

5

0

8.

_3o--

tll

25

2O

15

10

(a)

i JIIIIHI

L_al

tolerance, /{_

EPS

........ 1><10-2

.............. 1x10-3

lx10-4

(c)

pj

_ /

I

/

-/

!

/ /'_ifl

o /__

/./'\

_, l _ ........

-t.-_.:l_!_.....!'"['!_HJ!

(b)

40

35

30 m"

1

25 _ /_

20 _ r

!

j"

I

15 _ /.

10 _ /*

7

{d)

6

5

4

3

I'"_ 2

/ \

/

/ \ 0

i / ,

S_m,w ° --1 m

.g_L_4t'lJll...,.l..._-f'tl_. m

10-4 10-3 10 4

Re_-'tion time, t, s

I I IIIII11

tO-a

L/"-,,.

f °e

I Ii111111I I I1_1111

10-4 10 _

Figure 9.--Variation with time of percent rms error in (a) reactants, (b) intermediates, and (c) products, and of (d) percent error in temperature. Test problem

I results generated using CREKID.

16
a.

uJ

25

20

15

10

m

Local 2OO i{b)

tok_rance,

EPS 150

........ 1x10 -4

lx10 _

100

50 .,,7"

0

_11.,111li111l_ I IIIII111

'001''i (c)

80

60

40 _ _i

20 _/

o ! IIIIIIII

0-7 10 4

4 m

(d)

J

3

2

1

0

l_.t.d3.11111111111111_, I

10 _ 10 -4 10 -3 10-7

Reaction time, t, s

J o,oO_o°o_..

IIIIIIII I Ii111111I IIil1111I IIIIII11

10-s lIP s 10 -4 10-3

Figure lO.--Variation with time of percent rms error in {a) reactants, (b) intermediates, and (c) products, and of (d) percent error in temperature. Test problem

2 results generated using EPISODE-A.

17
25 F(a) _ Local

20_.- !_ tO_E_Snce,

g.......

u_'120 F----(c)

100 -- +- /;"i

,og./!i...... ;....

o1. II1111111II1|1111 I_JI, I IIIIilU

IO-T 10-6 lO.-S 10-4 10 3

200 --

m

150 --

100

50

o _11lllllll

4 !

(d)

m

3 --

2 --

i

1 --

0

m

-1 + I

10 -.7

Ill ! Itl.llllH IIIIII+

Reaction time, t, s

t

IIIii!11IIIIIIIII I._11111111 IIIII!111

10-_ 10_ 10-4 10-3

Figure 11 .--Variation with time of percent rms error in (a) reactants, (b) intermediates, and (c) products, and of (d) percent error in temperature, Test problem

2 results generated using EPISODE-B,

18
I1

_l I_L,a_,..,_WIlIII IIIIII I It

C

M,I

25

(c) .'D, (d)

20 _ _I_ "_;'j I 1.0

15

_iI .5_

10 -- j_

! !

,, _)

o

I11111111 II .':.I.4Hl_

V

I IIIIIU I IIIIIIII I IIIIIIII

10-8 lO-4 Io--3

_._ I I IIIIIII I

0 10-s

10-7 10 -6 10 -s 10 -4 10 -3 10-7

Reaction time, t, s

Figure 12.--Variation with time of percent rms error in (a) reactants, (b'Jintermediates, and (c) products, and of (d) percent error in temperature. Test problem

2 results generated using LSODE-A.

19
108 I (a) En'or ×10"--_

i'_ ..... 1: fi .....

50 m

(b)

!

30_ !

!

I i !'"''""

,o-_ /

0

EPS

........ IxI0 4

............. I_0 4

t_O 4

,-- Error x 10

/

Error × 100--,,

%

25 E(c)

2O

J• I

10

5

,,./-/ I*°_ i

1

I

I

I

1.5

(o3

1.0

.5

A

(.

i Error x 10 -"_

I _ .,-"'"'_",

_......./ _,.k_,.. ol

,._.:..__

10 4 104 10 4 10-4

:,I,.UJ!IIII.....t-_tG 11IIIIIII

0

10 4 10 4 10 4 10_

Reaction time, t, s

Figure 13.--Variation with time of percent rms error in (a) reactants. (b) intermediates, and (c) products, and of (d) percent error in lemperalure. Test problem

2 results generated using LSODE-B.

2O
.o

Q.

I,LI

35

30

25

2O

m

(a)

15_

Local

tolerance,

EPS

........ 1_10- 2

................. lx10- 3

_____._. lx10- 4

r- Error x 2

oI II IIIIII._ll_lll_.J [[IIIIIL...IJ.II_

350

300

250

(b)

- 1

1

200 _ I

f

150_ _

!

100_ f :

;

i "

/

50 _ / z

0

r- Error x 2

/

pl

_1111.I [IllllIL.I_11111111

2OO

150

(c)

100

\:

5om j

l'r,

L,

0

10-7 10-6

Error x 2

/

/

/

r',

ll_JlillL._l I.LII_

10 -5 10 -4 10 -3

5 m

(d)

4_

3_

2

1

-1

J

--.1-

_ _.__i

I,

rError 2

',,.,_/ "'£,.,,,"'-'-

L I ll[lllll I[1111111I!111111t I IIIHW

0-7 10-6 1 0-s 10-4 10-3

Reaction time, t, s

Figure [4.--Variation with time of percent rms error in (a) reactants, (b) intermediates, and (c) products, and of (d] percent error in temperature. Test problem

2 results generated using GCKP84.

21
30 30 "

"" Local

(a) tolerance, . 25 . (b) 10-2--.,.,,,,__/ .... EPS

25 I0-2 "-_ ,jr 10.-3

"",/_ (error x 10)--_.. (!

20 ---" 10_ 20 -" _ ,, I I

(error x 10) --,, _%/ _J

15 _ 15 -

10 10- ,o-, 4 Ill

I I IIIIIIL JJAJIIIIBI Jl-# |1 IPn

0 0

.5 r- lo-, &

(C) 10"2"",, / (error x 501 ---- ,, /I

I (d) ..,. J,.

30 m %%%

(errorx 10)-_r -"

15 ----_'_ -1.0

10 -- 10 -4

,s

-1.5

5

0 -2.0

10-7 10-6 lO-a 10-4 10-3 10-7 10-6 10-6 10 _ 10-3

Reaction time, t. s

Figure 15.--Variation with time ofpercent rms error in (a) reactants, (b) intermediates, and (c) products, and of (d) percent error in temperature. Test problem

2 results generated using CHEMEQ-A.

22
Local

tolerance,

EPS

(a) to _ -,,

15 f(error x 5 /,i1

lOl-- lo-, f I

, ,__x_/.//

0.

W

25

20 m

15 -

10 -"

5 -

m

(b)

m

I I innnmm

,'o.%x,_./...j

<:._x,,_./.. f

,,,,,,.,,,/..2'J,

20 C(c)

10-2"--"" _

16 _ ",,

12 (error x

8

o -I IIIII!11 I I111111l lk_ll

10-7 10-6 10 4 10 -4 10-3

2.0

1.5

(d)

0 I I IIIIII!

10-7 10-6 lO-e

10-3

(error x 2)-_

%

(error 1 2)---

.I I I lille i i'*_ I IIII m

10-4 10-3

Reaction time, t, s

Figure 16.--Variation with time of percent rms error in (a) reactants, (b) inlermediales, and (c) products, and of (d) percent error in temperature. Test problem

2 results generated using CHEMEQ-B.

23
|

2.0

1.5

1.0

m

(a)

Local

tolerance,

EPS

........ 1_0-2

/

.............. lx10 -_ I

1 xl 0 _ //N,, /

/

ore--

m

(b) ,/_ m

/

s--- /-_ Ic)

'- 1

32- , !

0 i n,

mS m,am,

-1.o 1i1111111

10.-7 104

,_"';_...,,

_ I - _a,mmm,a,_ h

\............. .y % I

o Ii11,11111 I IIIIIIg I IIillli

10-7 104 104 10- 4 10_ 10_ 10-4 10 _

Reaction time, t, s

Figure 17.--Variation with time of percent rms error in (a) reactants, (b) intermediates, and (c) products, and of (d) percent error in temperature, Test problem

2 results generated using CREKID.

24
...4 ,.I ._

o

o .E

_a

_._

I--

..../

<

.= ¢/

E _

e_

_ _ _,=

N

I I I I I I I I I I I I I I I I I I I I q

_- & _,'__, _ _ _ ___'_', _-_ _ _°" _ _ ---

_6 _6 _ o6 .# -- _d ,5 # _. ei -- ,.H ,5 _d ,,,'.;-- ,'.i _; _ ,-.i

o o

III IIIIIII ILl III III

...... _ _ _ _ _ ---# _ _ _ r-i -- --

? ? ??? ??? _?? 7_7 777 7??

_ _ _ _ _ _,_ _-

i , , ,

_o00 00o _ •

XXX XXX

i

X X X X X X X X X X X X X X X X X X X X X

88= " "

x_

2

dl

25
i

i

J

u ]

='=-

o .

,-.,

I

_ _ _ ..-:c4 _ -.6 '.-: " _ _ .,_ ,.,.4_ _ _;',,4

I [ I [ I I I I d r I I ! I i i 1 1 d

_ _,----; _ _ _ ......

_ I 1 I I

• . ...;_ , ._-- . . .

I I 1 I I I I I I I I

_ "a 6' ,. o _ _'_.__

_ _ Z

_ -_ .... _ ..... o ......

I I I

??? _?iI 7TT? 7?7 ??? ??7

,, , i _FF_ _77 777 T_

'-r' '-

l- 0 _ Z _ 000

_>

II Ill III III1 III I11 Ill

....... _.. ,_ _ _--M

_' _ 6,_ _'_,--+_,_,

III III IIIII tll

X XXX XXX XXXX XXX XXX XXX

_ _ _ _ _ -

]

,4

26
e,,i

ua

,,-I

©

,..,,

.u.1 ,-.

k'-' 0

o

0 "_

U.,,I

,_.]

i'-'

8

,-,,,

°

2

II 111111 IIII III III III III

,,--.ir-i ,',4.,'_ _ ,'4,,-4.,'4 ,,'.,i_--_ _,"4,",4 --0404 ',_d, ....

I I I I I I I I I I I I I I I I I I I I 1 I I I I

m e_l eq

XXX XXX XXX

X XX x

0 000

XXXXXX XXXXXXX XXX XXX XXX XXX

,,', ,,', < ,_ d_ 6,

0 0 _ _, ,_ :_ :E _"

27
,.,3E

=

, _ 71_ _I'1' '_?_', 1_1_1, 77T 7?? "7_7

. -- _ _ I I I I q i I I i I

I I I I I ) I

,=

o ,

_= _ _ _ _ _. _-._

r_ _ ---- _ t-i ...... ,,'4 --; N ...... ._: .....

E,

??? ??? 777 _777 ?_ 777 777 ?7?

5,_,_,_,4,

_-_ ©_,_ _ _-_ _,_ _8 o_

El II lit} ]J rlll

XXX XX× XXX XX XX XXX XXX XXX XXX

£

I I

oo_

IP P_

28
ATOLSP, the runs with LSODE-A and LSODE-B were

requircd to satisfy the accuracy criteria described in reference

2. Several of the runs with the other codes also satisfied these

criteria. For LSODE-B and EPS = 10-4 two ATOLSP

values (10--s and 10 -_2) satisfied both the accuracy and

execution time criteria lot P2. Table VI includes results

obtained with both values, to illustrate the effect of ATOLSP

on the accuracy. However, the error plot given in figure 13

was generated with ATOLSP = 10 _-_

The maximum percent errors in each species type and

tempcraturc are given in tables V and VII for PI and P2. For

each species type the species incurring the maximum error,

the reaction times at which the maximum errors occurred, and

the standard solution values for species mole fractions and

temperature at these times arc listed.

For test problem l the runs with EPISODE-A and

EPISODE-B and EPS > 5x 10 -6 predicted little or no

change in the composition and temperature after an elapsed

time of 1 ms. Hence, maximum errors of - 100 percent were

obtained for both intermediates and products (table V).

Because the temperature and the more active reactants H2,

CO, and O_ display monotonic behavior (fig. 2) the

maximum errors in reactants and temperature occurred at

t_,,_, the final time (>_1 ms) at which the solution was

generated. Both EPISODE-A and EPISODE-B required only

six steps to complete the problem, and, because a new step

size is considered after every successful step, t_,_0 was

significantly greater than 1 ms. For intermediates and products

maximum rms errors of 100 percent were obtained at several

time values: hence, reaction times and standard solution values

are not given in tables IV and V tbr these two species types.

Also, no species name is listed in table V for either type

because all intermediate and product species incurred the

shown maximum errors.

The solution returned by EPISODE wax also tound to depend

on the output stations specified by the user. For example, for

some combinations of output times, EPISODE-B (with

EPS = 1× 10 6) predicted no change in the composition and

temperature after an elapsed time of 1 ms. However. by

stipulating only one output station (at 1 ms), the correct

solution was obtained. For this problem the run with GCKP84

and EPS = Ix10 -2 exhibited serious instability and was

therefore terminated. For reasons just discussed, no error plots

for EPISODE-A, EPISODE-B, and the run with GCKP84

and EPS = 10 -2 are presented.

Similar remarks apply to the results obtained for P2 with

EPISODE-A and EPS >__5× 10 -_, and with EPISODE-B

and EPS >_ 5× 10 -3. The run with EPISODE-A and

EPS = 5 x 10-4 required only seven steps to complete the

problem and t¢,d was therefore significantly greater than 1

ms, For this EPS, for exactly the same reasons given for PI,

the |bllowing quantities are not shown in tables VI and VII:

reaction times at which the maximum rms errors occurred in

the intermediates and products, intermediate and product

species incurring the maximum errors, the reaction times at

which the maximum errors occurred, and the standard solution

values.

For EPS = 5× 10 -4 and 10 3 EPISODE-B successfully

completed P2 in that correct solutions were returned at t = 1

ms. However, during heat release they were significantly

inaccurate. For example, the run with EPS = 5 × 10 4 predicted little change from the initial composition and temperature

until t -- 40 p,s when heat release began. In contrast, the

standard solution shows that heat release is almost over by

this time (fig. 3). As a result, maximum errors of - 100 pcrccnt

were observed for the products (table VII). This error was

incurred by several product species at several time steps;

hence, table VII does not list the product species name,

reaction time, and standard solution value. Because of the

difficulties experienced by EP1SODE-A and EPISODE-B,

error plots for EPS > 5 x 10 -4 are not presented.

As discussed previously, all results with GCKP84 wcrc

obtained with h0 = 10 -_' s, although its current default value

is 5 x 10 -8 s. The effects of this change in/tq_ on the accuracy

and execution time were studied by generating results with

h0=5×10 -_ s. The maximum errors incurred by the

solutions produced with both h0 values are given in tables

VIII to XI. For this study new standard solutions using

h_ = 5 × 10 -x s were established to bias the results in favor

of the current default value lbr h_. Despite this bias, tables

ViII to XI show that in almost all cases h, -- 10 _'s produced

more accurate solutions than h, = 5 × 10 -s s. (No results are

given for Pl and EPS = 10 "_because the runs with both h0

wcrc terminated due to instability.) For PI the results with

h0 = 10 -6 s were significantly more accurate for all values of

EPS (tables VIil and IX). Surprisingly, for h0 = 5× 10 x s

the solution with EPS = 10 4 incurred substantially greater

errors than those generated with the larger EPS. For P2 the

differences in errors obtained with the two h, were small fi)r

EPS= 10 ' and 10 3 but for EPS= 10 z, lt_= 10 -_ s

produced significantly more accurate results than

h, = 5×10 -s s (tables X and XI). Finally, the execution

times required with the two h0 are comparable fi)r both

problems and all EPS (tables VIII and X).

Examination of figures 4 to 17 shows sudden increases in

the error plots for intermediate species and products. This

behavior is caused by species reaching values of 0.1 ppm or

greater (figs. 2 and 3) and introducing their contributions to

the rms errors. For example, for PI the intermediate species

producing the sudden increases in the error plots are H (at

t = 2 _s), O (at t -_ 4 p,s), OH (at t - 4 p,s), and N (at t =

20/zs). For products the pertinent species are H_O (at t =

2 #s), CO2 (at t = 4 p,s), and NO (at t -- 15 p,s).

EPISODE-A, EPISODE-B, LSODE-A, LSODE-B,

GCKP84, and CREKID all experienced difficulty tracking the

standard solutions during induction and early heat release when

the species and temperature change rapidly (figs. 2 and 3).

The essentially isothermal induction pcriod ends and heat

release begins when the temperature starts to rapidly increase

from its initial value. During induction, the reactants and

29
3O

F--. zoXz<©

,-¢

e,,,

e¢

e",

t-.-. <

t"" Zz

,<a

e,, m<

[..., ©Z©

,.7- z

,-1

,-.I

_<

[.... 0

[-..

r_

.1

I>

,,..1

e_<

..

o

_'7 7';' -_i

(

XXXXXX

,,q o-: ,h., _. o.r-:

11 71' ?'f X X X X X

X

XXXXXX

I'-,I ¢xl --

-- r_

?? ?? ??

)'(XXXXX

r*"* oo O', _'_ irq t',q

_N _s xm' X X X X X X _ _ -- X X X X X X

o .-

_D

P_o

g

e_•

o

,,,,,,

..Eo

"_,_i X X

X

_Jz0

n_

n.i zoz

,.J<Ez

r.

ee_ <

o

&

XXXXXX

E_

g ..... i _ i i L

r

_-8 _;S $8

_ 6' ----_ o_' & _ .'.'£'

.... _,_'

77 i''7 "i'

_ xx ;x x_

.... ¢'-!

_.._

&,= && o& X X X X X X © O© T T T " T 3'

xx_

•, ._ ,_,_

_o_ 11 'f'?

gO _o o=

_._. xx xx xx

e¢

",_/_" ,,-i _;

IIIIII

,,_ ._m

_,_i_

×XX
Z

0

¢,J

X

U_

Z

<

o

U_

U_

<

_D

Z

Z

o

_z

z_

ge

Z

0

Z

.-.1

e_

p.

m

<

Z

rt

o

rn

fT.

I

M

M

X X X X X X

i_, ¢'q _ elm

X X X X X X

f'l ¢",1 -- --

X X X X X X

r,itq ,'q c.i _ _

'-- 0..i

,,,,,1

P_

++

t_

2

b

g

e,,,,

,..,/

,..2

c_

_ 2

1_t ?_ ?t

× X X × X X

-- f"3 e'e_ r'e_ 1"-1 ¢",1

X X X X X X

i-,_ ¢+-i r-I f-I _ _

X X X X X X

o

_+L "_

× X X

.~ _

w i

&g

•-=' N

_.=-_l x_ _ ×x

,_ _,. m-, , . .

i+++......

<_ _ ¢+

' I

+

=I • im

_-. + ,i.+

_'=

" I

_.._ ,_=+2 "-

E

oo +6 66

M ,*'i t'-.+ e; e_ P_

X X X X × X

,.-i ,--i r.i eq +,-i ,,-J

-7-

-- + _, m ,_ +,-+"Ii

_; x× x×

x x × x x N

x; ;x ;_

+_ ,.+:+ ,,_ ,,+

_,.+4 _ m.,;,:,6

I I I _ I I

,,+ ,f -+ .,+

x x x
32

<

,.,,,

F- zz<©

la_o

<ca

O, Ua

ZZ

<<

z..._

o_

z_

0"

,,g<

e_

..1

OZ

<

0<0

h-I M<

g

_._ _-_.

_'=

ILl

T_'T

_XXX

g=

"r,

,-,,,

7???

XXXX

_-_

XXXX

_-_

XXXX

74

XXXX

T?

8£

p

L.,

._o _.

_._=0

i-<

II

r,, Z<

_4 0

e,l 0

e_

FP_0

L_z0

e_

e_0Z0

r_X<

,..1

,Y0m<

_=N

u

_-

_44_

_xx

xxxx

_'= R:_

._ o___.

_._ &as

"_

_._ o°oo 7_7

_=

TI

g _TT7

_ ,T_

__< ....
temperature remain fairly constant. EPISODE, LSODE, and

CREKID have virtually no errors in the reactants and

temperature until heat release begins at -9 #s (PI) and 3 p,s

(P2), when the errors in these quantities start to increase. For

GCKP84, however, the error increases start at earlier times. The

difference, more noticeable for P1, is due to the smaller

reaction times obtained with this code for the onset of heat

release. Although, for consistency, we used EPS = 10 -5 to

generate standard solutions with GCKP84, values of EPS as

small as 10 -_ (for PI) and l0 -_ (for P2) were found to be

necessary to achieve tolerance independence of the

temperature-time trace at early times.

During induction, the intermediate species and the product

H20 increase rapidly from negligible initial concentrations.

The errors in the intermediates and products in this regime

are, therefore, relatively large. During early heat release

(t _< 15 #s for PI and t -< 6/_s for P2), these errors continue

to remain large as more products are formed and the

intermediate species continue to change quickly. In this regime,

the reactants show a sharp decrease, and the temperature rises

significantly. Many of the ODE's are unstable (ref. 2), and

so errors introduced at any step will grow as the integration

proceeds (refs. 1 and 24). For PI the reactants and temperature

vary rapidly between 9 and - 15 _s. For P2 the temperature

rise is not as steep, but the reactants change sharply between

3 and - 6 #s. Between these times, the errors in the reactants

and temperature are relatively large (figs. 4 to 7, 9 to 14+ and

17). For P2 the errors incurred at early times are less for

CREK 1D than for LSODE because of the much smaller step

lengths used by the former code in these regimes (refs. 2, 3,

and 12). During late heat release and equilibration, however,

EPISODE, LSODE+ and GCKP84 incur much smaller errors.

In these regimes the ODE's are stable (ref. 2), and so the errors

decay as the integration proceeds, provided, of course, that

the numerical method is stable.

The error plots for EPISODE, LSODE, and GCKP84

illustrate the dangers of assessing the accuracy of a technique

(or of a run with a certain value for EPS) by comparing

solutions at the final time (_ 1 ms for both problems). Note

that, although all these codes have negligible errors at the final

times, the errors can be significant at early times. For example,

with GCKP84 and EPS = 10 -3 the maximum rms error in

products is over 500 percent for P1. These plots also indicate

that if the main objective of the calculations is to study postheat

release phenomena (e.g., NO formation), the use of large error

tolerances does not result in significant errors. The large errors

incurred at early times, however, have important implications,

especially in developing and validating reaction mechanisms.

A procedure commonly used for this purpose is to compare

ignition delay times (e.g., time required for the temperature

to increase by a specified amount) predicted by the mechanism

with those measured in a shock tube (e.g., refs. 25 and 26).

The temperature error plots show that caution must be

exercised in using some of the codes to develop reaction

mechanisms by applying the above procedure. If, for example,

we assume that the ignition delay time is the time required tor

a 25 K rise in the temperature, values of - 11 and 3.5/_s are

obtained for P1 and P2, respectively. At these times, the error

in temperature ranges from l0 to 25 K for EPISODE, 2 to

5 K for LSODE, 15 to 200 K for GCKP84, and 0 to 10 K

for CREK 1D.

In contrast to EPISODE, LSODE, GCKP84, and CREK I D,

CHEMEQ incurs virtually no errors during induction and early

heat release (figs. 7, 8, 15, and 16). Therefore, this code can

be used to generate accurate ignition delay times. CHEMEQ

is superior in these regimes because of the very small step

lengths that it selects (refs. 2 to 5). However, as pointed out

by Young and Boris (ref. 10), the continued use of the hybrid

method used in CHEMEQ results in the global errors

increasing with time. For example, with CHEMEQ-A and

EPS = 10 -2, the rms error in reactants has risen to almost

50 percent for P1 (fig. 7) and 25 percent for P2 (fig. 15). The

situation is worse with CHEMEQ-B (figs. 8 and 16). During

equilibration, for CREKID, also, the errors grow (figs. 9 and

17) because the formulation used by it in this regime is based

on that used in CHEMEQ. However, CREK 1D incurs smaller

errors than CHEMEQ for most of the species types and for

the temperature. For EPS _ 10 -3 both CHEMEQ and

CREK 1D either are more accurate than or compare favorably

with LSODE during late heat release and equilibration.

Figures 4 to 17 and tables IV to VII show the large variations

in the maximum errors for the different techniques. EPISODE

and GCKP84 experience the greatest difficulty tracking the

solutions at early times--rms errors in excess of 100 percent

are obtained with the two codes. In contrast, the errors incurred

by LSODE, CHEMEQ+ and CREKID are significantly less.

Comparisons of the runs with the largest EPS value show that

LSODE is the most accurate code for P1, and CREK1D for

P2. Comparing the errors in the different regimes shows that

CHEMEQ is the most accurate code during induction and early

heat release. During late heat release and equilibration,

however, the other codes are more accurate.

Examination of figures 4 to 17 and tables IV to VII shows

that all techniques are tolerance effective in the sense that a

decrease in the local tolerance generally results in decreased

global errors. We note, however, that with LSODE not all

plots show an error decrease with EPS. On the contrary, for

some runs the error increases with a reduction in EPS (figs.

4, 5, 12, and 13). This behavior can be explained by examining

the nature of the error control performed in LSODE. As

discussed in the section "Computational Procedure," the error

control selected to be performed by LSODE is mixed relative

and absolute for species mole numbers and pure relative for

the temperature. For pure relative error control, the estimated

local truncation error, di, in species i approximately satisfies

the inequality

di <- EPS ,oi (28)

For pure absolute error control, d, approximately satisfies

33
di -< ATOLSP (29)

These two inequalities are approximate because the code

controls only the rms norm of the estimated local truncation

errors in all variables and not the estimated local truncation

error in each variable.

Equations (28) and (29) show that, since oi << 1, relative

error control is more accurate for a given value of the local

error tolerance. Hence, relative error control is appropriate

for the two test problems. However, when oi = 0, relative

error control cannot be used. This problem is resolved by using

a mixed relative and absolute error control, and

di< EPS oi + ATOLSP (30)

Equation (30) shows that for the error control to be relative,

ATOLSP must satisfy the inequality

ATOLSP << EPS :oj_ (31)

If ATOLSP >> EPS Ioi the error control is absolute. In this

study, we have considered only species with mole fractions

(xi) >- 0.1 ppm. This value ofx i corresponds to oi = 3x 10 -9

and 4x 10 -9, respectively, for P1 and P2. Hence, for the error

control to be always relative, ATOLSP must be less than

3x 10 9 EPS and 4x 10 -9 EPS, respectively, for PI and P2.

Only the runs with EPS = 10 -2 for PI satisfy these requirements. Hence, they are the most accurate at early times when

the mole numbers of many intermediate and product species

have very small values. Note that for P2 even the standard

solutions do not satisfy the requirement on ATOLSP. To ensure

their accuracy, the standard solutions generated with LSODEA and B were checked, respectively, against the solutions

obtained with LSODE-A and B using EPS = 10 -5 and

ATOLSP = 10-i.s, which _tisfy equation (31). These comparisons showed agreement to three significant figures for all

species with mole fractions > 0.1 ppm. For LSODE-A agreement to three significant figures was obtained for all species,

even those with mole fractions significantly smaller than

0.1 ppm. But tor LSODE-B the agreement for mole fractions

< 0.1 ppm was not good for NO. For all other species,

however, good agreement was obtained for mole fractions

> 10-ix. This observation indicates that LSODE-A is more

accurate than LSODE-B.

Equation (30) also shows that for given values of EPS and

ATOLSP, oi must satisfy the following inequality

ATOLSP

oi >> - Crmi, (32)

EPS

to achieve relative error control. As oi increases from zero,

the error control becomes less absolute and is equally relative

and absolute at O'mi n. For o i _ Omin, the error control becomes

increasingly relative as oi increases. Hence, the quantity Omin

may be regarded as the value at which, for increasing oi, the

error control starts to change character from being more

absolute to becoming more relative. Because .ri = o,/o,,,, the

value of x i (=xmin) corresponding to o,nin is given by

xnli, = ATOLSP/(EPS o,,,) (33)

For P1, x.,,in --- 3×10 -7 and 3×10 -6 for EPS = 10 -3 and

10 -4 and the ATOLSP given in table IV. These values are

attained by most of the species at t = 5 and 7 txs, respectively

(fig. 2). Hence, until these times the solutions with EPS = 10-3

and 10 4 are expected to be worse than or, at best, as accurate

as the run with EPS = 10 2. Examination of figures 4 and

5 shows that the errors in intermediates and products for

EPS = 10-3 and 10 -4 are worse than those for EPS = 10 "_

until t = 6 and 5/zs, respectively, for LSODE-A, and until

t --- 7 and 6/xs, respectively, for LSODE-B. In addition, all

maximum rms and maximum errors, almost all of which occur

at t > 7tzs, exhibit reductions with decreasing EPS (tables IV

and V).

For P2 and LSODE-A, xmin has values of 2.5x10 -5,

2.5x 10 -4 , and 2.5x10 -3, respectively, for EPS = 10 2

10 3, and 10-4. Some of the species never reach these values

(fig. 3). Hence, the errors do not show much sensitivity to

changes in EPS (fig. 12 and tables VI and VII). For LSODEB, however, the values used for ATOLSP ensure comparable

levels of relative error control for EPS = 10-2 and 10 3: for

EPS = 10 4, the control is more relative in the sense that it

has a smaller value of ATOLSP/EPS. The errors, therefore,

display decreases with reductions in EPS (fig. 13 and tables

VI and VIl). The sudden increases in the product errors around

t = 10/zs were caused by the species NO, which LSODE-B

had difficulty tracking (table VII).

The above discussion should be regarded as strictly approximate because it applies only to the estimated local truncation

errors, whereas figures 4 and 5 give the estimated global errors,

which represent the cumulative effects of the local errors. The

number of integration steps required up to the relevant reaction

time should therefore also be considered. However, the global

errors accumulate in a complicated manner from the local errors.

Other factors that must be taken into account are that LSODE

controls only the norm of the estimated local errors and that

different species reach mole fraction values of 0.1 ppm at

different times. Finally, although we have ignored species with

x i < 0.1 ppm, they do incur errors whose magnitudes are

controlled by ATOLSP and which grow with reaction time in

the initial combustion regimes, for reasons previously given. It

is therefore difficult to draw definitive conclusions about the

ATOLSP to EMAX ratios required for combustion kinetics

problems. For example, for P1, EPS = 10 .3 is expected to

produce more accurate results than EPS = l0 4 for the intermediates and products at early times, especially between 5 and

7 _s (see eq. (33) and the discussion following it), but figures

4 and 5 show the opposite behavior for both LSODE-A and

LSODE-B. One conclusion that can, however, be made is that

care must be exercised in specifying ATOLSP.

34
Theeffectof ATOLSPon solutionaccuracyis further

illustratedforP2andLSODE-Bbytheresultspresented intables

VI andVII (EPS= 10-4) andXII andXIII(EPS= 10-5).

Notethesignificant errorreductionsobtainedbydecreasing

ATOLSP.ComparingtheerrorsgivenforLSODE-Bintables

VIandVIIwiththoseintablesXI1andXIII,respectively, shows

thattbr thesamevalue(=I0"-_or l0 9)of ATOLSP,

EPS= 10-5 doesnotproducesignificantlymareaccurate

solutionsthanthelargerEPS.

TablesXII andXIII shawthat,althoughtheuseof large

values of ATOLSP can result in significant errors for the

intermediate species and products, the effect on the temperature is small. Therefore, if the user is interested only in

temperature versus time traces at early times, as for example

in developing reaction mechanisms from ignition delay times,

fairly large ATOLSP values can be assigned.

The results obtained above indicate that the ATOLSP needed

to achieve acceptable accuracy depends as much on the nature

of the solution as on the value specified for EPS and the minimum mole fraction to be considered in the error analysis. The

estimate for ATOLSP given by equation (31) may not be small

enough, as for example, PI and EPS = 10-2 (table IV). On

the other hand, the estimate may be needlessly conservative.

For example, although the intermediate species increase much

more rapidly for P2 than for P1, larger ATOLSP produced

results that satisfied the accuracy criteria. Because the value

needed for ATOLSP is a function of the problem, it can be

obtained only after the problem is solved. The major problem

associated with using LSODE to solve chemical kinetic rate

equations is therefore the trial-and-error procedure necessary

to obtain the optimal value of ATOLSP, that is_ the value that

minimizes the CPU time while satisfying prescribed accuracy

requirements. Note that for P2, although the runs with LSODE-B

and ATOLSP = 10 s and I0-12 (EPS = 10 -4 ) required the

same CPU time, the latter is significantly more accurate

(table VI). In contrast, the runs with ATOLSP = 10 -9, 10 -_°,

and 10 -I1 required about 2.7, 1.7, and 1.7 s of CPU time,

respectively. The trial-and-error search for the optimal ATOLSP

can be time consuming, especially for large systems of ODE's.

The use of an extremely small ATOLSP to ensure solution

reliability can result in excessive CPU times. For example,

for P2 the run using LSODE-B with EPS = 10 -5 and

ATOLSP = 10 -II required about 3.4 s of CPU time; in

contrast, the run with ATOLSP = 10-15 required almost 20 s,

although the solution was not significantly different.

The error control used in EPISODE and GCKP84 is pure

relative for species with initially nonzero mole numbers and

for the temperature; it is, however, pure absolute for species

with initially zero mole numbers. Since most of the species

have zero initial mole numbers for both test problems, the error

control is mostly absolute. Hence, for the same value of EPS,

EPISODE and GCKP84 are not as accurate as LSODE for

solving chemical kinetic rate equations. To achieve comparable

accuracy, especially at early times when oi is very small,

small values have to be used lor EPS (ref. 2). The runs with

EPISODE and GCKP84 were therefore more expensive than

the ones with LSODE (refs, 2 and 3). Modifying EPISODE

to employ the same error control as LSODE produced

significant reductions in execution times. Preliminary results

with the revised EPISODE indicate that it is as fast as LSODE.

For example, for PI the runs with the modified EPISODE-A

(EPS-- I0 2 10-3, and 10 4 and the ATOLSP given in

table IV) required, respectively, 0.35, 0.41, and 0.61 s. The

execution times compare very favorably with those required

by LSODE-A: 0.37, 0.46, and 0.63 s (refs. 3 and 4). The

above observations indicate that the error control used in

EPISODE and GCKP84 is inappropriate far combustion

kinetics problems.

Examination of tables 1V to VII shows that temperature

calculation method A does not necessarily produce less accurate

solutions than method B. On the contrary, far most of the runs

method A is more accurate for all codes; this result is most

apparent for CHEMEQ and PI.

To provide a more comprehensive measure for comparing

the accuracy of the methods examined, we adopted the

following procedure: For each run, a mean integrated rms

error, _ ...... was defined as

=- 1 _t_ gCTld Gn,,(t) dt (34)

_rms lend 0

where ten d (- I ms) is the end of the integration time interval

and e_,,,_(t) is given by equation (27).

Equation (34) provides a single quantity that is a measure

of the average error incurred in solving the complete problem.

The integral in this equation was evaluated numerically using

Simpson's rule (e.g., ref. 23), modified for unequal step sizes.

For runs requiring an odd number of integration steps, the

trapezoidal rule (ref. 23) was used on the last two mesh points.

The effects of h0 on £_,_ for GCKP84 are given in tables

VIII and X for P1 and P2. Except for tbe run with

EPS = 10 -_ for P2, h0 = I0 -6 s incurred either comparable

or significantly smaller _;rm_ than ho = 5x10 -x s. Note,

further, for PI and ho = 5x 10 -s s the substantial increase

in E,N,_ when EPS is decreased to 10 -4 (table VIII).

The variation of _;r.... with ATOLSP is given in table XI!

for P2 using LSODE-B and EPS = 10-5. This table illustrates

the increasing accuracy obtained by reducing ATOLSP. It also

shows that ATOLSP must be chosen carefully, as discussed

previously.

The variations of Erm_ with the user-specified kx:al tolerance,

EPS, are shown in figures 18 and 19 for Pl and P2,

respectively. We have included the run with EPISODE-B and

EPS = 5× 10 _ in figure 19 because it was successfully completed. The E..... given in figure 19 for LSODE-B and

EPS = 10 -4 was that obtained with ATOLSP = 10 J2.

These figures show that all methods are tolerance effcctive

(i.e., decreasing EPS results in reduced _;rn,_)" For both test

problems temperature calculation method A is as accurate as

35
10 0

_.- I"1 LSODE-A

• LSODE-B

X _KP_ CHEMEQ-A

• CHEMEQ-B

10 "1 CREK1 D

70

10.2

10 .3

m

10 4 illl l I

lO_

! I lilili I I 1 lilili I i I I

10 .3 10 .2 10 1

Local tolerance, EPS

Figure lB.--Variation of the nlean integrated rn]s error with the hx:al tolerance

for lest problem I.

method B. In many cases it is significantly more accurate,

especially with CHEMEQ and EPISODE. For P2 and

EPS = 10 -4, LSODE-B is more accurate than LSODE-A

because it used a smaller ATOLSP (table VI).

For the same value of EPS, EPISODE and GCKP84 are

significantly less accurate than LSODE (figs. 18 and 19)

because the error control used in the two codes is inappropriate

for chemical kinetics rate equations. For all techniques, note

the significant discrepancies between the values specified for

the user-supplied local tolerance and the errors actually

incurred. With CHEMEQ-B, a value of EPS = l0 -2 (1 percent) has resulted in an average error of almost 50 percent.

The relatively large E,,1,_ incurred by CREKID and CHEMEQ

is due to the difficulties which these codes experienced tracking

the standard solutions during late heat release and equilibration.

With LSODE, especially the runs with EPS = 10 -2, the

correspondence between EPS and {;r,,,s is better. These plots

show that for a given value of EPS, LSODE is the most

accurate code currently available for solving chemical kinetic

rate equations. However, for P2, especially with the smallest

EPS examined, GCKP84, CHEMEQ-A, and CREK1D

compare favorably with LSODE (fig. 19).

Figures 20 and 21 present the variations of the computational

work (expressed as the CPU time in seconds) with the mean

integrated rms error for problems l and 2, respectively. Note

the large differences in the CPU time required by the different

codes to achieve comparable accuracy. For P1 and a I/2 percent

mean integrated global error, the CPU time varies from about

0.4 s for LSODE-A to over 40 s for CHEMEQ-A. In general,

to produce an order of magnitude reduction in /_r,,l_ approx100 F- 4_ EPISODE-A imately doubles the computational cost. For both test problems

l- O' EPISODE-B LSODE is the most efficient code in the sense that it requires

[.- I'1 LSODE-A the least CPU time to attain a specified accuracy level.

[- • LSODE-B Figures 20 and 21 show that the CPU times required by

I'- V GCKP84 A

10 -1 L A CHEMEQ-A j_ _ temperature calculation method A are less than, or compare • c.E.E . / .,,I/

..p.o 10 _

I 101 V GCKP84 in A CHEMEQ-A

• CHEMEQ-B

i 0 CREK1D

il i I i 10 -1

10-4

104 I Illili I t I Itltlt I i I I

10 -5 10 -4 10 -3 10 -2

Local tolerance, EPS

Figure 19.--Variation (/1"the mean integrated rms error with the I(x:al tolerance

t_)r test problem 2.

I i l ilihl I i l ilihl I i I Illhl I

lO-3 lO-2 lo-1

Mean integrated rms error, %rms

Figure 20.--Variation ot the CPU time with the mean integrated rnls errt_r

tor test problenl I.

I I llihl

100

36
10 2 =-

101 _.

d

E

a.

O

100 --

10 "1

10-4

EPISODE-A V GCKP84

• EPISODE-B A CHEMEQ-A

[] LSODE-A • CHEMEQ-B

• LSODE-B O CREK1D

u l uldt[ I i l lldnl 1 i l lllhl I i I nlnhl

10-3 10-2 10"1 100

Mean integrated rms error, _rms

Figure 21. Variation ¢_t"Ihc CPU time ',vilh the mean inlcgraled rms error

for lest problmn 2.

favorably with, those required by method B. This difference

is most pronounced lor CHEMEQ and EPISODE. For

example, for P2 and a 1 percent ,g...... CHEMEQ-A required

only about half as much CPU time as CHEMEQ-B (fig. 21).

Note that for EPISODE B the computational work increases

with increasing error.

EPISODE-A compares very faw_rably with LSODE for P2

(fig. 21). However, the solution generated by EPISODE can

be strongly dependent on the value selected by the user for

h0 and a poor guess can result in incorrect and unstable

solutions (refs. 2 to 5). It can also result in excessive CPU

times. For example, the run using EPISODE-A with

EPS = 10-a and ho = 10 -s s required about 129 s for P2; in

contrast, the run with ho = 10-7 s required only 0.59 s.

Conclusions

The accuracy of several codes (EPISODE, LSODE,

GCKP84, CHEMEQ, and CREKID) in solving combustion

kinetic rate equations has been examined in detail. The

accuracy studies were made by applying the codes to two

practical combustion kinetics problems. _th problems dcscrit_,.d

adiabatic, homogeneous, gas-phase chemical reactions and

included all three combustion regimes: induction, heat release,

and equilibration.

During induction and early heat release, when the species

mole numbers and temperature change rapidly, EPISODE.

LSODE, GCKP84, and CREKID had difficulty tracking the

solutions. The errors incurred by EPISODE and GCKP84 in

these regimes were significantly larger than those incurred by

LSODE and CREKID. In contrast, the solutions generated

with CHEMEQ displayed virtually no crrors during induction

and early heat release. However, during late heat release and

equilibration, the errors obtained with CHEMEQ increased

significantly. In these regimes, the other codes were more

accurate.

Among the codes examined, LSODE was the most accurate

for solving chemical kinetics problems. This study has also

shown that LSODE is the most efficient code, that is, it

required the least execution time to attain a specified accuracy.

The major difficulty associated with its use is the trial-anderror procedure necessary to obtain optimal values for the local

absolute error tolerances R_r the variables. A poor guess for

the absolute error tolerance can result in excessive execution

times or in seriously inaccurate solutions.

An important conclusion is that calculation of the temperature

by solving the algebraic enthalpy conservation equation can

be more accurate and efficient than integrating its differential

equation.

37
Appendix A

Description of Codes Studied

The ordinary differential equations (ODE's) (2), (9) and (12) Ka

describing homogeneous gas phase chemical reactions can be Y_," : E OLJ'nYi'n-J "[- h,,

generalized as follows: j= 1

dyi

3;,-=-_-t =f((yk)) i,k = 1..... N

yi(t o) = given

(AI)

where for temperature calculation method A (see the section

"Evaluation of Temperature")

Yi = °i i : 1..... N s

N=Ns

for temperature calculation method B

(A2)

Yi = oi i= 1..... N S )

YNs+ l = T

N=Ns+ I

(A3)

and for the code GCKP84

Yi = oi i = 1..... Ns )

YNs+I := V YNs+ 2 =p YNs+3 = T _ (A4) N=Ns+ 3

In vector notation equation (A1) becomes

dy

) ---:_tt =f(Y) y(to) = given (A5)

where the underscore is used to denote a vector quantity. A

matrix is denoted by a boldface letter. This notation is used

throughout this appendix. In equation (A5) the N-dimensional

column vectors y and f contain the dependent variables and

their temporal derivatives, respectively.

The initial-value problem is to determine values for [Yi] at

one or more times in a prescribed integration interval, given

_/] and the values [Yi(to)l at the initial time to. We now

describe the codes studied in the present work and how they

solve the above problem.

EPISODE AND LSODE

Both these codes use linear multistep methods of the form

(refs. 6 to 9)

g 2

E

j=0

i= 1..... N

(A6)

where Y,.,, is an approximation to the exact solution yi(t,,),

f,,, (equal to f([Yk.,,l)) is an approximation to the exact

derivative );i(t,) (equal tof(Lvk (t,)/)), and the/%,,,] and l_j.,l

(30,,, > 0) are associated with the particular formula selected

by the user. The options include a variable-step, variable-order

implicit Adams method (suitable for nonstiff problems) of

orders 1 to 12, and a variable-step, variable-order backward

differentiation formula (BDF) method (suitable for stiff

problems) of orders 1 to 5. As discussed in the section

"Methods and Codes Examined," the BDF method was more

efficient for the problems examined in this study. Therefore,

the discussion is restricted to this method. For a BDF method

of order q, Kl = q, 1(2 = 0, and equation (A6) reduces to

q

Yi,. = E c_J,"Yi,"-J + hd3°,'_," i = 1..... N (A7)

j=l

The step length hn can vary from one step to the next in

EPISODE but is held constant for q + 1 consecutive successful

steps in LSODE. Hence, for EPISODE, [oq.nl and [3j.,,I can

vary from one step to the next, but in LSODE they are

predetermined constants corresponding to the order used.

Both codes use a predictor-corrector process to solve for

Y,. An explicit method generates a predicted value, _01

which is then corrected by iterating equation (A7) to convergence,

that is, the improved estimates yJ,ml (m = 1..... M) are produced

until the magnitude of the difference (_,,1_ _,,-tl), in

EPISODE, or (h,,_,, ''1 -h,_ m- Jl), in LSODE, approaches

zero within a specified accuracy. Here, _yl,,,I and __nml are,

respectively, the approximations generated for Y, and _, on

the m th iteration, the integer M is the number of iterations

required for convergence, and _ml is accepted as the numerical

solution at t_, provided it satisfies a prescribed local accuracy

requirement. At each iteration m, h,,_,,"1 is computed in LSODE

from y},ml via the relation

q

= /_ _ l)'[m ]

+....0.,,_.o

j=l

(A8)

so that the pair (_ylml, hn'_ytnml)satisfies the BDF method (eq.

(A7)) exactly. The predicted values of Y, and h,_,, denoted

by h,,_, °1, also satisfy equation (A8).

The predicted quantities EI,,°l and h,,_Yl,, °l are obtained by a

qth-order Taylor series expansion as follows: The history of

the solution is maintained in the Nordsieck array (which is

38
a Taylor series array) z,, of size N x (q + 1) (e.g., ref. 1).

The ith row (i = l ..... N) contains the q + 1 elements Yi.,,,

h,,_ ,,, h_,/2! Y,, ...... h[[/q: _'/,), where ytJ, is the approximation

to (dJY,./dt:),. The (q + 1) columns of z, are numbered from

0 to q, and the .fib column (/' = 0,1 ..... q), which will be

denoted by the vector z,, (j), contains the vector h,,Y_,,?j! of the

jCh-order scaled derivatives. If z, i has been obtained, the

predicted history matrix, z,I,°1, at t,, is given by (ref. 1)

z),°l = z,,_ iA(q) (A9)

where A(q) is a (q+ 1) x (q+ 1) matrix, with element Ai_(q)

given by

A_k(q) = I O(Jk) j>_k j<k 1 j,k = 0,1 ..... q

The binomial coefficient, (_,), is defined as

(_) - j.'

k' U - k).'

Thus, the predicted array z},°l is obtained by a simple qlh

order Taylor series expansion by using equation (A9). The

matrix z,I,°l contains predicted values of Y,, and its scaled

derivatives up to order q, the current method order. Note,

however, that because a qth-order Taylor series expansion

method is used z,l,°l(q) = z,,_ l(q).

The estimates y,/,,,l and, in LSODE, h,,Y,l'''1 (m = l ..... M)

are generated, as described below, until the iteration

converges. The local error test is then applied and, if passed,

the Nordsieck history matrix z,, is constructed by using the

relation

z,, = z},°l + e,,l,,(q) (A10)

where

e,, = El,,'vtl- 1A,,, °l

in EPISODE, and

e,, = h,,#,,, Mjin LSODE. The (q + l)th-dimensional vector l,,(q)

I,,(q) = (Io.,,(q), It.,,(q) ..... Iq.,,(q)) (All)

contains the method coefficients for the Nordsieck history

formulation of the q_h-order BDF method. Because EPISODE

and LSODE use different calculation procedures, the [li.,/_

values are, in general, different in the two codes. For

EPISODE, l,,(q) depends on the method order and the step

length history, satisfies Io.,(q) = 1 and/i. , = 1/t30 .... and has

to be recomputed at the start of each step. For LSODE, I,,(q)

is a function of only q, satisfies Io.,,(q) =/3o., and Ii.,,(q) = 1,

and has to be recomputed only when the method order is

changed.

To correct the initial estimate _,01 (i.e., to solve equations

(A7)), both codes include a variety of iteration techniques.

For combustion kinetics problems the most efficient is the

Newton-Raphson iteration (ref. 2), which is given by' the

recursive relation

q

p(yl,,,+,l_ 1A,,,,I)= E %.,,E,-, + h,,_o.,j (_;,,I) _ IA,,,,,I

.i= 1

(AI2)

for m=0,1 ..... M-I. The NxN iteration matrix P is

given by

P = ! - h,,_0.,,J (AI3)

where I is the identity matrix and J is the Jacobian matrix,

with clement J,j given by

J0 = Of,/O3) i,j = 1..... N

For this method, much computation time is required to form

the Jacobian matrix and to perform the linear algebra necessary

to solve equation (AI2). To reduce this computational work,

P is not updated at every iteration. For further savings, it is

updated only when it has been determined to be absolutely

necessary for convergence. Hence, the iteration matrix is only

accurate enough for the iteration to converge, and the codes

may use the same matrix over several steps of the integration.

In any case, both EPISODE and LSODE update P at least

every 20th step. The linear algebra required to solve equation

(AI2) is performed by the LU method (e.g., ref. 27), rather

than by explicitly inverting the matrix, which requires

prohibitive amounts of computer time (ref. 23).

Convergence of the estimates is ascertained as discussed

below. EPISODE constructs a vector Y,,,ax, which depends

on the user-specified value for the local error control IERROR

as follows:

IERROR = 1 (absolute error control):

Ymaxi= 1 i = I ..... N

IERROR = 2 (pure relative error control):

Y,,,,,x, = lYe,,, I I i= 1..... N

3q
IERROR = 3 (semirelative error control):

Ymax, = max[ IY/.,,_t I,IY/.,-2[I for i that satisfy Yi(to) ;_ 0

= max[i, IYi.,,_l II for i that satisfy yi(to) = 0

(A14)

The test for iteration convergence is based on the successive

differences (_,,,I _ _,,-II), as compared with _Ym_xand the

user-supplied local error tolerance parameter EPS. Convergence is said to occur if

N /'El,,0 _ y Im-II

6,,, -- Yn'7 -< GEPS (AI5)

In LSODE an error weight vector EWT is constructed as

follows:

EWTi = RTOLi IYi.,,_ I I + ATOLi i = 1..... N

where RTOLi and ATOLi are, respectively, the user-supplied

local relative and local absolute error tolerances for the i th

component. Both RTOL and ATOL can be specified either

as a scalar or an array, as discussed in the section "Computational

Procedure." The value of the user-supplied parameter ITOL

indicates whether RTOL and ATOL are scalars or arrays.

ITOL has four possible values which correspond to the types

of RTOL and ATOL as follows:

ITOL = 1: scalar RTOL and scalar ATOL

ITOL = 2: scalar RTOL and array ATOL

ITOL = 3: array RTOL and scalar ATOL

ITOL = 4: array RTOL and array ATOL

The convergence test is based on the successive differences

(h,,_,, mt- h,,_, m-II) as compared with EWT, and is given by

i=1

,\t/2

/h y.[m]_ h IZI,,-II\-\ 9

(AI6)

The factors CE and CL in equations (AI5) and (AI6) are

chosen to make the convergence tests consistent with the local

truncation error tests. In particular, CE = 0. 1 3E(q) and

CL = 3L(q)/2(q + 2), where 3E(q) and 5L(q) are the test

constants used, respectively, in EPISODE and LSODE for the

local error test (eqs. (A27) and (A28)) and where the variable

q indicates the method order.

If convergence is not achieved after the first iteration, the

codes anticipate the magnitude of 6,, one iteration in advance

by assuming that the estimates converge linearly. Thus, 6re+l,

which does not yet exist, is estimated by

=6 6,. =6,,,c,,,, _) t?l + [ ill-- --

6m- 1

where C,,, ( = 6,,,16,,,_ i) is the convergence rate. This assumption is used to modify the convergence tests (eqs. (AI5) and

(AI6)) as follows:

for EPISODE

,)

6/. _ CE EPS

where

6,i, = 6,,, min (1, C,I,)

C,_, = max (0.1 C,,,_ I, C,,,)

and for LSODE

where

a,;, = 6,, min (1, 1.5 C,',,)

C,I, = max (0.2C,,,_ I, C,,,)

Now, at least two iterations are required to compute C,,,. For

the first iteration, C,',, is set equal to 1 in EPISODE and equal

to the last value of C,, from the previous step in LSODE. For

the first iteration of the first step and after every update of

the Jacobian matrix, LSODE sets CA equal to 0.7.

If the corrector iteration fails to converge in three iterations,

h, is reduced by a factor of four if P is current and the step

is retried; otherwise, P is updated at y = El,,°1, and the step

is retried. The same corrective actions are taken by LSODE

if C,, > 2 after the second iteration. In the event of a singular

iteration matrix, both codes reduce h, by a factor of four and

attempt the solution with the new step length. The integration

is abandoned if either the step size is reduced below a minimum

value (both codes) or 10 convergence failures have occurred

(LSODE).

If the corrector converges after M (___3) iterations, an

estimate of the local truncation error is made, as described

below. For a BDF method of any order k, the local truncation

error, d,,(k) at t,, is given by

= b,k+ I

_d,,(k) ek+ t,, ,, .yl_+ J_(t,,) (A17)

where the variable k denotes the method order and the constant

ek+_ depends on the method formulation. For the variablestep method used in EPISODE (ref. 6),

k

j=|

ek+ I = (AI8)

(k+ 1)! l_.,,(k)

4o
where andforLSODE

_JJ-t"- t,, j j = 1..... k (AI9)

h,

and Ii.,(k) is the second element of the (k + l)th-dimensional

coefficient vector l,,(k) for the k'h-order method (see

eq. (AI 1)). For the fornmlation used in LSODE (ref. 28)

1

C_+ i - (A20)

k+l

The error d,,(q) in the qtCorder method (i.e., k = q) used

on the nth step (i.e., (t,,_ t,t,,)) is estimated as follows: As

discussed previously, the last column of z,,, z,,(q) contains the

q (q) vector h,,Y, /q. and that of z,[[)1, z,l,°l(q) contains the vector

h,lv_q) ,,,,,-=/q!. The difference of z,(q) and zl°l(q) gives

z,,(q) - zl01(q) =--h_[+l y,}q+t) + O(h_[+2) (A21)

q.t

by using the mean value theorem for derivatives. From

equation (A10) the above difference is seen to be equal to

lq.,,(q) e,,, which, upon substitution into equation (A21), gives

h;[+ly, lq+ t) = q! /q.n(q)e,, (A22)

if higher-order terms are neglected. This approximation is used

in LSODE. EPISODE, however, takes into account errors in

the past values and uses the following expression (ref. 6):

1

hq+ly(q+L) = _ e,, (A23)

l! = It

where

'l ,( ) Io,,(q)

d,,(q) - q: q'' q e,, - ' e,, (A26)

q+l q+l

because q!lq.,,(q)= lo.,,(q) for the formulation used m

LSODE (ref. 28).

The local error tests used in the two codes are as follows:

For EPISODE

,)

< EPS

which, upon using equation (A25) can be written as

Oq

\ "__ I/'_

I el,i I - -

3L(q) EPS

,)

_< 1 (A27)

where the test constant 3_:(q) is given by

3E(q) = Ii.,,(q)

And for LSODE

j= 2 \t_t -- t,, i/ j

(,, ....

i=l \EWTJ/ _< 1

tt --

q

I-[

j=l

(q+l)! I + H 4,,-, - t,,_j/j

j=2

(A24)

and where _.i is given by equation (AI9).

By substituting the above expressions for h,,Y,} q+ll and the

appropriate equation (A 18) or (A20) into equation (A 17) (with

k = q) and simplifying the resulting expressions, we obtain

the following estimates for d,,(q):

For EPISODE

By using equation (A26), the above inequality can be expressed as

l / \2"_ I/2 \EWT,/ /

Dq - . < I (A28)

/

3L(q)

where the test constant 3L(q) is given by

1 1 + H -t, -_t,,_) e,, (A25)

d,,(q) - l=.,,(q) j=2 "t"-I - t,,_)/j

q+l

3t_(q ) -

In.,,(q)

41
If the error test tails, the following corrective actions are

taken. In EPISODE h,, is reduced so that equation (A27) is

satisfied (see eqs. (A29) and (A30)), and the step is retried.

If the results with the new step length do not pass the error

test, h,, is reduced by a factor of five. After three and more

error test failures EPISODE reduces the step length by a factor

of 10 and reduces the method order by one if it is greater than

one. If an error test failure occurs with q = !, the Nordsieck

history matrix z,,_ i is reconstructed from Y,,_ i andf(Y,,_ t).

After the first error test failure, LSODE reduces h,, and/or

q by one and then retries the step. If the error test is again

not satisfied, h,, is reduced by a factor of five. After three and

more such failures, the method order is reduced to one if it

is greater than one, the step size is reduced by a factor of 10,

and the step is retried with a new Nordsieck history matrix

z,,-i, which is constructed from Y,,-t and f(Y,,_0. Both

codes abandon the integration if h n is reduced below a

minimum value. The maximum number of error test failures

allowed is seven in EPISODE and 10 in LSODE, after which

an error exit is taken.

If the error test is passed, the step is accepted as successful,

and the entire Nordsieck history array, z,,, at t,, is updated by

using equation (A10).

Periodically, both codes attempt to change the step length

and/or the method order to minimize computational work while

maintaining prescribed accuracy. After every step on which

no convergence test or local error test failure occurs,

EPISODE attempts to use a larger step length at the same

method order The new step size h'(q), where the variable

q denotes the order to be used on the next step, is chosen such

that it exactly satisfies the local error bound (eq. (A27)) by

assuming that the highest derivative remains constant. Then,

because d,, varies as hq+l (eq. (A17)),

I

h'(q) (_q) q+l r_,me = - (A29)

h,,

where r is the ratio of the step length to be attempted on the

next step to its current value and the subscript "same"

indicates that the current order (q) is to be used on the next

step. To allow for inaccuracies in the error estimate, certain

safety factors are built into the calculation procedure for h '(q)

to produce a smaller value than that given by equation (A29).

The formula used in EPISODE for rsame is

(A30)

rsame = 1

1

(5Dq) q+ + 10 -6

To increase the efficiency, both codes consider changing

the method order to q - 1or q + 1 at periodic intervals. After

an unsuccessful step or when the current order equals the

maximum order, q,,_x, the choice q + 1 is not considered.

Also, if q = 1, the choice q - 1 is rejected. For each method

order q' the step size h' (q') that will exactly satisfy the local

error bound is obtained by using the procedure outlined above

for q' = q (eq. (A29)).

For the case q' = q - 1, d,,(q- 1) varies as h_() '(q) (t,,) (eq.

(AI7)), which is equal to q!z,,(q). The local error test for

q' = q - 1 is as lbllows:

For EPISODE

zi ,,(q) "

i=1 9

Oq -1 =- <_ I

3E(q--1) EPS

where Zi,,,(q) is the ith element of z,,(q) and (ref. 6)

And for LSODE

I1(q -- 1)

3_:(q- 1) -

q-I

I-I _J

j=l

Oq-I _

(1 _ (zi.n(q)_2_ '/2

<_1

3L(q--1)

where (ref. 28)

1

3L(q- 1) - --

(q-1)

The step length ratio, if the order is to be reduced to q - 1,

is then given by

1

h'(q-1)- (Z) q (A31) rd°wn = h,,

where the subscript "down" indicates that the order is to be

decreased. If q = 1. rdo,,. is set equal to zero because q

cannot be decreased.

For the case q' = q + 1, d,(q+ 1) varies as hq+2_,(q+2J(t .., _ ,.,,,,

which is estimated by differencing the quantity h q+ I_y(q+ 1)

over the last two steps and then using the mean value theorem

for derivatives. For EPISODE equation (A23) gives

42
/ _ \q+l

Iy(q+l) _ hq+ly(q+l)

%, \h,._ 1/ _,,-i

=-.,hq+ 2y(q+ .,_ 21 + O(hu q+3)

where %, is given by equation (A24). For LSODE, equations

(A22) and (A26) show that the approximation for hu+lyIq+ 11

is given by lo.,,(q)e,,. Because the methods used in this code

are based on a constant step size. the quantity lo.,,(q)[e,, - e,,_ i]

gives

Io.,,(q)Ie,, - e,, i] -= ]/q+2y_lq+2) + O(hff +3)

The local error test for q' = q + 1 is given by

For EPISODE

V)'.- f/',.h._ _/

Ynlax i 9

<_1

3E(q+ 1) EPS

where (ref. 6)

3E(q+ 1) =

(q + 2)ll,,(q+ 1)

_q+ I

And for LSODE

Dq+ 1

] t "7\1/2

_<1

3t(q+ I )

where (ref. 28)

q+2

3c(q+ I ) -

to.,,(q)

The step length ratio, if the order is to be increased to q + I,

is then given by

= h'(q+l) - ( -- 1 _ q+2 (A32)

rup h,, \Oq + 1/

where the subscript "up" indicates that the order is to be

increased. If q = q,,,_ and in LSODE after a failed step, rup is

set equal to zero to prevent an order increase.

For reasons given previously certain safety factors are built

into the step length ratios (eqs. (A29), (A31), and (A32)). The

formula used in EPISODE for r_,,, c is given by equation (A30);

the other two ratios are computed as follows:

rdo_n I

(59<t_1) q+ 10-6

rup = I

)q+2 (10Dq+ i + 10 -6

The formulas used in LSODE to calculate the step length ratios

are

rdo_n -----

1

1

1.2 l(Dq)q+l + 10-61

#'up =

I

1.2 I(Oq+l)q+2 + 10-61

The order corresponding to the maximum step length ratio

r = max(r0o_ ., r_u_,rop) and the step length ratio r are selected

to be attempted on the next step if, after a successful step, r _> 1.3

(EPISODE) or i.I(LSODE); otherwise, both changes are

rejected. After a failed step, the order is decreased in LSODE

if rd,,,_.n> r_me; however, r is set equal to one if it is greater

than one. Several additional tests are performed on r before the

step length to be attempted next is selected. These tests may be

summarized as follows:

For EPISODE

r_min Ihm_ ( hu,m )1 t.h. , r,,_. max r, h-_u' rmin

43
where the arrow denotes the replacement operator.

And for LSODE

r -- min /I,.

h,nax

h,, ' rmax,

max (r, ++)1

In EPISODE, hmi n and hnm X are set equal to, respectively, ho,

the user-supplied value for the step length to be attempted on

the first step, and 10(tom - tout.ore) where tom is the current

time at which the solution is required and toot,,,m is the

previous value of tout. On the first call to EPISODE tom is set

equal to to, the initial value of t. In LSODE, however, hmi,

(default value = 0) and h,,,,x (default value = oo) are usersupplied optional input parameters. The quantity rmin (used

only in EPISODE) is set equal to 0.1, and r,,_x depends on

the code. In EPISODE, r,,,,x is set equal to 10 for the first

10 integration steps; thereafter, it is set equal to 1.5. In

LSODE, rmax is normally set equal to 10; for the first step

length increase following either a convergence or local error

test failure, it is set equal to two. In both codes for the first

step length increase for the problem rm__ is set equal to 104

if no convergence or error test failure has occurred.

After the step length ratio r has been computed, the step

length h' to be attempted on the next step is given by

adds a column of zeros, representing =,(q + 1), to z,,.

LSODE sets Z, (q + 1) equal to h_]+ i y,lq+ I I/(q + 1)?, which

is equal to lq.,,(q)e,,/(q + I) (see eq. (A22)). Both codes then

rescale all q + 2 columns of z,, by powers of h '/h,, to account

for any change in the step size.

The solution values at prescribed output times t,,m.i, t,,m._....

are obtained quite easily from the history array. For each

output station t.m, the codes continue the integration until the

first mesh point n for which t,, -> tom and then compute the

solution at t,,ot by a q,;+lm-order Taylor series expansion

about t,,:

Y(tout) =

q,;+n q'_+ (t,,_, - t,,)_ /t,_ - t,,\ j .

j=0 j=0

(A33)

where q,i+n and h,i+l are. respectively, the method order and

step size to be attempted on the next step.

Both codes start the integration with a single-step, first-order

method because information is available at only the initial

point, to. The Nordsieck history matrix z0 at to is constructed

from the initial conditions y(t0) and the ODE's as follows:

z0(0) ----Yo = y(t0): Zo(I) -= h0_'u = hal(Y0)

h t _ rh91"

Changes in method order (and step length in LSODE) are

attempted only after S successful steps with the same order (and

step length in LSODE), where S is normally set equal to q + 1.

However, if an unsuccessful step occurs, the step length and/or

order may be reduced. Following a failed error test or a failed

convergence test_ if P is current, EPISODE resets S equal to

2 if it is less than two, but LSODE resets S equal to q + 1

irrespective of its current value. If three or more error test failures

occur on any one step S is set equal to five in LSODE and either

q + 1 (ifq > 1) or 10 (ifq = I) in EPISODE. Following a step

for which the method order is not changed EPISODE sets S equal

to 2. If method order and step length changes are rejected because

r < 1.1, LSODE sets S equal to 3.

After every S - 1 successful steps, if q < q,,m_, EPISODE

saves e and _, and LSODE saves e, so that rup can be

computed. To minimize storage requirements, the vector e is

saved as the (qmax + 1) th column of z.

If the step size and/or method order is changed on the nth

step, z,, has to be modified. For the case q' = q, h' _ h,,, the

j;h column (j = 0,1 .... q) is multiplied by (h'/h,,) j. For the

case q' = q - 1, h' _ h,,, the last column of the old z,, is

ignored because it is not needed on subsequent steps. In

addition, EPISODE adjusts the first q columns to reflect the

reduced set of data represented by z,, (ref. 6). In both codes

the above scaling by powers of (h'/h,,) is performed on the

first q columns. For the case q' = q + 1, h' _ h,,, EPISODE

where ho, the step size to be attempted on the first step. has

to be supplied by the user to EPISODE. In LSODE. however.

h0 is an optional input variable and is computed by the code,

unless the user has specified a value for it.

GCKP84

GCKP84 is a general-purpose chemical kinetics code

designed to solve a wide variety of problems (ref. 15). It uses

the integration technique developed by Zeleznik and McBride

(ref. 18). As implemented in GCKP84 the integration

algorithm is an extensively modified version of the GEAR

package (ref. 19), which is similar to LSODE. In particular,

GEAR includes the two linear multistep methods discussed

previously. The methods are based on a constant step length,

and the method coefficients i/,,I (eq. (A11)) have the same

values as in LSODE. Hence, l_,(q) is a function of only the

current method order q, satisfies/0. , (q) =/3o.,, (see eq. (A7))

and ll.,,(q)= 1, and has to be recomputed only when the

method order is changed. GCKP84 uses the same two linear

multistep methods but the maximum method order is different:

11 for the implicit Adams method and 8 for the BDF method.

The methods are also implemented differently as discussed

below. For reasons given previously we restrict discussion to

the BDF method (eq. (A7)).

As in EPISODE and LSODE, GCKP84 maintains the

solution history in the form of the Nordsieck history array,

z. The array z,, at the current time t,, is obtained by using a

44
predictor-corrector process. The prediction step is performed

in two stages. First, an initial estimate for z),°l is computed

via equation (A9): that is, the result of the prediction step used

in LSODE and EPISODE (and GEAR) serves only as an initial

estinlate for z,l"1 in GCKP84. Second, the above result is

modified by means of an expression similar to equation (A 10),

as follows: The difference (=,,(1) -:,l,I)l(1)) (equal to e,, in

GEAR because li(q) = 1, eq. (AI0)), that is.

P is used fi_r a maximum number of 20 steps. At each iteration

the approximation h,,£1,'''1to h,,f, is computed by using equation

(AS). If any Y,!_I'l < 0, the iteration is abandoned. The step

length is then reduced by a factor of two, and the step is retried.

Convergence of the estimates is said to occur if any of the

following three tests, which are applied in the order they are

given here. is satisfied. The first test involves the magnitude

of the successive differences (}>,/,,,I_ y,/,,, II):

e,, = h,,};',, - h,,_',l,Iq (A34)

may be regarded as the error in the qth-order predictor relative

to the converged array z,,. Equation (AI0) gives the history

matrix z,, by adding the remainder term associated with using

the qm-order predictor. Of course, since e,, can be computed

only alter the converged solution is produced at 6,, the above

procedure cannot be used. However. since e,, i is available.

it can be used to improve the initial estimate given by equation

(A9). However, for additional accuracy improvement.

GCKP84 uses the quantity E obtained by accumulating the

errors {e,,l (see eq. (A35)). The quantity E,, may be regarded

as the estimated global error in h,,£1,,°1. Since E,, is not known

at the start of the step, GCKP84 uses E,,_ i to improve the

estimate given by equation (A9) as follows:

/_ \q+ I

z,l"l- z,l'u+ i'"'q &_J,,(q)

where hL. which is normally equal to h,, _. is the step size

that E,, i is based on and the term (h,,/h__) q+l accounts for

this fact: the exponent q + 1 arises because the current order

is q and the local error varies as tFj+l (see eq. (AI7)). On

the first step, E,, i (equal to E0) is set equal to zero because

_ (equal to f(y(t0))) is known exactly.

After the prediction process is performed, the code checks

the [ yl_,_l]for negative values. Because it is physically imtx_ssible

for species concentrations, temperature, density, or velocity

to be less than zero, the results of the predictor step are rejected

if any _!_]l < 0. Also, for each variable i lbr which the above

condition is obtained, &.,,_ i is reset to zero if it is less than

zero. The step length is then reduced by a factor of two and

a new z},"1 is generated. The above procedure is repeated until

either all predicted solution components are nonnegative or

the step length is reduced below a miniraum value, h,,,i,, in

which case an error exit is made•

To correct the initial estimate GCKP84 includes a variety

of iteration techniques. For reasons given previously the

discussion is restricted to the Newton-Raphson method. The

procedure used to generate the improved estimates y},,,,I

(m = 1,2 .... ) is exactly the same as that described for LSODE:

solve equation (AI2). The iteration matrix P (eq. (AI3)) is

only accurate enough to achieve convergence, but the same

(_ ttt -=X i=1

--]elm 11t2tl/2 ')

t,iJ I

where Y,,a,.i is given by the expression used in EPISODE for

semirelative error control (eq. (A 14). The second test is based

on the size of the current estinaate for e,, relative to the size

of the current estimate for E,, (see eqs. (A34) and (A35)):

h,,g!;;'J- h,,r,);,'J - ,

- <0.1

/_, \q+ I

- V"q

• V'LI

If for any i the denominator in the above summation is less

than 10 5o, it is set equal to 10 2. The third criterion is

based on how rapidly the iteration is improving the solution

and is given by

"_<10

I 6,. I

which can be applied only after two iterations. However. the

third test is applied only alter five iterations and that too only

if 6,,, _< 5.

If convergence is not achieved after lour iterations, the

iteration matrix P is updated at v = _oj and the correction

process is retried. This procedure is repeated four times, alter

which, if the estimates have not converged, the step length

is reduced by a factor of two and the step is retried. The same

corrective actions are taken if on the filth or subsequent

iteration 6,,, > &,,,_ 2. The above cycle of updating P every,

four iterations and then reducing h,, by a factor of two after

h)ur such updates is repeated until either convergence is

obtained or the step length is reduced below h,,,,, in which

case an error exit is taken.

After corrector convergence the local error test is applied.

This test is based on E,,, which is estimated by using

(_ \q + I

E,, = e,, + ,,,q E,, t

V',../

(A35)

45
and can be written as

A,_ i = I \ Y,n,_×,// 9"

--- < 5 (A36)

Dq - EPS 5_(q) EPS 5c(q)

/)q+l

_

E,.,, \hJ

i = I Ymaxi

EPS 5c(q+ 1)

where EPS is the user-supplied local error tolerance and

5_;(q) (equal to 2/lo.,, (q)) is the local error test coefficient for

order q.

If the error test fails, the error vector E,,_ 1 is updated by

using equation (A35), and hE is set equal to h,, because E,_ i

is now based on h,,. The code GCKP84 then reduces the step

size and/or the method order by one and retries the step. After

three and more error test failures, the method order is reduced

to one if it is greater than one, and the step length is set equal

to hmin. A new Nordsieck history matrix at t,,_ i is constructed

from Y,,_ I and f( Y,,_ I), E,,_ z is set equal to zero, and the

step is retried. After seven such failures or if h,, is reduced

below h,nm, the integration is abandoned and an error exit is

made.

If the error test passes, the step is accepted as successful,

the Nordsieck history array z,, is updated by using equation

(A10), E,, is computed by means of equation (A35), and hE

is set equal to h,,.

To increase the efficiency of the integration, the code

periodically considers changing the method order to q - 1 or

q + I. Of course, ifq = 1, the choice q - 1 is not considered.

After an unsuccessful step or if either q is equal to the

maximum method order, qmax, or Dq > 4 tG(q), the choice

q + 1 is rejected. For each method order q' the step size

h'(q') is computed from an estimate of the local error in a

manner similar to the procedures used in EPISODE and

LSODE (eqs. (A30) to (A32)). For each method order q'

GCKP84 computes the step length ratio r(q') as follows:

h'(q') 1

r(q') - ---- = (A37)

h,, 1

4 '+1 -6

-_(Dq,D:) q + lO

where

Oq_ I

YFIla_. t

EPS 5c(q-1)

and

The local error test coefficients 5_;(q-I) and 5d(q+ 1) for

orders q-1 and q+l, respectively, are given by

3c(q-1) = 2

and

5_(q+ 1) = 2(q + 2)/Io.,,(q)

The quantity D: in equation (A37) is set equal to 10, unless

,,x,, i >- 10 25/xfN, in which case it is set equal to 0. t. If

a,, is also greater than 10-25/x/N, D- is set equal to a,,/

a,_ i. Finally, if D_ is less than (0.25) q+l, it is set equal to

this quantity.

The order corresponding to the maximum step length ratio

r = max(r(q- 1), r(q), r(q+ 1)) and the step length ratio rare

selected to be attempted on the next step if r _> 1.1 after a

successful step; otherwise, both changes are rejected. After

a failed step, q is decreased if r(q- 1) > r(q); however, r is

set equal to 1 if it is greater than 1. The following additional

tests are pertbrmed on r before the step length h' (equal to

rh,,) to be attempted next is selected:

r--min I hm_ - ( hmi'_)l (. hn , /max, max r,

The minimum, hmin, and maximum, h,,,ax, step sizes are,

respectively, set equal to h0, the user-supplied value for the

step length to be attempted on the first step, and 10(to_,¢-

tout,old). On the first call t,,_t is set equal to to. The quantity

rr_ is set equal to 10. For the first step length increase

following either a failed convergence test or a failed error test,

it is set equal to two. However, after three or more error test

failures, it is set equal to min (104, h,,/hmin), thereby ensuring

that the new step length equals hmi n. For the first step length

increase for the problem, r,._,_ is set equal to 104 if no convergence or error test failure has occurred.

Changes in method order and step length are attempted only

after S successful steps with the same order and step length,

where S is normally set equal to q + 2. However, if an

unsuccessful step occurs or if D u > 4 5c (q), the step length

may be changed, and the method order may be reduced.

Following a failed convergence or local error test. S is set equal

to q + 2. After three and more error test failures, S is set equal

46
to three. If method order and step size changes are rejected

because r < 1. I, S is set equal to 10. Finally, the successful

step counter is increased by one only if convergence is obtained

in eight or fewer iterations.

If the step size and/or the method order are changed on the

;?h step, z,, has to be modified. For the cases q' = q - 1 and

q' = q, the modifications are made exactly as in LSODE

(described previously). For the case q' = q + 1, z,, is first

augmented by a column containing the vector /q.,,(q)

E,,l(q + 1), which is approximately equal to h_[+tY, lq+tl

/(q+ 1)!, and then all (q + 2) columns are scaled by powers

of h'/h,, to account for any change in the step size.

The solution values at the prescribed output times Gut.i,

Gt,.2 .... are obtained from the Nordsieck history array by

using the Taylor series expansion meth(xt (eq. (A33)) described

tbr EPISODE and LSODE. The same procedure used in these

two codes to start the integration is used in GCKP84: the step

size, ho, to be attempted on the first step must be supplied by

the user.

CHEMEQ

In this technique, developed by Young and Boris (ref. 10),

the species rate equation (A1) is expressed as a difference

between two positive-definite terms as lbllows:

(/3';

-- = f = (Pi - :D_ i = 1 ..... N s (A38)

dt

where, for species i, the production rate (Pi and the destruction

rate 33 i can be derived from equation (3):

(Pi = p

NR

' j:_lNk (uo'R-i + vi7 Ri) 1

' E (u,;R j+ vi;'R_j)]

j=l ../

i= 1 ..... N s

(A39)

When the temperature ODE (eq. (9)) is required (method B),

it can be cast in a similar form by combining equations (9)

and (A38)

N_. N s

- -

dy% + 1 dT k = 1 k = I

dt dt N s N,,;

k=l k=l

= 0P','s+,. I -- 3).vs+ I.

= (Pr -- :Dr,

where

N5

E33 a*I_

(p/, _ k= I

N_

E Y#?I,,I,

k=l

and

, ),

X)T_ k=l

E ),k(,l,,_

k=l

The objective of this decomposition is to enable factorization

of Yi from X),

_c), = "_iYi = y,/ri

where £;, the loss coefficient for species i. is obtained simply

by dividing 33, by yi (i.e., £i = X)i/Y;). With this new

notation, equation (A38) can be written as

(t3' '

-- (P, - 32iy i = _; - yi/ri (A40)

dt

which, fi, r constant 07 i and £i, can be solved to give

Vi(t,,) Vi(t,,-i +h,,) (Pi (P; = = -- "+- Yi(tn - I) -- exp(-£,h,,)

Expressed in this way, it can be seen that 1/£ i ( = r_) describes

how quickly the variable y, reaches its equilibrium value.

In advancing the solution from time t,,_ t to time t,,, all of

the equations are separated into two classes, stiff and nonstiff,

according to the criterion

h,, ( >_ I stiff

< 1 nonstiff

Ti,n- I

where ri.,, _ denotes the value of ri at time t,, i. The two

types of equations are integrated by separate predictorcorrector schemes. For equations classified as nonstiff, the

improved Euler method (with the Euler method as predictor

and the modified Euler method (or trapezoidal rule) as

corrector) is used. For equations classified as stiff, a simple.

stable, asymptotic formula is used.

47
Nonstiff predictor:

ylOl I, II _ Y,. ,, -- l + hJi.,, _ I (A41 a)

Stiff predictor:

y lOl _ Yi.... i (27-i.,,_1 - h,,) + 2h,, ri.,,_l(Pi.n_ I (A41b)

,tl

Nonstiff corrector:

2"ri.n- I + h,,

After each step n the step size h,i+ Lto be attempted on the

next step is computed from the converged integration cycle

as follows:

h,',+l = h. [' ]

+ 0.005

(o/EPS) b'2

The step size, ho, to be attempted on the first step is

determined such that none of the variables will change by more

than a prescribed amount. The fornmla used for ho is

h,, [ + fl,,,I] yl,,,+tli.,,= Yi.,-t + _- f.,,- t a,.,, t

(A42a)

Stiff corrector:

,.,, = ri.,, + ri.,,-i [ . + (Pi.,,-i

7"I''1 -- h. + r/,._ 1 + hn + Y/,n- I i. + "ri,,- 1 .I_ _.n

(A42b)

In equations (A41) and (A42), m + I is the current iteration

number. The zeroth estimate is the result of the predictor step.

Also, fJ_,,l =f(l_',_ll). Convergence is ascertained by comparing .Y!''+,,,, ]1 with y),,,I,.,, for all N components using the relative

error criterion

I y l,,+ I] ylml I 9

max ,, - i.,

o = < EPS (A43)

i min 1y,l_,l . ' yl,,,+ i., tI [

To avoid numerical difficulties with the use of equation (A43),

each estimate is constrained by a minimum value. In the

present work, a variable that is less than 10 -2o is set equal

to 10 -2°. Thus, for a species with decaying concentration,

convergence is obtained trivially once Y/ < 10 -2°, and its

equation is decoupled from the equation set.

If convergence is not achieved after ITMAX iterations, the

step length is halved and the step repeated. In this study, a

value of ITMAX = 5 was used because it minimized the

execution time for both test problems (ref. 2). If the corrector

converges after M iterations (M < ITMAX), the step is

accepted as successful, and the solution is updated

Yi ,, = yIMI i = 1 ,N.

, I,tl ' • • "

No attempt is made either to estimate or control the local

truncation error.

h0 EPS min [yi(t0) 1 ]

= --, or- if yi(t0) --< l0 -2f_

' Lf(t0) £i.o

The solution at each output station t,,_t was computed by

linear interpolation between the computed approximations at

t,,__ and t,,, where t,, is the first mesh point that is >tout:

Y(t,,u0 = Y,,-t +

tout -- In- 1

htt

(r',,- E,-_)

CREK1D

In CREK1D, attention is paid to the distinguishing physical

and computational characteristics of the induction, heat release,

and equilibration regimes (refs. 11 to 14). This code consists

of two algorithms developed for the two distinctly different

regimes: (a) induction and early heat release, when the ODE's

are dominated by positive time constants and (b) late heat

release and equilibration, when the ODE's are more stable

(ref. 2). Both algorithms are based on an exponentially fitted

trapezoidal rule, but they use different iterative methods for

convergence.

The code CREKID solves a mixed differential-algebraic

system of equations: ODE's for the species mole numbers

and the algebraic enthalpy conservation equation (8) for the

temperature. The ODE's and algebraic equation are solved

simultaneously; however, in the following discussion the

variables y and Y refer only to the species mole numbers.

The solution method used for the species ODE's is a generalized, tunable, single-step, implicit procedure:

i= 1..... N s

(A44)

where Ui.,, is a degree-of-implicitness factor. This parameter

is obtained by "exponentially fitting" it to a locally exact

solution of equation (A l) as follows: The species rate expression

f is expressed in a locally linearized form such that

f =f.,,-I + Oi.,,(Y, - Yi.,,-]) i = 1..... N s (A45)

48
wherethechoiceof 0j.,,, a suitable linearization constant, is

discussed shortly. Equation (A45) assumes that in the interval

[t,, ). t,,] (i.e., locally) each species mole number varies exponentially. Integration of this equation gives

Y,.,, = Y,.,,-i + h,_.,,-t

exp(h,,0i.,,)- 1

hnOi,n

(A46)

To exponentially fit U,.,, we first replacer.,, in equation (A44)

with the expression obtained from equation (145):

f.,, = fi.,, i + Oi.,,( Yi.,, - Yi.,- t) (A47)

and then equate the resulting expression for Y,., with equation

(A46). These operations give

I 1

U_.,,- + i= 1..... N s (A48)

h,,Oi.,, 1 - exp (h,fli.,,)

In order to maintain absolute A-stability of equation (A48)

(i.e., to keep errors introduced into the numerical solution at

any one step bounded as h,, is increased indefinitely), Ui.,,

must be restricted to the interval (0.5,1.0). For values of

0,.,, > 0, equation (148) gives Ui.,, < 0.5. CREKID resolves

this problem by setting Oi.,,= 0 whenever it is greater than zero.

This value of Oi.,, gives Ui.,, = 0.5, so that equation (A44)

defaults to the second-order-accurate trapezoidal rule. However,

for Oi.,,< O, equations (A44) and (A48) together are equivalent

to the locally exact or exponential solution, which has an

equivalent polynomial accuracy of order six to eight (ref. 11).

Thus, equations (A44) and (A48), with the constraint

(0.5 _< Ui.,, < 1), constitute an exponentially fitted trapezoidal

rule, a method which is A-stable and has a polynomial-order

accuracy of at least two and as great as six to eight.

The linearization constants lOi.,,t are obtained in one of two

ways. In the first, called functional linearization (see refs. 11

to 14), equation (A47) is solved explicitly for 0,.,, to give

fi,n -- f,tt - I

0_.,, - = Z,.,, (149)

Y,.,,-Y,.,,-,

In the second approach, called formal linearization (refs. 11

to 14), the net formation rate of each species is expressed as

a difference between two positive-definite terms, as described

in the previous section (see eqs. (A38) and ('A40)). Comparing

equations (A47) and (140) gives

Oi.,, = - _i .... 1 (A50)

for this procedure.

At each integration step, equation (A44) must be solved for

Yi.,,. The solution is accomplished by Newton-Raphson (NR)

iteration in regime b and Jacobi-Newton (JN) iteration (ref. 29)

in regime a.

A Newton-Raphson functional F,.I','1 (i = 1..... Ns) for each

species mole number is defined from equation (A44) by

( )

Flml _ ,.. - Yi,,-, i 1 -- Ui,'! " £1ml

.... h,,Ui,,, U,_ f., , --Ji.,,

(A51)

F [ml for i = 1.... N s. For temperature the functional r.,, is

defined from the enthalpy conservation equation (8) as

F//,,I vt,,,I

.... = E "_." I'_(TI,,'''1) - Ho(To) (152)

k=l

where m is the iteration number. _,,,I is the re'h-approximation

to the exact value T(t,,). J{k(ZI,,''1) is the molar-specific enthalpy

of species k at temperaturc T,I'''1, and H0(T0) is the initial

mixture mass-specific enthalpy at the initial temperature To.

Newton-Raphson corrector equations with log variable corrections (for self-scaling of the widely varying mole numbers)

are given by

N>; OFtl;;,i off _l;n]

k=l am IL, , OIn _t,/,,I

for i= 1..... N s.

Aln _})m+ I I = _FI.d I,H

(A53)

Ns _[m] .:.q_ [m]

E .,7- _'_l',nqT;,I Aln i/tin , _.,,+ I1 + --u'_r,,i

k= I oin xi,,,, Oln TI'''1.,;

A]n 7_nm+ I1 _,1.,1

(A54)

where

AIn Yi!;;'+'1 = In Y,.[;;'+'}- In _i;; I i= 1..... N s

and

Aln TI,/'+ll = In _/,,+,1_ In TI,,''1

The partial derivatives in equation (A53) are obtained from

equation (A51) and are as follows (with the step and iteration

numbers suppressed for clarity):

0In Yk -- Y_ \hU, cOYaJ

(A55)

4O
where 6it, the Kronecker symbol, is

(_ik =0 i¢k =1 i=k

and OftOYk and OflOT can be derived from equations (3) to

(8), and (10). In evaluating Of/OY_, the partial derivatives

with respect to o,,, are assumed to be negligible compared

with the other terms. The required partial derivatives are then

given by

NR

- (plr_) -I E (v_ - v_) (%R/ t - %tt R_j) (A56)

OY_ j=l

_ _f,

OT T p-Tj:_ _'/j - "_/) " + 3T_ .j

-R _ N_j+T-n _

where

Ns Ns

n/= E v,_ n,"= E vo"

i:l i=1

The partial derivatives of Fr are obtained by differentiating

equation (A52) and are as follows:

OFT

- Yk_k,

01n Yk

NS

3Fr

- E Y_)"kT' 0In T

i=[

(A57)

where, again, the step and iteration numbers have been

suppressed. The Ns + ! equations (A53) and (A54) are

solved simultaneously by LU decomposition and backsubstitution (ref. 27). The resulting log variabJe corrections

are used to update the current estimates I_,'_ + Ill and _,'"+ _1

by the approximate equations

yil,,,+11 =,i,, 1 +Aln i= 1..... Ns

(A58)

The solution procedure does not use a predictor; instead, the

converged results In.,,-tt and T,, i from the previous step are

used to start the iteration.

The JN iteration technique can be derived from the NR

iteration procedure by neglecting the off-diagonal elemcnts of

the Jacobian matrix for the mixed differential-algebraic system

of equations. With this simplification, equations (A53) to (A55)

reduce to

Oln y.l,,,I

I J1

-- A[n El,,,+ i.n II = -F .,,,,,I i = 1..... Ns (A59)

o_.l,,,I _/,, + tl awIml )

--I'.,_ Aln

0In _,,d = --'T.n

aln

aiO'"l

_"

Y,!l;

i,n

'l

_

h,,U,.,,

yl,,,i

i.n yH

*'" aY_l;;

0(..i,,,i

a_.n

'l. 3

(A60)

The iteration procedure is further simplified by approximating

Of/OY,., (eq. (A56)) as follows:

af,

an

NR

--- -(PYi) ' E (vi;Ri + vi'/R .i)

)=1

which, when combined with equation (A38), gives

Of 33i

or',. r,

With these simplifications equation (A59) can be solved

explicitly for the iterative corrections

F.Iml

Aln Yt''+ II = ',"

- "" ,,d ff)l,.I _.,, //,,,u,.,, + _,.,,

i = 1..... Ns (A61)

The temperature correction is obtained by substituting equation

(A57) into equation (A60):

FI,,I

 ̧]',,1 (A62) Aln _1,,,,+Jl =

Ns

TI[.I E v[.,I

k=l

where cp,_(El,,'''1) is the constant-pressure molar specific heat

of species k at temperature _,,,I. The current estimates are

updated by using equation (A58).

To start this iteration process, the predicted values for the

species mole numbers are given by equation (A46):

_o,)= n.,,-, + h,Z,,,_,exp(h,,0,.,,) -- 1

httOl,n

i-- 1..... N s

(A46a)

5o
The predicted temperature is obtained by a single NR iteration

of the enthalpy conservation equation (8)

El,,°l = T,, i +

N_

H,,(r,i - Y'J

k=l

N3

YL' (r,, ,)

k=l

(A63)

For both NR and JN iteration schemes the test fi)r

convergence of the estimates {vl,,,l_ ,,.,, is based on the magnitudes

of the log variable corrections, and is given by

Ns _'2

(AIn y,,,,I,,,I)2

6,. = i= I '_ _< EPS (A64)

Ns /

This test is used only with variables whose magnitudes are

greater than 10 2(_; that is, the summation does not include

species with mole numbers _< 10-2°. At each iteration the

estimated convergence rate, C,,, defined as

is also computed. If convergence is not obtained after ITMAX

iterations, where ITMAX is the user-supplied maximum number

of corrector iterations to be attempted, or if C,,, > 0.8, the

step length is decreased. The new step length is calculated as

follows:

h,,- h,, rain 10.5, max(0.1, 0.5/C,,,)I

and the step is retried with the new step size. At least two

iterations are required to define C,,,; on the first iteration 6,,,

is set equal to 10. A value of ITMAX = 10 was used for both

problems examined in this study.

If convergence is achieved in M iterations (M _< ITMAX),

the step is accepted as success[ul, and the solution is updated:

Y,,, = v>II i= 1 ,Ns

, - I,It ' • • "

T,, = T,, 1''41

After corrector convergence, the step length, hirer, that

would produce a convergence rate in the range (0.4,0.5) is

estimated as follows:

biter = h,,(0.4/C,,,)i/2 C,,, < 0.4

=h,, 0.4 _< C,,, _< 0.5

=h,,(0.5/C,,) _/-_ C,,, > 0.5

If convergence occurs on the first iteration. C,,, is set equal

to 0.1.

At each step an average weighted local truncation error

estimate, d,,, is computed by using the approxinmtions

d,, = 6

i=

for the JN iteration, and

d, = ._

i= t max (Y,,,_ I, g,,,))

for the NR iteration. The above summations include only

species whose mole numbers are greater than 10 2o. For both

iteration techniques the step length, h_ccv, that would exactly

satisfy the user-specified local relative error tolerance, EPS,

is calculated from

h,,.._.? = h,, ( EPS/d,, ) I/3

The step length h,i+l to be attempted on the next step is

taken to be

h,i+ I = min(hit_., h_c_-y. 10h.) (A65)

However, if convergence difficulties forced a reduction in the

step length on the current step. h,_+l is restricted to

h,',+i - rain(h,,, h,i+ I) (A66)

to prevent a recurrence of the problem.

CREK 1D automatically selects the linearization method and

the iteration scheme to be used for solving equation (A44).

During induction and heat release, when small step lengths

are required for solution accuracy, the JN iteration is used

to minimize computational work. During late heat release and

equilibration, when the ODE's are more stable and larger step

lengths can be used, the NR iteration is preferred since it has

better convergence properties than the JN iteration. The regime

identification test exploits the tact that during equilibration

many reactions achieve a condition in which the forward and

reverse rates are large but with vanishingly small differences

(refs. 13 and 30). The actual test used at the beginning of each

time step is

,)

Ill _ 10 3 ((P, + 3)_) (A67)

where (Pi and X), are given by equation (A39). If any two

species satisfy equation (A67), regime b is obtained, and the

NR iteration is used for the step. If fewer than two species

51
satisfy equation (A67), regime a is obtained, and the JN

iteration is used for the step. Once the NR iteration is selected

for any one step, the above test is no longer applied, and the

NR iteration is used for the rest of the problem.

Whenever the reaction rate for any species satisfies equation

(A67), that species is considered to be in "quasi-steady state"

and the "L-formulated" equation (A50) is used. For all other

species the "Z-formulated" equation (A49) is used. To minimize

computational work, the [Z,•,,I are evaluated only once per

step: at the beginning of the time step, using equation (A49).

However, since Y,.,, and f.,, are not known at the start of the

step, the [Z,.,,I are approximated using values from the

previous step:

Zi,tt _ Zi.n--[ _-

f,n-I -- f,n-2

Y_.n- I - Yi,n- 2

CREK1D also includes an algorithm for filtering the initial

conditions that may be ill posed. These ill-posed conditions

may arise, for example, in multidimensional modeling because

of the averaging of mole numbers over adjacent grid nodes.

CREK 1D therefore "filters" the initial conditions to provide

physically meaningful initial mole numbers and net species

production rates. For purposes of this filtering CREK 1D uses

the decomposition performed in CHEMEQ (eqs. (A38) and

(A40)). On the first call to CREK1D it uses this formulation

over one time step of length hi, given by

1

h_ - (A68)

max c_i,0

i

The predictor-corrector algorithm uses equation (A46) (with

0i.i = -£i.0) as the predictor

Yi.Ii°l = Yio + hlfo (]1 exp (-hl£i.o)_ i= 1..... N s

An implicit Euler corrector is then iterated to convergence

yl,,,+ II b, ¢l m+ II i,l : Yi,o 2¢- "Idi, l

In the above two equations, the subscript 1 indicates that this

is the first step. Using equations (A38), (A40), and (A58),

col,,,+It =- (pl,,,I and o_tm+II = together with the approximations ,.,,. i ,._ ,.I

£I_'_I, the preceding corrector equation can be rewritten to

provide the following expression for the log variable

corrections lain Yil,'_'+ Ill:

yIm+ II = Yl_o -- Yi!'_ I -t- hlftl'[ ']

Aln ,. I yl,,,I

,., + hL:_151

Equation (A69) is iterated until converged; that is, the criterion

given by equation (A64) is satisfied• If convergence is not

obtained after 10 iterations or C,,, > 0.8. the step length is

halved, and the step retried. If convergence is obtained after

M iterations (M < 10), the step is accepted as successful, the

solution for the mole numbers is updated

1I,.i = y_!]ll i= 1..... N s

and the temperature T_ is obtained by a single NewtonRaphson iteration

TI = T0 +

Ho(To) - ]_ Yk,l_'k(T0)

k=l

]_ Yk._q,,k(To)

k=l

The step size, h_, to be attempted on the next step is

determined from the maximum loss coefficient at t_ by using

an expression similar to equation (A68). For this step, the JN

iteration (eqs. (A46a) and (A61) to (A63)) is used, with all

0i.,, set equal to zero, so that all U,.,, = 0.5 (see eq. (A48)).

The predictor step (eq. (A46a)), therefore, reduces to the

explicit Euler method, and the corrector (eq. (A44)), to the

trapezoidal rule. For the next and subsequent steps the step

size is adjusted according to equation (A65) or (A66), and the

iteration procedure and linearization constants are selected as

described previously. If NR iteration is used, the Jacobian

matrix for the mixed differential-algebraic system of equations

is updated at y = Y,,_ 1, T = T,,_ lThe solution values at the prescribed output times t,,_t._,

to_t.2.... are obtained by adjusting the step length so that the

internal mesh points coincide with these times. Thus, the step

size h,',+t is given by

h,;+l -- min (h,;+l, t,,_, - t,,),

where tou t is the current value of the output time, and the

results at /out are generated by solving the governing

equations. To continue the integration past each output time,

the procedures described above for the second and subsequent

steps are used.

To reduce the computational cost, the use of exponential

functions is minimized by replacing them with rational function

approximations. For example, the term (e' - 1)/x in equation

(A46a) is evaluated by means of a (2,2) diagonal Pad6

approximation, e_2,2), for exp x:

X X 2

1+-+--

2 12

e-_2,2 _ -- x < 0

X X 2

i= 1..... N s (A69) l --+--

2 12

52
which gives

e'--I 1

x<0

x 1 -x(l/2-x/12)

Similarly, the tuning factor U,.,, (eq. (A48)) is evaluated by

using the approximation

Ui.,_ _ [ - _

I ex p (.t_ /.rl_

2 \ 12/

This equation requires six operations to evaluate and does not

exhibit the singularity at x = 0 of the exact expression (eq.

(A48)).

Although log-variable corrections are used in the code.

evaluation of logarithms of the variables is not required. Also.

the use of the approximations given by equations (A58) awfids

the cost of computing the exponentials of the log-variable

corrections m obtain the new estimates.

Another technique used in the code to reduce computational

work is to locally linearize the expressions for the thermodynamic properties of the species and the rate coefficients.

In particular, during the course of iterative convergence of

the equations, the thermodynamic properties and rate

coefficients are not reevaluated while the current temperature

is within a local window (T, T + AT), where ATis specified

by the user. Use of this strategy has been shown to reduce

the computational work (refs. 2, 3, and 5).

53
Appendix B

Description of Temperature Calculation

In this method, the temperature T,, at each discrete time 6, is given by

is computed from the solution for the species mole numbers

returned by the integrator by using the algebraic enthalpy

conservation equation

N S

i=1

Oi.n_i( E,) = Ho (8)

Equation (8) is solved for the temperature by using the

Newton-Raphson (NR) iteration technique (e.g., ref. 23). This

equation is rewritten as

NS

F(T,,) = y]_ o,,,,f.i(T.) - Ho (BI)

i=1

so that solving equation (8) is equivalent to finding the zero

of F. The quantity F( 7J,,''1) is the amount by which the mixture

mass-specific enthalpy at the m _happroximation for T,,, T,I'"I

(m = 1,2.... ), fails to satisfy equation (8). A new approximation,

T,I,,,+ II for the temperature is obtained from equation (BI)

by locally linearizing F at T,["I:

F(_, '')) F(_, "I)

(OF/O T) T= l']ff d N,,;

i=1

The test for convergence of the estimates is based on the magnitude of the corrections 87J,,'''+ II (equal to T,I'''+ll - T,I'"I) and

Method A

fiTl,.+_,,i II [ _,) ERMAX

where the vertical bars denote absolute value and ERMAX

is the local relative error tolerance. If convergence is not

obtained after MAXITS iterations, where MAXITS is the usersupplied maximum number of corrector iterations to be

attempted, an error exit is taken. If convergence is achieved

in M iterations (M N MAXITS), the solution TI,,MI is accepted

as the temperature at t,,:

T. : TI,)Mt

The NR iteration will converge if the initial guess (i.e.,

_01), is sufficiently accurate (ref. 23). The present work did

not utilize a predictor to generate ZI,°I; instead, the most recently

computed temperature was used to start the iteration. Now,

the temperature was evaluated at the end of each integration

step and whenever the species derivatives and Jacobian matrix

were computed. Hence, the converged value obtained either

at the end of the previous step or from the previous estimates

for the mole numbers was used as the initial guess for the current

temperature. For the very first temperature computation for

the problem, the initial temperature, T0` served as the predicted

value. The above procedure was found to be satisfactory in

that the iteration converged for all integration methods and

EPS values used in this study. In addition, the converged temperature was not significantly different from that obtained by

integrating the temperature differential equation (ref. 2).

54
References

I. Gear, C.W.: Nunmrical Initial Value Problems in Ordinary Dift_:rential

Equations. Prentice-Hall. Englewood Cliffs. NJ. 1971.

2. Radhakrishnan. K. : Comparison of Numerical Techniques lor Integratkm

of Stiff Ordinary Dillerenfial Equatkms Arising in Combustion Chemisl_.

NASA TP 2372. 1984.

3. Radhakrishnan. K.: Nev. Integration Techniques for Chemical Kinetic

Rate Equations. 1. Efficiency Comparison, Combusl. Sci. Technol.. vol.

46. no. I 2, 1986. pp. 59-gl.

4. Radhakrishnan. K: A Comparison of the ETficiency olNumerical Meth_ds

lilr Integrating Chemical Kinetic Rate Equations. Computational Methods.

K.L. Strange, ed.. Chemical Propulsion Information Agency', CPIA

PUBL 401, 1984. pp. 69-g2 (Also, NASA TM 83590).

5. Radhakrishnan, K.: Fast Algorithms Ibr Combustion Kinetics Calculations:

A Comparison. Combustion Fandamentals Research, NASA CP-2309.

1984, pp. 257-267.

6. Byrne. G.I)." and Hindmarsh, A.C: A Polyalgorithm Ior the Numerical

Solution of Ordinary Differential Equations. ACM Trains. Math.

Software, vol. I. no. I, 1975. pp. 71 96.

7. Hindmarsh, A.C.; and Byrne, GD.: EPISODE: An Effective Package

for the Integration of Systems of Ordinary Differential Equations.

UCID 30112, Re','. I, l,awrence Livermore Laboratory. 1977.

8. Hindmarsh, A.C.: LSODE and LSODI, Two New Initial Value ()rtlinary

Difli:rential Equation Solvcrs. ACM SIGNUM Newsletter, ,,,_fl. 15. no.

4, 1980, pp. I0 II.

9. Hindmarsh. A.C.: ODEPACK: A Syyslematized Collection of ODE

Solvers. Scientific Computing. R.S. Stepleman, ed., North Holland

Publishing Co.. Amsterdam, 1983, pp. 55 e,4.

10. Young, T. R. ;and Boris, J. P. : A Numeri,,:al Technique for Solving Stiff

Ordinary Dift;,:rential Equations Associated with the Chemical Kinetics

of Reacti,,e Flow Problems. J. Phys. Chem., vol. gl, no. 25, 1977,

pp. 2424-2427.

II. Pratt. D.T.: CREKID: A Computer Code lk*r Transient, Gas-Phase

Combustion Kinetics. Paper WSCI g3-21, Combustion Institute, Apr.

11 12, 1983.

12. Pratt, D.T.: and Radhakrishnan. K.: CREKID: A Computer Code for

Transient, Gas-Phase Combustion Kinetics. NASA TM 83806, 1984.

13. Pratt. D.T.: Fast Algorithms for Cornbuslhm Kinetics Calculations. Still"

Computation. R.C. Aiken. ed., Oxford University Press, New York,

1985, pp. 217 229.

14. Pratt. D.T.: and Radhakrishnan, K.: A Fast Algorithm for Calculation

of Combustion Chemical Kinetics, Paper 8 A 1. Combustion Institute,

May' 5-6, 1986.

15. Bittker, D.A.- and Scullin, V.J.: GCKP84--General Chemical Kinetics

Code lbr Gas Phase Flow and Batch Processes Including Hal Translier

Effects. NASA TP 23211, 1984.

16. Pratt. G.L.: Gas Kinetics. Wiley. London, 1969

17. Gear, C.W.: The Automatic Integration of Ordinary Differential

Equations. Commun. ACM. vol. 14. no. 3. 1971. pp. 176-179.

18. Zeleznik, F.J.: and McBride. B.J.: Modeling the Internal Combustion

Engine, NASA RP-I094, 1984.

19. Hindmarsh. A.C : GEAR: Ordinary Differential Equation System S_d,,er.

UCII)-30OOI, Re',. 3, Lawrence Li'.ermore Laboratory. 1974.

20. Liniger, W.; and Willoughby. RA.: Efficient Integration Meth_ds for

Stiff Systems of Ordinary Difli_rential Equations. SIAM J. Numer Anal.,

vol. 7, no. I, 197(I. pp. 47 66.

21. Brandon, I).M., Jr.: A Ne,,,, Single Step hnplicit Integration Algorithm

with A-Stability and Improved Accuracy,. Simulation. vol. 23. no. 1,

1974, pp. 17 29.

22. Babcock, P.D., Stutzman. L.F. and Brandon, D.M., Jr.: hnprovenlenls

in a Single Step Integration Algorithm. Sinmlation. vol. 33. no. I, t979.

pp. I - I0

23. Forsythe, G. E. : Malcolm, M. A. ; and Moler, C B. : Computer Methods

lot Mathematical Computations. Prentice-Hall, Englewood Cliffs. N J.

1977.

24. Shampine, L.F.: and Gordon. M.K.: Cmnputer Solulion of Ordmarv

Differential Equati_ms. The Initial Value Problem, W.H. Freeman and

Co., San Francisco, 1975.

25. Burcat. A.: and Radhakrishnan. K.: High Temperature Oxidation of

Propene. Combust. Flame, vol. 60. no. 2, 1985. pp. 157-169.

26. Radhakrishnan, K.; and Burcal. A. : Kinetics of the Ignition i_l Fuels in

Artificial Air Mixtures. I1: Oxidation of Propyne. Combust. Sci.

Technol., vol. 54. no. 1 6, 1987, pp. 85-101.

27. Forsylhc, G.E.: and Moler, C.B. : Computer Solution of Linear Algebraic

Systems. Prentice Hall. Englewood Cliffs. NJ. 1967.

28. Hindmarsh, A.C.: Construction of Mathematical Software. Part Ill:The

Control of Error in the GEAR Package for Ordinary Dilli:rential

Equations. UCID-30050, Part 3, Lawrence Li,,ermore l,aboratory,

1972.

29. Ortega. J.M.: and Rheinbolt, W.C: lteralive Solution of Nm_linear

Equations in Several Variables. Academic Press, Nc'* York, 1970.

30. Pratt, D.T.; and Radhakrishnan, K.: Fast Algorithms Ibr Combustion

Kinetic Calculations. Paper WSCI 81-46. Comhustitm Institute,

Oct. 20, 1981.

55


REPORT DOCUMENTATION PAGE I Form OMB Approved NO. 0704-0188

Public repot'-ing burden for lhlS COlleCtion o| inlormatlon is estimated to average 1 hour per response, including the t{me Ior rewewlrEj instructions, searching eXlshng data sources

gathenng and maintaining the data needed, and completing and rev=ew_ng the collection ol information Send comments regarding this burden estimate or any other aspect of this

colFeCtlon ol in|ormatJon, _ncludlng suggestions for reducing this burden, to Washington Headquarters Services, Directorale for mformahon Operations and Reports, 1215 Jefferson

Davis H_jhway, Suile 1204. Adlngton, VA 22202-4302, and to the Ofhce of Management and Budget, Paperwork Reduction Prolect (0704-0188), Washlnglon, DC 20503

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE I 3. REPORT TYPE AND DATES COVERED

January 1993 ] Technical Paper

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

A Critical Analysis of the Accuracy of Several Numerical Techniques

for Combustion Kinetic Rate Equations

6. AUTHOR(S)

Krishnan Radhadrishnan

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Lewis Research Center

Cleveland, Ohio 44135 - 3191

9. SPONSORING/MONITORING AGENCY NAMES(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, D.C. 20546- 0001

WU-505-31-42

8. PERFORMING ORGANIZATION

REPORT NUMBER

E-5861

10, SPONSORING/MONITORING

AGENCY REPORTNUMBER

NASA TP-3315

11. SUPPLEMENTARY NOTES

Krislman Radhakrishnan, National Research Council--NASA Research Associate at Lewis Research Center;

presently at Sverdrup Technology, Inc., Lewis Research Center Group, 2001 Aerospace Parkway, Brook Park,

Ohio 44142.

12a. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 07

12b.

DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

A detailed analysis of the accuracy of several techniques recently developed for integrating stiff ordinary differential

equations is presented. The techniques include two general-purpose codes EPISODE and LSODE developed for an

arbitrary system of ordinary differential equations, and three specialized codes CHEMEQ, CREK ID and GCKP84

developed specifically to solve chemical kinetic rate equations. The accuracy study is made by application of these

codes to two practical combustion kinetics problems. Both problems describe adiabatic, homogeneous, gas-phase

chemical reactions at constant pressure, and include all three combustion regimes: induction, heat release and

equilibration. To illustrate the error variation in the different combustion regimes the species are divided into three

types, reactants, intermediates and products, and error versus time plots are pre_nted for each species type and the

temperature. These plots show that CHEMEQ is the most accurate code during induction and early heat release.

During late heat release and equilibration, however, the other codes are more accurate. A single global quantity, a

mean integrated root-mean-square error, that measures the average error incurred in solving the complete problem is

used to compare the accuracy of the codes. Among the codes examined, LSODE is the most accurate for solving

chemical kinetics problems. It is also the most efficient code, in the sense that it requires the least computational

work to attain a specified accuracy level. An important finding is that use of the algebraic enthalpy conservation

equation to compute the temperature can be more accurate and efficient than integrating the temperature differential

equation.

14. SUBJECT TERMS

Stiffness; Reaction kinetics; Combustion chemistry; Numerical

method; Accuracy; Efficiency

17. SECURITY CLASSIFICATION "1"18.

OF REPORT

Unclassified J

NSN 7540-01-280-5500

15. NUMBER OFPAGES

6O

16. PRICE CODE

A04

SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRACT

OF THIS PAGE OF ABSTRACT

Unclassified Unclassified

Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-1B

298-102

NASA-Langley, ]993