A Numerical Method for Chemical Kinetics ... - ACS Publications

Chemical Kinetics Modeling Based on the Taylor Series Expansion

simpler than, other formulations. Because of its simplicity it should be useful to those who wish to write their own chemical kinetics programs or to incorporate chemical kinetics into a reactive flow calculation. The basic equations are eq 6, 10, 11, and 12, together with a subroutine to calculate the aiand bi and their derivatives, and a linear equation solver to carry out the iterative solution of eq 10. For computers with small word length, roundoff error can be reduced by replacing yn+land ykn+' by the backward differences Aynfl = yn+' - yn and Aykn+'= ykn+' - yn. In this case eq 6 becomes

Downloaded by UNIV OF TEXAS SAN ANTONIO on September 13, 2015 | http://pubs.acs.org Publication Date: December 1, 1977 | doi: 10.1021/j100540a016

so that the extrapolation is done with backward differences. Then the extrapolants are recovered from the relations ykn+' = yn + Aykn+'and the rest of the calculation proceeds unchanged. References and Notes (1) C. F. Curtiss and J. 0. Hirschfelder, Proc. Natl. Acad. Sci., 36,735 (1952). (2) C. W. Gear, Inf. Process., 68 (1969).

2413 (3) R. K. Brayton, F. G. Gustavson, and G. D. Hachtel, Proc. I€€€, 60 (l), 98 (1972). (4) A. M. Winslow, "Numerical Integration of Chemical Rate Equations", Lawrence Livermore Laboratory Report UCRL-73911 (Jan 1972) and "A New Formubtionof Backward Differenced Linear Mutisteo Methods for Ordinary DifferentialEquations", Lawrence Livermore Laboratory Report UCIR-681 (Jan 1973). (5) T. Rubner-Petersen. "An Efficient Alaorithm Usina Backward Time-Scaled Differences for Solving stiff Differencal-Algebraic Systems", presented at the Seminar in Numerical Analysis, Royal Institute of Technology, Stockholm (Oct 10-12, 1973). (6) G. D. Byrne and A. C. Hindmarsh, ACM Trans. Math. Software, 1, 71 (1975). (7) W. M. G. Van Bokhoven, I€€€ Trans. Circuits Syst., 22, No. 2, 109 (1975). (8) F. T. Krogh, "Changing Stepsize in the Integration of Differential Equations using Modified Divided Differences", in Proceedings of the Conference on the Numerical Solution of Ordinary Differential Equations, University of Texas at Austin (Oct 1972). (9) To be submitted to the J. Comput. Phys. (IO) A. C. Hindmarsh, "GEAR: Ordinary Differential Equation System Solver", Lawrence Livermore Laboratory Report UCID-30001, Rev. 3 (Dec 1974). (1 1) B. R. Bowman, D. T. Pratt, and C. T. Crowe, Symp. ( I t . ) Combust., [Proc.], 14th, 824 (1973). (12) G. D. Byrne, A. C. Hindmarsh, K. R. Jackson, and H. G. Brown, "A Comparison of Two ODE Codes: GEAR and EPISODE", University of California, Lawrence Livermore Laboratory Report Preprint UCRL-79141, January 1977; to appear in Comput. Chem. Eng. (13) R. J. Field and R. M. Noyes, J . Chem. Phys., 60, 1877 (1974).

A Numerical Method for Chemical Kinetics Modeling Based on the Taylor Series Expansion John

P. Kennealy and William M. Moore*

Department of Chemistry and Eiochemistry, Utah State University, Logan, Utah 84322 and the Air Force Geophysics Laboratories (OPR), Bedford, Massachusetts 0 1730 (Received May 27, 1977) Publication costs assisted by the Air Force Geophysics Laboratory

The sets of non-linear, first-order differential equations which describe the time behavior of chemical systems are ordinarily integrated with numerical techniques that utilize the analytic form of only the first derivative of species concentrations. However, these differential equations are somewhat unique in that the analytic form of all higher order derivatives can also be generated quite readily. A numerical integration algorithm based on the Taylor series expansion has been developed which exploits the specific form and symmetries of the differential equations governing chemical kinetics to provide an ease and efficiency of operation not found with other techniques. The algorithm is self-starting and a local convergence test is applied to vary the step size through the integration. In addition, a stiff equation technique is folded into the solution in the form of perturbed derivatives for any component which, by transient pseudo-steady-state behavior, forces the integration to small step sizes inconsistent with the time constants of the reaction set.

