2424
T. R.
Young and J. P. Boris
A Numerical Technique for Solving Stiff Ordinary Differential Equations Assoclated with the Chemical Kinetics of Reactive-Flow Problems T. R. Young” and J. P. Boris Plasma Dynamics Branch, Naval Research Laboratory, Washington, D.C. 20375 (Received April 26, 1977) Publication costs assisted by the Naval Research Laboratory
Among the many reactive-flow problems of interest are those associated with the disturbed and naturally occurring ionosphere, with combustion, with laser-beam-plasma interactions, and with many plasma physics and fluid-dynamic phenomena. An efficient and flexible numerical integration technique has been developed to calculate the time-dependent solutions of the stiff ordinary differential equations which help describe chemically reactive flows. A significant improvement in efficiency over the classical methods has been demonstrated using a second-order predictor, iterated-corrector scheme which involves an “asymptotic” treatment for equations which the technique determines are stiff. The method is self-starting (Le., does not depend on previously calculated values of the functions or derivatives). Since the error can be estimated, an appropriate time step is automatically chosen to give accurate results.
Introduction Mathematical descriptions of chemically reacting phenomena reduce to sets of simultaneous coupled ordinary differential equations. The solutions of these rate equations describe the time history of the various ionized or neutral species and temperatures of interest. Including the spatial dependence in reactive-flow problems often necessitates solving these sets of equations for many thousands of points in time and space. Classical methods for solving such sets of numerically stiff equations can be subject to very severe instabilities due to step size limitations which often render the solution of the problem impractical. Hence, a suitable numerical method that is stable for practical step sizes and has a high level of accuracy is sought. Rate equations for species densities and temperatures can normally be written as a set of first-order, coupled, nonlinear ordinary differential equations. Given a set of initial conditions and chemical reactions, the physical problem can be cast in the following form: dn,/dt = Q , - L,n, or
dn,/dt = Q , - (n,/.i.,) (1) where n, is the density of the ith species, Q, is the production rate, and L, (or 7,) is the inverse loss time (or the characteristic loss time). If Q, and T , are constants, eq 1can be solved analytically to yield
ni(t)=
Qi.i.i
+
[ni(0)- Qi.i.i]e-t’Ti
The characteristic loss rate (7i-l) describes how quickly the concentration of the ith species reaches its equilibrium value, &pi. Of course Qi and Li are normally functions of time. The differential equations often have time constants for the different species which differ by many orders of magnitude. The shorter time constants encountered are almost by definition much less than the time step required for performing numerical intergations in a reasonable amount of time. Moreover, instabilities result from the classical (polynomial) methods if the product of the loss rate Li and integration time step, 6 t , is greater than some number typically near unity (Le., 6tL; 5 1 for stability). The Journal of Physical Chemistry, Vo/. 8 1, No. 25, 1977
When stiffness occurs, the problem of integration becomes more difficult and nonclassical methods (i,e., free of the constraint 6tLi < 1) must be employed to obtain meaningful solutions within a reasonable amount of computer time. Nonclassical methods confront two obstacles, accuracy and stability. Nonclassical methods fall into two classes: those which treat all the differential equations identically and those which recognize the stiff equations or terms and integrate them by some very stable method while the rest of the equations are treated by a faster but less stable, classical method. This latter class of methods, known as “hybrid methods”, may gain large advantages over classical approaches by making use of the particular structure of eq 1. In hybrid methods two important problems must be addressed. The first is to establish a criterion determining by which method specific equations are to be integrated. The second is how the chemical balance and charge balance is to be maintained and how much imbalance will be allowed. Hybrid methods are fast, but conservation errors may occur after many time steps therefore care must be exercised when using them. The main ideas behind the classical methods along with a survey and discussion of the literature have been given by Byr0n.l He compares predictor, predictor-corrector, and Runge-Kutta methods and lists 75 references. Also, he gives a few comments regarding nonclassical methods. The specifics are far too numerous to give here. Among the other nonclassical methods which have been applied to reactive flow atmospheric calculations are the following: the hybrid method of Keneshea,2 the algorithm of T r e a n ~ rthe , ~ explicit method of Fowler and Warten,4 the selected order, multistep predictor-corrector method of Gear,5 the matrix method of Kregel.6 Solving complicated reactive-flow rate equations when coupled to a hydrodynamic calculation is a problem very different from detailed “one-point’’ chemical kinetic problems. First, the reaction-rate equations must be integrated for many grid points, perhaps tens of thousands, for relatively short periods of time between hydrodynamic time steps. Thus, self-starting techniques with low start-stop overhead are required. Second, the large size of the hydrodynamic grids usually precludes storing and auxiliary information for the spatial grid points between time steps. Since much of the computer time would be
Numerical Technique for Solving Stiff Differential Equations
2425
spent starting a high-order method, under these circumstances low-order methods invariably are n tore flexible and run faster. Third, hydrodynamic errors are rarely smaller than a few percent a t best, therefore application of detailed high-accuracy chemical kinetics is usually a waste of valuable time. In any case, reaction rates are at best poorly known so low-order methods again are the logical choice over complicated high-order approaches. Fourth, reactive flow codes generally get changed repeatedly so a general integration technique, rather than a complicated highly specialized program for a specific set of rates, is usually preferable. Finally, fifth, the quantity of computation involved in a multidimensional reactiveflow calculation is usually so vast when using even the very best methods that a viable technique should not contain computer evaluation of transcendental functions. Whenever analytic solutions are to be approximated a polynomial or rational-function approximation should be considered instead. Within this framework of large, multidimensional reactive-flow problem^,^-^ a special hybrid integrator called CHEMEQ has been written using a method termed selected asymptotic integration method developed by the authors. This method takes the above five special aspects of large reactive-flow problems into account. The next section describes the method and the last section describes a test problem which has been solved by several different methods and gives the conclusions.
The Asymptotic Integration Method Employed In CHEMEQ
6t
/~~> ( 0c2) stiff
< e2 normal
-
where c2 1. The equation is treated as stiff throughout the step, if it is considered stiff a t the start. The two classes of equations are then integrated by the same iterated predictor-corrector scheme but using a simple asymptotic formula to replace the usual second-order corrector equation in all those equations which were determined to be stiff by condition (6) above. The asymptotic integration method applied to the stiff equations best treats the situation where the solution is slowly varying or nearly asymptotic but the time constants are prohibitively small. This occurs when Q, and L,nrare large and nearly equal and is the result of coupling between the equations. Thus we treat the stiff equations with a very stable method which tends to damp out the small oscillations caused by the very small time constants. If, however, Q, and L,n, are small compared to n,, the simple classical methods can be utilized for these equations to give the combined method. The predictor part of the step is performed as follows:
ni(1) = ni(0)+ 6 tgi(0)
(normal)
and
integrates a set of coupled, stiff-differential equations of the form (1)which can be rewritten as CHEMEQ
d n i ( t ) / d t= g i ( n ( t ) ) (i = 1, . .., N ) where n is a vector with the N components ni.
takes place very rapidly after which much longer time steps can be taken. After a time step has been chosen, all of the N equations are separated into two classes, stiff and normal, according to the criterion:
n j ( l )=
n i ( O ) [ Z ~(0) i - 6 t ] + 26 t ~ i ( o ) Q i ( O ) (stiff) (7) 27i(O) + 6 t
(3) CHEMEQ
uses a one-step algorithm with no previous values of the variables n, or the derivatives gi stored from one cycle to the next. Furthermore, the method is based on a second-order predictor-corrector method which takes special notice of those particular equations which are determined to be stiff. The derivatives g, are evaluated from (4) or equivalently (1)as
(4) The initial trial step size is chosen independently of the stiffness criterion and is determined such that none of the variables will change by more than a prescribed amount which is typically of the order of a few percent. The actual 6t is chosen so that
6 t = c i min {ni(0)/gi(O)} (5) where 6 is a scale factor typically the same value as the convergence criterion described in eq 9 and where the minimum is taken over 1 Ii IN. This time step may be larger than some or all of the equilibration times in which case the corresponding equations would be called stiff. Nevertheless, when solved by asymptotic integration, this time step ensures that accuracy can be maintained. When a stiff equation is close to equilibrium, the changes in ni over the step will be small even though the adjustment rate toward equilibrium can be very much shorter than the time step. When the stiff equation is far from equilibrium, the time step should in fact be less than or comparable to the equilibration time t o ensure that the transition t o equilibrium of this stiff equation can be followed accurately. This readjustment, because of the very fast rate, generally
Here we are assuming that we start a t t = 0 and wish to find (ni(6t)j. In the following, we will let the integer in parentheses denote the iteration number. Thus, ni(m)is the mth iterated value of n;, an approximation to ni(6t). The zeroth iterate, n;(O),is the starting value at t = 0 and n,(l) is the result of the predictor step. Also note that gi(m) 3 gi(n(m)). The corrector formula for the two types of equations are
ni(m + 1 ) = ni(0)
+
and
+ Tj(0) - 6t]
ni(O)[Ti(m)
+ ~ j ( 0+) s t ] (stiff) By comparing ni(m + 1) with ni(m) for all of the N [Tj(m)
equations using the relative error criterion
Inj(m -t 1)- ni(m)l
min [ni(rn+ l), nr(rn)]
E3
(9)
convergence of each of the individual equations in each iteration can be ascertained. In practical applications e3 is typically The integration step is not complete until all of the equations, whether treated as stiff or normal, have satisfied eq 9 on the same iteration. In eq 9, it should be noted that n, is constrained by a minimum value in the denominator. Thus, when n, is decaying
-
The Journal of Physical Chemistry, Vol. 8 I , No. 25, 1977
T. R. Young and J. P. Boris
2426
exponentially toward zero, as it does for some species in many cases, we never divide by a number less than some lower bound. This number is chosen so small that it in no way affects the physics but eventually decouples the equation when n, gets small enough. Convergence for these equations is then trivial as should be expected. We have found that maximum speed is realized by keeping the maximum number of iterations very small; we typically use two to four. If convergence of all equations has not been obtained before or during the last iteration, the step is started over with a smaller time step. By keeping the maximum number of iterations small, a minimum amount of time is wasted on an unstable or nonconvergent step to find out that the iteration procedure has gone bad. By the same token, we have found it best to reduce the time step sharply (a factor of 2 or 3) when nonconvergence is encountered rather than to reduce it slowly. Less time is wasted this way getting down to a sufficiently small step for convergence if the initial estimated step size is found to be too large. On the other hand, when increasing the time step, as for example when convergence is achieved on the first or second iteration, we have found it best to only increase by 5-1070 each step. During the integration of several successive reaction rate steps, we use the appropriately modified time step from one converged integration cycle as the trial time step for the next integration cycle rather than using eq 5 . Once convergence of all the equations is achieved, the new values of the ni(St) are set equal to the values of n,(m 1). One can obtain convergence and completion of an integration step after only two derivative-function evaluations even when some or all of the equations are stiff. Thus the method is exceptionally fast. In fact, it is regularly utilized for massive chemically reactive flow calculations.7-9 All methods, such as our selected asymptotic integration method, which do not conserve chemical balance and charge balance automatically may be forced to do so by a t least two techniques. In one technique, balance can be restored by adding the various concentrations to find the errors and then by distributing these errors throughout the densities in a number-conserving manner. The major fault with this is that a portion of the errors is incorporated into concentrations from which the errors may not have arisen. The second and better method is to reduce the frequency of the asymptotic treatment to the point where errors due to nonconservation are within tolerable limits. Significant improvement in efficiency can still result. The complete CHEMEQ algorithm as described above is written in a simple transparent Fortran and has as an argument the name of the rate equation derivativefunction evaluator. It is intended as a very fast but moderate accuracy integrator which can be used a t each grid point of a large hydro- or magnetohydrodynamic code. When a single, more accurate, calculation is required, another method which is higher than second order may be superior. Many of the methods listed in the references combine slow speed with high accuracy. We have found the best method, as far as high accuracy is concerned, is a variable-order deferred-limit integrator of the type proposed by Bulirsch and StoerlO or Boris and Winsor.” The next and concluding section presents a comparison of CHEMEQ with other more or less standard algorithms. A direct comparison of methods written by several people (perhaps for different computers) is complicated by the fact that all methods may not be optimized equally in the coding or the layout or for the test machine. The version of CHEMEQ used in the following test problem and the
+
The Journal of Physical Chemistry, Vol. 8 1, No. 25, 1977
TABLE I : A List of the Seven Species Together with Their Initial and Accepted Final Concentrations for the Test Problem Number densities Initial
i 1 2 3 4 5 6 7
Species e0,CS+
cs cso, N, 0,
- Final
y;, cm-3
1.0 X 5.2 x 6.2 X 1.0 X 0 1.4 x 3.6 x
y i , cm-3
10, 10’ 10,
lo’,
1015 1014
4.9657897283 X l o 4 2.5913949444 x lo4 7.5571846728 X l o 4 1.5319405460 X l o 3 1.000 x lo1, 1.400 x 1015 3.590 x 1014
TABLE 11: A List of the Seven Reactions and Reaction Rates Through which the Seven Species of the Test Problem Interact Rate constant or No. Reaction frequency
a
1 2 3 4 5= 6
0’- t cst + cs t 0, Cs+ t e‘-. Cs t hv Cs t hv Cs+ t e0,- t hu 0, t e0 , t Cs + M CsO, t M 0, t e- t 0, -+ 0,- t 0,
7
0, t e- t N,
--).
-+
M = [Cs]
A
-
.+
0;
+ N,
[CsO,] t [ N , ] t
5 x cm3 1x cm3s-l 3.24 x 10-3 s-1 4 x 10-l s-l l x l o - ” cm6 s-l 1.24 x crn6 1 x l o - ” cm6 s-l
[O,].
function evaluator have not been close tuned or optimized in any way. Performance on a Test Problem a n d Conclusions The test problem involves the integration of seven rate equations which describe the time evolution of seven species of interest. This set of rate equations is stiff and not well suited for numerical integration by classical methods. The seven species together with their initial and accepted final (1000 s) concentrations are shown in (Table I). The set of chemical reactions through which these species interact along with their corresponding rate coefficients are shown in (Table 11). This problem was suggested by Kregel12(this particular set of equations was originally due to Edelson of Bell Laboratories)13who also provides a simple Fortran program which solves this problem and thereby provides a timing standard for comparison of different methods on different computers. Various methods reported at a numerical workshop14 are reported for comparison with the library CHEMEQ routine. Using the timing standard program, variations in computer speeds were eliminated by scaling. The vertical axis of the graph in Figure 1 gives running time in units of the standard benchmark. The variable 4 is an estimate of the error between the final calculated densities (nil and the accepted final densities (&I.
It appears from the slope on the graph (Figure 1)that Kregel’s method have about the same order of accuracy (second order), whereas Gear’s method is higher order, as indicated by the smaller slope. Since CHEMEQ and Kregel’s method have nearly the same slope, CHEMEQ is faster than Kregel’s method by a factor of 4 or 5 for all 4. Gear’s method is higher order than CHEMEQ, so CHEMEQ is more efficient than Gear for 4 greater than T o demonstrate the increased performance of CHEMEQ and
2427
Numerical Simulation of Complicated Flow Fields
t
IO 0 ,
I
I
I
I
I
I
TEST PROBLEM ~
Concerning accuracy, we feel that 4 less than lo4, in the current test problem, gives no useful additional information and could be a waste of valuable computer time in the context of a coupled chemical kinetic-hydrodynamic reactive flow problem. The hydrodynamics are not sufficiently accurate to deserve more than a 1% calculation of the reaction rates. This is reinforced, in any case, by rate coefficient uncertainties.
References and Notes
aFigure 1. A graph that gives a comparison of efficiency for various numerical integration schemes used to solve the test problem. The graph gives normalized run time vs. (the error criterion).
CHEMEQ over
classical methods a couple of points from a second-order predictor-corrector scheme without asymptotics has been plotted. CHEMEQ is approximately 25 times faster. There is clearly room for a substantial further improvement in CHEMEQ since the calculations were performed by our library integrator and no hand “tuning” has taken place. For example, improvement could be gained by performing the asymptotics on the equation(s) known to be stiff without testing.
P. R. Byron, Simulation, 11, 219 (1968). T. J. Keneshea, Research Report AF CRL-0221 (April, 1967). C. E. Treanor, Math. Comput., 20, 39 (1966). M. E. Fowler and R. N. Warten, I5M J . , 11, 537 (1967). C. W. Gear, Commun. ACM, 14, 176 (1971); DIFSUB for Solution of Ordinary Differential Equations (D2), bid., 14, 185 (1971). M. D. Kregel, J . Atmos. Terr. Phys., 34, 1601 (1972). J. P. Boris, 6. E. McDonald, T. P. Coffey, and T. R. Young, Proceedings of the DNA High Akitude Nuclear Effects Symposium, Vol. 2, DASIAC SR-130 (Dec, 1971). J. P. Borls, Proceedings of the 2nd European Conference on Computational Physics (also NRL Memorandum Report 3327 and references therein). A. W. Ali, T. Coffey, P. Kepple, S.Piacsek, A. Warn-Varnas, and T. R. Young, Trans., Am. Geophys. Union, 53, 458 (1972). R. Bulirsch and J. Stoer, Numer. Math., 8, 1 (1966). J. P. Boris and N. K. Winsor, MATT-652, Princeton Plasma Physics Laboratory Report (Nov, 1970). M. Kregel, private communication. D. Edelson, J. Chem. Ed., 52, 642 (1975). Proceedings of the DNA 1973 Atmospheric Effects Symposium (U), Vol. 4, DNA 3131P-9 (VI 75-442), pp 571-591.
Parallel Computation of Unsteady, Three-Dimensional, Chemically Reacting, Nonequilibrium Flow Using a Time-Split Finite-Volume Method on the Illiac I V Walter A. Reinhardt Ames Research Center, National Aeronautics and Space Administration, Moffett Field, California 94035 (Received May 12, 1977) Publication costs asslsted by the Ames Research Center, National Aeronautics and Space Administration
The system of unsteady, three-dimensional, partial differential equations used to simulate the inviscid flow of air in chemical nonequilibrium is approximated by a set of factored, finite-volume difference operators where the effect of chemical production is also contained in the factorization. The method is similar to that of Rizzi and Bailey, except for the emphasis on vector-matrix re-formulation designed to be suitable for the special architecture of modern advanced computers (e.g., the Illiac IV, CDC 7600, or Star). The systematic application of the operators yields a second-order accurate numerical algorithm. The method is programmed in the vector Fortran-like language called CFD developed in the CFD Branch at the NASA Ames Research Center. Translators are available for this language which allows debugging on serial computers, but all results for the examples given were obtained from the Illiac. The problem described is a numerical simulation of the flow in the high temperature stagnation region of a reentering space shuttle orbiter flying at large angles of attack (-40O). Capability for treating arbitrary geometry in a flow containing subsonic, sonic, and supersonic regions is demonstrated by this method. The air chemistry is described by a five-reaction model which includes the three dissociation reactions for Nz,02, and NO and the two rearrangement reactions involving NO. The vector-matrix formulation and the unique disk-memory mapping results in extremely efficient data management for the architecture of the Illiac and makes maximum use of the Illiac’s “data crunching” capability. Comparative running times are given for the Illiac IV and the CDC 7600.
Introduction The new generation of very fast, special purpose vector computers (e.g., Illiac Iv, CDC Star, CDC 7600, TI ASC, Gray 1, and IBM 370/195) has made possible the numerical simulation of complicated flow fields, including chemical reactions, about geometrically complex b0dies.l The need for these solutions results partly from the
continuing interest and usefulness of more sophisticated atmospheric entry vehicles such as the space shuttle. TO obtain such results, the split finite-volume method discussed in this paper is a viable numerical method. The equations that are approximated using this scheme are quite general and, with the exception of the easily modifiable chemical reaction model, are applicable for studies The Journal of Physical Chemistry, Vol. 8 1, No. 25, 1977