Introduction Over the last decade or two the increased availability of good kinetic data and the rapid improvements in digital computer technology have made the numerical modeling of large scale.chemica1 systems a most fertile area of activity. For example, numerous models have been developed and are being applied in the study of atmospheric chemistry problems, ranging from local air pollution effects through stratospheric ozone depletion to ionospheric perturbations produced by aurora. During the same time period, many other disciplines have developed their own requirements for sophisticated numerical modeling of a broad spectrum of phenomenology and descriptive systems. Much of the spectrum of modeling requirements shares

* Address correspondence t o this author a t U t a h State University.

the same general concern, i.e., the numerical integration of sets of coupled differential equations. Consequently, these apparently common requirements have been a major incentive for a large body of work directed toward the development of better general numerical integration techniques. Today one sees a profusion of elegant techniques being applied to the large variety of modeling problems. Whereas not long ago the chemist concerned himself with a relatively few extrapolative and predictor-corrector methods (e.g., Eluer,' Runge-Kutta, Kutta-Merson,2 etc.), today he sees numerical modeling being increasingly dominated by rather complex, general purpose methods such as the Gear a l g ~ r i t h m . ~ , ~ These multistep implicit numerical integration methods, which usually achieve significant reductions in execution time for typical "benchmark" problems, require the construction, handling, and inversion of matrices; consequently, storage requirements for the matrix methods The Journal of Physical Chemistry, Voi. 81, No. 25, 1977


2414

J. P. Kennealy and W. M. Moore

scale as N2,and execution times are proportional to some linear combination of N2 and N3 (when N is the total number of differential equations, i.e., the number of time-dependent species). In addition, it is required that the user provide as input to the integrator not simply the reaction set, but some explicit form of the differential equation set.4 For both the occasional user of an integration method (e.g., a chemical kineticist needing to parametrically integrate a small set coupled equations in order to 'estimate a rate constant) and the large-scale modeler (whose model often incorporates several hundred differential equations), a compact, efficient, and accurate explicit method which can be easily implemented on medium-sized computers and time-sharing systems is highly desirable. In following sections, a method and algorithm will be described which exploits the specific form and symmetries of the differential equations governing chemical kinetics. Called KEMO, not only can it and has it been meaningfully used with an interpretive language on a small minicomputer with only 4K of 12-bit core memory, but it is also routinely used in a Fortran version to integrate large sets of rate equations for modeling a variety of atmospheric phenomena. The method to be described is based on the Taylor series expansion, and its efficiency of execution is largely brought about by taking advantage of the very specific nature of the differential equations describing chemical kinetic^.^-^ Thus, it is not a general differential equation integrator, but it is configured and optimized for chemical systems; its stability and efficiency of operation are demonstrably superior to the more traditional but outmoded methods still being used, and it also compares favorably with the more "modern" techniques developed in recent years.

The Method and Algorithm To illustrate the technique most clearly, the description will initially be restricted to sets of bimolecular reactions. In a later section the extension to third and mixed-order reaction sets will be discussed, as will the treatment of reaction sets which produce stiff equations. Consider a set of J bimolecular reactions involving I different species, X , , where i = 1, 2, . ., I. Let N, be the time-dependent concentration of the ith species; i.e., N, = [X,]. The j t h reaction, with rate constant h,, can be written as

Rj

=

kjNjiNj2

(3)

and the general differential equation describing the time behavior of the concentration of X, becomes J

J

d N i / d t = j =E:1 Yij E S i , j 3

~ijRj= j= 1

+ Si,j4-

yijhjNjlNj2

Si,jl-

(4)

Si,j~

where yi = 1if Xi is a product of the j t h reaction and yil = -1 if is a reactant. (It is worth noting here that the I X J dimensional matrix y,while useful in the description of the method, is not required in the actual coding of the algorithm.) Clearly, the nth time derivative of N, is just

ii

(5) which, using Leibnitz's theorem, becomes

and Bqn-lare the binomial expansion coefficients

Bq"-l =

( n - l)!

(7)

( n- q - l)!q!

Obviously, given these recursive relationships, at any point in time one can easily calculate exact values for as many time derivatives as desired, for all of the species concentrations in a problem set; Le., given all the Ni(tl) one can calculate all the Ni(l)(tl); then using all the Ni(tl) and Ni(l)(tl), one can calculate all the Ni(2)(tl), etc. The ready availability of exact values for all the time derivatives immediately invites application of the Taylor series for stepwise numerical integration of the set of differential equations. Integrating from time tl to t2with the Taylor series

N i ( t 2 ) =N i ( t l )+

1- '

2

n--1

dtn

t,

n!

and, applying the recursive relationships in eq 5 and 6 to eq 8 yields

k*

xj, t xj21,xjs + xj,

(1)

where jl, j2, j3, and j 4 are simply the specific values of the index i associated with the reactants and products of the j t h reaction. For example, for a set of three reactions among five species (i.e., I = 5 and J = 31, a table can be constructed which displays the values of jl, j 2 , j3, and j 4 for each reaction: j l j 2 j3 j4 k

+ x,-t,

j = 1

X,

j= 2

X, t

j= 3

X,

k

x,-t,

t X+

k

X,

+

X,

+ X,

1 2 3 5

X, t X,

3 4 1 2

X,

6

3 2 4

(2)

Utilization of j l , j2, j3, and j 4 is introduced here because it not only simplifies the description of the technique, but it also becomes useful in the implementation of the algorithm as a device which helps minimize computation overhead and storage requirements. Continuing with the description of the technique itself, the time-dependent rate of the j t h reaction can be written as The Journal of Physical Chernistty, Vol. 8 I, No. 25, 1977

(9)

A careful examination of each term in the summation over n recommends some regrouping of factors and the definition of a new set of variables, which will exploit the symmetry of eq 9 and thereby optimize the coding and execution of the algorithm. The general term in the expansion is

n-1

L'

q=o

(n - l)! dqNjl dn-'-lNj2 (n - q - l ) ! q ! d t q dtn-q-l

which, by regrouping the factorials and powers of A t , can be rewritten as

2415


and all else would remain as described in the previous section. This assumes implicitly that the concentration of the third body does not change enough over each interval A t to affect the solution significantly. The more precise treatment for the catalytic reaction is totally equivalent to the method for the noncatalytic termolecular reaction (i.e., where j 5 # j 6 ) . For these cases, eq 3 is extended and becomes

By defining new variables

Rj = kjNjlNjzNj5 which eventually leads to the three-body extension of eq 15, i.e.

eq 11 is simplified to

and now a J X 6 dimensional integer table is used for the accumulation of the Mc(n)from the


and the Taylor series in eq 8 is reduced to

N , ( t z )= N j ( t l )+

5

n=l

1

Convergence and Error Analysis for a Truncated Series

t,

This reformulation does much more than simply make the equations more compact. Equation 13 is now a generating function which will efficiently generate any term in the Taylor series expansion from all the preceding terms; moreover, it is self-starting because Mi(o) N,. In addition, explicit evaluations of and multiplications by the various factorials and powers of At have been eliminated, and the total number of multiplications and divisions required to construct each series expansion has been reduced. The final, and perhaps most obvious, saving of computation time is best effected and demonstrated by the introduction of one final new variable - n-1 S j ( " - l )=

kj C q

=o

Mj1(q))Mj2(n-q-')

(15)

which is analogous to the Rj(n-l) of the original formulation and permits the rewriting of the generating function of eq 13.

This eliminates some redundancy which would otherwise result due to recalculation of the Sl(n-l) for each reactant and product of the j t h reaction. By referring to the J X 4 dimensional integer table, exemplified by eq 2, numerical values for the Mc(n)are accumulated by adding (or subtracting) the values for to the proper locations in the M array as the S array is evaluated. The algorithm meticulously generates the correct numerical coefficients for the stoichiometry, since the reaction set is scanned term by term, and the complete procedure is repeated in steps for each derivative.

Extension to Third-Order Reactions Retaining the notation of the previous section, the general termolecular reaction can be written Xjl t

xj* + xj5 xj, t xj4t +

Xj6

(16)

Three-body reactions can be divided into two principal classes, although the most general treatment will handle these two classes identically. In the first case, j 5 = j6, and the third body is just catalytic. For many problems, it is perfectly adequate to treat this by modifying the definition of kj (eq 12) so that -

k j ( 3 )= kjNjs A t

(17)

In theory, as already shown, any number of derivatives are available for the Taylor series expansion, and there are two approaches to form solutions. One method is to collect sufficient terms to produce a truncated solution that has an error less than a predetermined limit over a fixed step size At. However, with a system of differential equations such as found in chemical kinetics, the stringency of the solutions varies during the course of the reaction. The best method for these systems is to fix the number of terms in a truncated Taylor series, and then use a local convergence test to vary the step size and form solutions which satisfy the error/convergence test. As will be illustrated later, there is a relationship between the optimum length of the series and the accuracy desired in the solution. Due to the simple nature of the recursive relationships developed earlier, the KEMO algorithm can be and has been successfully operated as high as twelfth order to produce extraordinarily precise solutions. However, it is only very occasionally that this is required or justified. We follow the definition given by Gear* ''.. . the concept of convergence to mean that any desired degree of accuracy can be achieved for any problem, satisfying a Lipschitz condition by picking a small enough time step". The differential rate equations for chemical mechanisms seem to always satisfy a Lipschitz condition under the convergence test that compares the value of the last term in the truncated series to its sum as shown in (20) using the symbolism from eq 14, where the fractional convergence P

N j ( t l )t

n=1

value is V-l. As a reaction sequence proceeds, different chemical components limit the step size. Early in the reaction sequence, chemical products have convergence difficulties. This shifts to chemical intermediates and finally the reactants as their concentrations become small. In consequence, the total relative truncation error, e , is smaller than the convergence value V-l (Le., t V I 1). The relationship between the order of truncated series p , the local convergence, V-l, and the number of time steps, I , has been determined for many chemical mechanisms. In the absence of a quasi-steady-state period during the reaction sequence, the relationship shown in eq 21 holds

I = Vl'PCp

(21) The Journal of Physical Chemistry, Vol. 81, No. 25, 1977


2416

1

2ooc

P

I

i k.15 k = 1.0

I

A+A-+B A+C

The Taylor series solution will not be stable in this region, since a zero or near zero value for the first derivative does not guarantee small values for subsequent derivatives. The higher derivatives generally become relative large, which in turn demands a small time increment for convergence. A single example which illustrates the problem is the consecutive reaction sequence

x, k x, -k1, x, k,


1 e

I

IO

lo"

3

IO

IO

lo'

5

IO

Convergence Test V Figure 1. The relationship between the number of derivatives, the number of time steps, and the convergence test.

TABLE I: Comparison of Computer Time to the Number of Time Steps for a Given Number of Derivatives in the Taylor Series for the Chemical Mechanism A t A B and A --t C --f

Time increments No. of to constant deriv- truncation atives error

-

2

3 4 5 6

7 8 9

363 64 28 18 13

10 9 8

Computer relative time for operation

Relative computer time per time increment

3.65

1.00 0.59 0.49 0.44 0.42 0.44 0.46

-~

1.00 1.6 2.1 2.7 3.4 4.2 4.9 5.8

quite well, as shown in Figure 1 and Table I for the test mechanism used by deTar and deTarl and in a slightly different form by Gear? where the total time span = IAtav, and C p is a constant slightly dependent on the order of the truncated series. Such a relationship is not surprising, since it can be shown for many cases of ordinary differential equations that, when the Lipschitz condition is satisfied, the relative truncation error for a time span with a fixed At is given by9Jo

Through a suitable choice of hl, h2, and k3,the system can be made extremely stiff. Since there is an analytical solution, the true values of N2(t)and its time derivatives can be evaluated.ll A t any to in the stiff time regime, the instantaneous time derivatives will be large compared to the true time derivatives. The quasi-steady-state solution

N ~ ( s s=)Pi/Di

assumes that the first derivative for N , is exactly zero. Farrow and Edelson4 have discussed some of the pitfalls of using eq 25 in an uncontrolled numerical integration technique. Keneshea2 has corrected the difficulties of an uncontrolled quasi-steady-state approximation and our approach has been similar. At any time point in the stiff time regime, the instantaneous steady-state approximation is an extremely good indication of the concentration of the stiff component if all other components have accurate concentration values (i.e., the steady-state value represents good limiting value). The technique of partial integration employed by Adams and MegilP was adopted. It is based on the assumption that P, and D;are changing slowly and can be assumed constant over a time increment At. The differential equation dhTi/dt = Pj - DjNi

The Journal of Physical Chemistry, Vol. 8 1, No. 25, 1977

(26)

can be solved to give a time-dependent function for N,. Our modification is to equate this solution to a truncated Taylor series as shown in

N j ( t o )+ N i " ) ( t o ) A t = ( N j ( t ~) Pi Di ' I ) exp(-DiA t ) + Di A new first derivative is thereby synthesized from the old first derivative

N j ( ' ) ( t l=) N j ( ' ) ( t l ) ( -l e x p ( - D j A t ) ) / ( D i A t ) ( 2 8 ) new

Stiff Equation Modification All numerical methods for solving sets of ordinary differential equations encounter difficulties a t the onset of a stiff equation condition. In chemical kinetics, the problem arises when one or more of the chemical components achieve a quasi-steady-state condition. This happens when the sum of rate terms for production of the stiff chemical component Xi, symbolized by Pi, and the sum of the rate term for the removal of Xi, symbolized by DiNi, are large compared to the difference between them, the first derivative dNi/dt, as shown in

(25)

old

The new first derivative is used to evaluate all second derivative terms involving Xi which minimizes mass losses. In the limit, eq 28 generates a term which produces the quasi-steady-state value of eq 25. A stiffness test then can be constructed on the convergence test of eq 20. The quasi-steady-state value is on the left and the limiting value of the right side of eq 28 multiplied by At is on the right as shown in Pi - DiNi

Di

The test guarantees that the change in integration method will not be used until the stiffness problem is of the same magnitude as the convergence problem. The complete test


2417

TABLE 11: Cesium Problem with Initial and Final Conditions Chemical R e a c t i o n s React i o n s

No.

1 02-

-+

02

-+

E-

K = 4.00-01*(T)**

.OO*EXP

(0.000 / T )

2 CSf

E-

+

cs

K = 1.00-12*(T)**

.OO*EXP

(0.000 / T ‘

3 02

E-

+

92-

K

= 1.40-16*(T)“’

.OO*EXD

(0.000 / T )

4 92-

CS+

+

cs

f

02

K

= 5.00-08*(T)**

.OO*EXP

(0.009 /i)

+

CS+

f

E-

K

3.24-03*(T)**

5


Rate C o e f f i c i e n t s

cs

6 (32

cs

7 02 E 02

cs cs

9 32

cs

+

02

10 32

E-

t

02

cs cso2

+

cs

+

cs02

f

f

cso2

+

cso2

+

-+

cs92

+

cs02

+

02

02-

t

02

+

.OO*EXP ( 0 . 0 0 0 /T)

K = 1.00-31*(T)**

.OO*EXP

0.000 /T)

K = 1.00-31*(T)**

.OO*EXP

0.000 / T )

K = 1.40-16*(T)**

.OO*EXP

0.000 / T )

K

1.00-31*(T)**

.OO*EXP

0.000

K

= 1.24-30*(T)**

.OO*EXP

0.000 / T )

/TI

...................................................................................................... THERE H A V E BEEN 10 REACTIONS INPUT I N V O L V I N G 6 S P E C I E S .

Total Relative Error (E) f o r three derivatives w i t h w i t h local convergence a t V 500 With s t i f f equation No s t i f f equation modification modification

Species

Number Densities (cm-’) t = O t 1000

E-

1.0000 E + 02

4.9658 E

04

1.4 E

02-

5.2000 E + 02

2 . 5 9 1 4 E + 04

4.0 E

CSf

6.2000 E

f

02

7.5572 E + 04

CS

1.0000 E

f

12

1.5319 E + 03

cs02

0.0

02

3.6000 E

-t

14

f

-

8.0 E

04

1.0 E - 03

ni 1 4.0 E

9.1 E

-

-

04

05

-

ni 1

ni1

3.5900 E + 1 4

ni 1

ni 1

4750

4 10

14

1

Relative Computer Time demands that the component Xi becoming stiff is also the component limiting the size of the time increment by the regular convergence test. Major concentration components which are near equilibrium values will not satisfy both conditions. During the period of the integration when chemical components are stiff, mass and charge conservation are tested to prevent these errors from becoming propagated. Conservation of mass testing has been used successfully for algorithms which never conserve mass or charge.12 In the present algorithm, these problems have been minimized, as will be shown in the next section, and further improvement is still possible. The advantage of mixing integration methods is the dramatically reduced computation time for stiff equation conditions.

A Sample Reaction Set Edelson13originated a small 10-reaction problem, shown here in Table 11, which has an interesting practical origin;

04

1 . 4 E - 03

1.0000 E + 12

Number of Time Steps

04

it is a small, self-consistent model which illustrates the effects of releasing a modest amount of cesium vapor into the sunlit atmosphere a t an altitude of 70 km. Graphical output from the KEMO numerical integration of the 6 species, 10-reaction set shows that, neglecting the effects of diffusion, the local electron density in a small region of the atmosphere is enhanced for well over 1000 s by as much as eight orders of magnitude. Although not yet cited frequently in the literature, this has become an interesting and useful “benchmark” problem in parts of the chemical modeling community. Edelson’s cesium photoionization problem has a number of attributes which highly recommend it as a simple, first-order measure of performance for the integration algorithms which abound today. The first of these is its simplicity itself; it should require a rather trivial amount of effort to prepare and code the problem for treatment by any method, and the computation cost, except for the most inefficient or ineffective techniques, will be tolerably small. Second, the time-dependent species concentrations The Journal of Physical Chemistry, Vol. 81, No. 25, 1977



2410

span twelve orders of magnitude which should fairly well establish the immunity or vulnerability to propogation of round-off and other errors for most methods, and thereby also provide some measure of stability. Third, for approximately 0.1 < t < 100 s the pseudo-steady-state behavior of 02-, which is the minor species, makes the set of differential equations moderately stiff in that time interval; any algorithm which does not deal capably with stiff equations will spend a greatly disproportionate fraction of its total execution time integrating from t = 0.1 to t = 3 00 s. Fourth, achieving the complete and proper solution will demand that mass and charge conservation be well maintained throughout the solution. Edelson's problem has been treated with the KEMO algorithm on both a UNIVAC 1108, using a Fortran version, and on a 4K DEC PDP8/E, using a Focal version. (Focal is an interpretative language, similar to but generally less efficient than the better-known Basic, and only 1000-2000 of the 4K 12-bit words are actually available for a Focal program; the remainder are required for the resident Focal interpreter.) On the 1108, an integration of the reaction set which produced solutions with a maximum error of 0.14% (predicted 0.20%) required 410 time increments. Table I1 shows how easily the reaction set is alphanumerically coded for input to the Fortran version; in addition to the reaction set, of course, there is one line (card) of input for each species, which contains the alphanumeric formula of the species and its concentration a t t = 0. Different algorithms must be coded and run on the same machine under the same operating system in order for execution time to become a meaningful comparison criterion. The cesium problem described here was used in a 1973 workshop test of numerical integration algorithms on a CDC 6600 computer.14 We have been provided with the results of that test,15 and we have had access to the same computer under the same operating system. The KEMO algorithm performed a t least as well as any of the other algorithms tested, including a version of the Gear algorithm.16 For larger problems, it seems clear that KEMO would operationally outperform any algorithm requiring matrix multiplication and matrix inversion. For any specific order the number of machine operations for the Taylor series algorithm is proportional to the sum of the number of species and the number of reactions. The matrix methods require a number of operations proportional to the square of the number of time-dependent species for multiplication and proportional to the cube of the number of species for matrix inversions. For a large chemical reaction set containing over 100 species and reactions, we estimate the machine operations for the matrix inversion problem alone would equal the total number of operations for the KEMO algorithm. This estimate includes an allowance for the use of sparse matrix methods and a factor of 10 for the number of time steps per i n v e r ~ i o n . ~ It J ~ has been reported elsehwere that implicit algorithms such as the Gear are "prohibitively costly" due to their utilization of matrix meth0ds.l' Furthermore, for large problems the preparation of differential equations for input to the integrator becomes a not insignificant task and generally involves development and utilization of a rather specialized "differential equation writer" p r ~ g r a m which , ~ ~ , is ~ ~usually not portable from one type of computer to another.18

Summary of the Taylor Series Integration Algorithm The entire algorithm, as developed in the preceding sections, can now be totally described and implemented The Journal of Physical Chemistry, Vol. 81, No. 25, 1977

with five equations, which are repeated here in order of their usage:

-

kj

kjAt

M i ( ' ) ( t ) = Ni(t)

where At is such that the magnitude of the last term in the truncated series is less than a preselected fraction (usually 0.01 or 0.001) of Nifor all i. If, for a given value of At, a series fails to satisfy the convergence test, then At is reduced by a factor of 2; this is easily done by reducing the stored values for all of the Mi(n)by 2n and only resumming each series. Division by two, in this manner, is repeated until all series satisfy the test. Periodically (e.g., after every five successful steps), an attempt is also made to increase At by a factor of 2, which keeps the size of the time step very near optimum all the time. The experimental evidence from quite a few years of experience with this algorithm indicates that the method is absolutely stable with respect to truncation and round-off errors; i.e., the maximum fractional error in any species concentration at any time does not exceed the maximum fractional error permitted during any previous step. Without exception, the authors's experience has been that this Taylor series method is an easily implemented, efficient, and accurate numerical integration technique for both large and small sets of reactions.

Acknowledgment. This work has been supported in part by the Defense Nuclear Agency and the U.S. Air Force Geophysics Laboratories through Contract No. F1962872-C-0317. References and Notes D. F. deTar and C. E. deTar, J . Phys. Chem., 70, 3842 (1966). T. J. Keneshea, Air Force Cambridge Research Laboratory Report AFCRL-67-0221 (1967). C. W. Gear, Commun. Am. Comput. Mach., 14, 176 (1971). L. A. Farrow and D. Edelson, Int. J . Chem. Kinet., 6, 787 (1974). I. G. Darvey, S. J. Prokhovnik, and J. F. Williams, J . Theor. Biol., 13, 90 (1966). P. A. Adams and J. G. Sheppard, J . Chem. Soc., faraday Trans. 7 , 70, 130 (1974). W. M. Moore and T. K. Eccies, Air Force Cambridge Research Laboratory Report AFCRL-TR-73-0749 (1974). C. W. Gear, "Numerical Initial Value Problems in Ordinary Differential Equations", Prentice-Hall, Englewood Cliffs, N.J., 1971. F. Scheid, "Theory and Problems of Numerical Analysis", McGraw-Hill, New York, N.Y., 1968. P. K. Henrlci, "Elements of Numerical Analysis", Wiley, New York, N.Y., 1964. H. S. Johnston, "Gas Phase Reaction Rate Theoty", The Ronald Press, New York, N.Y., 1966. G. W. Adams and L. R. Megill, ESSA Technical Report, ERL 177-SDL15 (1970). D. Edelson, J . Chem. Ed., 52, 642 (1975). Numerical Methods Workshop (March 1973) at Visidyne Inc., Burllngton, Mass. T. C. Degges, private communication. A. C. Hindmarsh, "Gear: Ordinary Differential Equation Solver", UCID3000 1, Lawerence Livermore Laboratory, Universityof California, Livermore, Calif. (Aug 1972). R. J. Gelinas, J . Comput. Phys., 9, 222 (1972). D. Edelson, Comput. Chem., 1, 29 (1976). R. P. Dlckinson, Jr., and R. J. Gellnas, "Numerical Methods for

2419

Solution of Kinetic Rate Equations Differential Systems”, L. Lapidus and W. E. Schiesser, Ed., Academic Press, New York, N.Y., 1976.

Discussion


D. EDELSON (Bell Laboratories). The observation that the method violates conservation laws should be sufficient t o warn against its use. In a problem having a rising solution, the errors may well be swamped by the solution itself. In a falling solution, the results can be disastrous. In the test problem cited, such errors are known to lead to an unbalance of many orders of magnitude in charge neutrality. Ad hoc adjustment of the balance may appear to work, but it is really changing the problem parameters.

CAKE.

In cases having extreme sensitivity to particular rates or concentrations, e.g., switching or oscillating reactions, the results will be completely misleading. W. M. Moore. We agree with Dr. Edelson that uncontrolled use of a technique which does not conserve mass is a dangerous procedure. However, we have never encountered a situation where applying the conservation of mass as a boundary condition to a nonconservation integration technique has generated incorrect results. With KEMO, the stiff equation method can be turned off or on with a flag change, so that any mechanism can be tested initially to satisfy the user of the desirability of including the stiff equation solution. A reduction in computing time by a factor of 10 with only a small loss in relative error as shown in Table I1 of the paper can be very useful for multiple application to a particular chemical mechanism.

A Computer Program for the Solution of Kinetic Rate Equations P. F. Ridler Depadment of Physics, University of Rhodesia, P. 0. Box MP 167, Mount Pleasant, Salisbury, Rhodesia (Received May 27, 1977) Publication costs assisted by the University of Rhodesia

The paper describes some of the features of a computer program designed to solve differential rate equations. The input is user-oriented, free of format restrictions, and is interpretive, so that the input data may be checked and diagnostic messages printed for the more usual errors. The data are stored using list processing techniques, so that the information defining the reaction topology can be recovered efficiently. The Jacobian matrix, necessary for the solution of the nonlinear equations, is stored using a sparse matrix technique to ensure that a minimum number of arithmetic operations is necessary for the solution of the equation set. The integration method is a modification of the backward differentiation formula which is fast, accurate, and extremely robust. The range of problems on which the program has been tested is very comprehensive, covering many of the stiff problems quoted in the literature as well as some oscillating reactions which are extremely demanding. The program has also been subjected to student use for several years without damage. CAKE is written in standard Fortran and is fully portable.

group, vpi is the stoichiometric number of the p t h species in the ith group, s is the number of distinct chemical species in the entire set of reactions, and g is the number of groups of species in the reaction set. As xp may occur in more than one group, its rate of decay in the complete reaction set must be obtained by summing (1)over each group in the set, so that

1. Introduction The simulation of chemical reactions has received considerable attention in the last few years’-3 and is proving to be a valuable tool in the better understanding of many chemical processes. This paper describes some of the features of a simulation program which is essentially user oriented and which has proved, over a period of 5-6 years to be extremely rugged. It is suitable for small machines and requires only about 64K bytes of core when overlaid in three sections, or about 100K when not overlaid. It uses sparse matrix techniques for processing and, while not necessarily optimally fast is reasonably so. It has interpretive input routines and is thus able to produce error diagnostic messages for the commoner input data errors. It has, so far, resisted the attempts of all users to make it produce spurious results. The kinetic rate equations may be written in completely general form by considering the rate of decay of the p t h chemical species in the ith group of species, a “group” of species being the collection of chemical species on one side of a reaction. This rate of decay will be

This set of equations, for p = 1, . . ., s, must be solved to give the behavior of each species as a function of time. It is possible to formulate (2) explicity and then use any available differential equation solving package to obtain a solution, but this method involves considerable mathematical and programming effort on the part of the user and can give rise to misleading results. A better approach is to design a program which accepts, as input data, a description of the reaction scheme which is close to that used in normal practice and process this into whatever form is required for the final solution.

where x p is the concentration of the p t h species a t time t , hi; is the rate coefficient from the ith group to the j t h

2. The Input Language The input language of CAKE is a little different from some others, but has advantages in flexibility and compactness of notation. Edelson4 has chosen to require The Journal of Physical Chemistty, Vol. 81, No. 25, 1977

A Numerical Method for Chemical Kinetics ... - ACS Publications

Recommend Documents