To identify the fastest algorithm currently available for the numerical integration of chemical kinetic rate equations, several algorithms were examined. Findings to date are summarized. The algorithms examined include two general-purpose codes EPISODE and LSODE and three special-purpose (for chemical kinetic calculations) codes CHEMEQ, CRK1D, and GCKP84. In addition, an explicit Runge-Kutta-Merson differential equation solver (IMSL Routine DASCRU) is used to illustrate the problems associated with integrating chemical kinetic rate equations by a classical method. Algorithms were applied to two test problems drawn from combustion kinetics. These problems included all three combustion regimes: induction, heat release and equilibration. Variations of the temperature and species mole fraction are given with time for test problems 1 and 2, respectively. Both test problems were integrated over a time interval of 1 ms in order to obtain near-equilibration of all species and temperature. Of the codes examined in this study, only CREK1D and GCDP84 were written explicitly for integrating exothermic, non-isothermal combustion rate equations. These therefore have built-in procedures for calculating the temperature

Radhakrishnan, K.

English

NASA Technical Reports Server

_ FAST ALGORITHMS FOR COMBUSTION KINETICS CALCULATIONS: A COMPARISON*
Kr f slman Radha_n" f shnan**
NI_Sh Lewis Research Center
Many practical problems arising in chemically reacting flows require the
simultaneous numerical integration of larg_ sets of =hemical kinetic rate
equations of the type shown in figure I. The initial value problem is that of
finding the composition and temperature at the end of a prescribed time
interval, given the initial mixture composition and temperature, the pressure,
and the reaction mechanism. Multl-dimensional modeling of reactive flows
requires the integration of the system of ordinary differential equations
(ode's) given in figure 1 at several thousand grid points. To make such
calculations practicable, it is necessary to have a very fast batch chemistry
integrator.
Tn identify the fastest algorithm currently available for the numerical
integration of chemical kinetic ra%_ equations, several algorithms have been
examined. In _he present paper, we summarlze our findings to date -- details
are available in references (I) and (2). The algorithms examined in this work
include t, 9 general-purpose codes EPISODE and LSODE (refs. 3 and q), and three
special-purpose (for chemical kinetic calculations) codes CHEMEQ (ref. 5),
CREKID (refs. 6 and 7), and GCKPSq (refs. 8 and 9). In addition, an explicit
_mge-Kutta-Merson differential equation solver (ref. 10) (IMSL Routine DASCRD)
is used to illustrate the problems assJciated with integrating che "c_l kinetic
rate equations by a classlcal method. These methods are summarized in figure 3.
The algorithms summarized in figure 3 were applied to two test problems
drawn from combustion kinetics. These problems, summarized in figure 4,
included all three combustion regimes: induction, heat release ar"
equil_bratlon, Figures 5 and 6 presert variations of the temperature and
species mole fractions wlth time for test problems I and 2, respectively. Both
-- test problems were integrated over a %ime intec"al c _ I ms in order to obtain
: near-equilibration of all species and temperature.
%
Of the codes examined in this study, only C_EKID and GCKPSq were written
- explicitly for integratlrg e_othermic, non-isothermal combustion rate
equations. These therefore have built-ln procedures for calculating the
temperature (T). For the other codes, two different methods, labeled as Methods
A and B, were used to compute T. The following convention was adopted zn naming
these other codes: those using temperature method A were given the sufflx-A
(e.g. LSODE-A) and zhose using temperature method B were given the sufflx-B
(e.g. LSODE-B). In Method A, T was calculated from the mole numbers and the
_ Inltial mixture enthalpy using an algebraic energy conservation equation (glven
in figure 7) and a New_on-Raphson iteration technique. In this method, the
__ temperature is not an explicit independent variable, so the number of ode's is
Q
• Work partially funded by NASA Grantz NAG3-1q7 and N_G3-294.
•, NRC-NASA Research Associate: on le]ve from The University of Michlgan, Dept.
of Mechanical Engineering m_ Applied Mechanics, _mn _rbor, |II 48109.
257
_%__ |N_TIO_,[[_ g_ PRECEDING PAGE BLANK NOT FILME_
®,
1984012457-249
https://ntrs.nasa.gov/search.jsp?R=19840012486 2020-03-22T08:16:15+00:00Z
:0
equal to the number, NS, of distinct chemical species in the mixture. The
integrator therefore tracks only the solutions for the species mole numbers. In
; Hethod B, the temperature was treated as an additional independent variable and
evaluated by fntegratin Z its time-derivative given in figure 7. In this method,
: the number of ode's is equal to NS+I, and the tnLezrator tracks the solutt . ;
for both the temperature and the species mole numbers.
All codes were run on the NASA Lewis Research Center's IB)I 370/3033
; computer using single-preclslon accuracy, except GCKPS_ which was in double-
: precision. A typical computational run consisted of initializing the species
mole numbers, temperature, and CPU time. The integrator was then called with
" values for the necessary input parameters. On return from the integrator, the
¢o¢aI CFU tlme required ¢o solve the test problem was calculated. Other
performance parameters were also recorded -- see reference (2) for details.
Figures 8 and 9 present the computational work (expressed as the CPU tlme i
in seconds) plotted against the local error tolerance, EPS, for test problems 1
• and 2, respectively. For all codes except EPISODE, EPS is the local relatlve
error tolerance. For EPISODE, EPS is a mixed error tolerance ---relative fo_
species wlth initially nonzero mole numbers and for the temperature (method B)
and absolute for species with in_tially zero mole numbers. Also shown onl
figures 8 and 9 are the CPU times required by the explicit Runge-Kutta method
" for one value of EPS. Note the excessive CPO times required by ¢hls technique.
Its use would make multidimensional mode]ing of practical combustion devices
* prohibitively expensive.
For test problem I, very small values for EPS had to be used for EPISODE
(figure 8). For values of EPS 2 S x I0-', EPISODE predicted llt¢le or no change '_
In the composi_lon and temperature after an elapsed time of I ms. Similar
remarks apply to test proble_a 2 (figure 9), for which values of 10-" and 10-_
had ¢o be used for EPISODE-A and EPISODE-B, respectively. Although the runs
wlth EPISODE-B and EPS _ S x I0-_ were successfully completed, the solutions
(especially for minor species) were si_nlflcantly different from those given in
figure 6. With GCKPOq and _PS = I0", the solution for test problem i exhibited
serious instabillty and so was termlnated. A more detailed dlscusslon of the
accuracy of the codes tested in thls study can be found In reference (2).
Examination of figure 8 shows that the difference in computational work
required by methods _ and B is small for test problem 1, with method B being
more efficient. For test problem 2 (figure 9), the difference is small for
,_ large values of EPS. But for small values of EPS the difference is more marked,
with method k bein; stgnificanlty superior to method B.
Figures 8 and 9 show that LSODE and CREKID are superior ¢o the other
oodes. EPISODE is an attractive alternative, especially for _est probl,',_ 2. ,
However, in usinz EPISODE, a word of caution is in order. The compu%_tloPal
work tan be strongly dependent on the value for the initial steple,_th (HO) )
selected by the user. k poor _uess for HO can make EPISODE prohibitively i
expensive to use. Figure 10 illustrates this behavior fo_ test pr_.]_m 2. Note
In order of _nitude increase in the CPU time for a change in H¢ from 10 -_ to !'
10 "e s. hlthough not shown here, a poor guess for HO also resulted in ,
Inaccurate and unstable solutions. In addition, as discussed in refe, ence (2),
the error control performed by EPISODE is unsatisfactory for problems of the
£
1984012457-250
b{
g F
type eyamined in this study.
A simple method for incraasinz the efficiency of the algorithms as applied
to the present problem was explored. This involved updating the rate constants
kj and k_j (which was calculated from kj and the concentration equilibrium
constant) only for temperature chan_es greater than an amount AT. To avoid a
trial and error search for the optimum value of &T -- defined as that value
which results in minimum co,_putational work -- an approximation for it was
derived and is presented in figure 11. Comparlsons of figure 12 with figures 7
and 8 show the significant reductions in computational work realized by use of
" the above approxi_ition for AT.
References
1) Radhakrishnan, K., "A Comparison of the Efficiency of Numerical Hethods for
Inte_ratin Z Chemical Kinetic Rate Eqau_ions," NAS_ TH, 198q. Presented at '_
the 198t_ J_/_N_ Propulsion Conference, Specialists Session, New Orleans,
: LA, 7-9 February, 198zi.
2) Radhal:rlshnan, K., "A Comparison of Numerical Techniques for the
Integration of Stiff ODE's arlsin Z in Combustion Chemistry," NASa% TP, 198_.
J
3) Hindmarsh, A.C. and Byrne, G.D., "EPISODE: _n Effective Packa[Je for the
Integration of Systems of Ordinary Differential Equations," Lawrence i
Liwermore Laboratory Report UCID-30122, Rev. 1, 1977.
4) Hlnd_larsh, A.C., "ODEPACK, A Systematized Collection of ODE Solvers," .I
Lawrence Livermore Laboratory Report UCRL-88007, 1982.
1i
S) Youno", T.R. and Boris, J.P., "A Numerical Tecl_aique for Solving- St_ff
Ordinary Differential Equation,s Associated with the Chemical Kinetics of
Reactive-Flow Problems," J. Physical Chemistry, 81, 1977, pp. 2q24-2q27.
6) Pratt, D.T., "CPJ_K-ID: _ Computer Code for Transient, Gas-Phase Combustlon
Kinetlcs," Paper 83-21, Proceedin_ of _he Spring Technical Meeting of the
_estern States Section of the Combustion Institute, April 1983.
7) Pratt, D.T. and Radhakrlshnan, K., "CP.EK-1D: /%Computer Code for T_anslent,
" Gas-Phase Combustion Kinetics," to appear as a NASA Contractor Report, 198_
t/
,_ 8) Zeleznik, F.J. and McBride, B.J., "tfod(llln: the Internal Combustlon
Engine," NASA RPIO9g, 198_J.
9) Bittke_, D.A. and Scullln, V.J., "GCKP84-General Chemical Kinetlcs Code for
Gas-Phase Flow and Batch Processes Includin Z Heat Transfer," NASA TP, 1984.
10) Fox, L. (ed.), "Numerical Solutlon of Ordinary and Partial Differential
Equations," Addison-_esley Publishing Co., Inc., Readln_, MA, 1962.
11) Pratt, D.T., "New Computational Al_orithms for Chemical Kine_cics, '' National
Bureau of Standards Special Publica_ion $61, 1979, pp. 1265-1279.
.., 259
i
1984012457-251
• _ r 'i
w
: OF pOOR QUAL!'i_ /,
i Adiabatic, Constant-Pressure,Gas-Phase Chemical Reaction
I
D
)
dn.
; _- = f.(nk,T) ; i,k = I,NSdt ]
:' -I jJ
_: fi = " p i=IX(vi_- _i")(Rjj - R_j)
_: NS Vkj
' R. = k. II(Pnk)
J J k=l
:' NS v "
kj
R = k ]I(Pnk)
"J "J k:l
J,
! N.
"il kj. : A.Tj J exp ( -Ej/RT ) i
N !
: "J exp(-E jlRT) )k_j A_jT _
In the above equations, !
P
f. = molar rate of formation of species i per unit .lassof mixture,
i kmol-i/kg-mixture s
-' k., k = forward and reverse rate constants for reaction j
nJ_ "J = mole number of species i, kmol-i/kg.mixture
A'.j, A_j : pre-exponential constants in forward and reverse rate equations for
reaction j i
- Ej, E.j = activation energy in forward and reverse rate equations for reaction j, I
calImole
_' JJ = number of distinct elementary reactions in mechanism
NS : number of distinct species in gas mixture i
R : universal gas constant, 1.987 cal/mol K _
R.j = forward and reverse molar reaction rates per unit volume for reaction j, ' i
-" Rj, kmol/m3 s
: T = temperature, K I
p = mixture mass-density, kg/ms i! i!
vij _ij = stoichiometric coefficients of species i in reaction j as a reactant,
and as a product, respectively
Figure l Governing Ordinary Differential Equations
260
1984012457-252
,4
!
g r
"_" Given
/' Initial Hixture Composition and Tet....2rature
Pressure
. Reaction Hechani sm
Find, at the End of a Prescribed Time In'cerval
Hixture Composi¢ion and Temperature
- Figure 2 Problem Statement !
q
'I
He thod Description I
CC_LP84 Details not yet available.
.r,ZKtD Varlable-step, predictor-correcter method based on an
exponentlall_-fltted trapezoidal rule; includes filterln_ of
Ill-posed initial conditions and automatic selection of
Oacobi-Newton iteration or Near,on iteration.
l
L LSODE VariaLle-step, varlable-order bac|_ard-dlfferentlation i
EPISODE method with a generalized New_on iteration*. !
-_ CHEMEQ Vat lable-step, second-order predictor-correcter method wlth
. an asymptotic intezration formula for stiff equations. I[
DASCRU Variable-step, fourth-order, explicit Runge-Kutta-Merson I
r solver.
! *Other options are included in these packages.
If
": Figure 3 Summary of Hethods Studied
261
i
1984012457-253
ORIGINAL PAGE [g _'
OF POOR QUALITY
Two Problems Descrfbfn_ Adiabatic, Constant Pressure Chemical Reactions
Test Problem 1:
Combustion of a Mixture of 33 _ CO and 67 % H 2 with 100 % Theoretlcal Air
(taken from reference II)
12 Reactions
11 Species + Temperature
Pressure = i0 arm.
J
Initial Temperature = 1000 K
, Reaction Duration: 1 ms
1
! i
i Test Problem 2:
Combustion of a Stolchlometrlc Mixture of H and Air i
(taken from reference 9) 2
30 Reactions
15 Species + Temperature
Pressure = 2 arm.
Initial Mixture Temperature = 1500 K
Reaction Duration: 1 ms
Both Problems Include All Three Regimes of Combustion:
Induction, Heat Release and Equilibration.
u
Figure q Test Problems
262
1984012457-254
OR1O_,ltiEPP.,q_|9
OF POOR QUALITY
- 10-3___ 2500C, 10-3 3000
>" _ I " _ _
o ! / z _"
0 e_
_- 1o-4_ ,.= io-4#lli ! No/ / _ _ =E
7- itll" __,_o:
o ?500
- !!H '-_
,o-,_i/j// ,°-'j/ l,,. -106 -H 10"6
.20;'Ill II N I
I S/l;Co"2'o.ILL ./i/ i/ I
10-7 # " 1000 10-7 1_00
10-6 10-5 10-4 lO-] 106 10-5 10"4 10"3
t.S t,S
f_gureS - Variationwllh limeo! temperature Figure6 - Variationwithtimeof temperature
andspeciesmoletractionsIor testproblem1. andspeciesmolefractl_sfortestproblemZ. '
SolutiongeneraleclwithLSODE-BandEPS• SolutioogeneratedwithLSOOE-BandEPS•
105. I0-5.
263
- ®
1984012457-255
ORIGINAL PAGE Ig
OF POOR QUALITY|
: Hethod /_
t
_! For adiabatic, constant-pressure combustion reaction, energy conservation gives
!
l NST. n.h. = h = constant (1)
i=l 1 I o
t where, h i : molal-spectfic enthalpy of species i (J/kmol)I :
: and h 0 = mass-specific enthalpy of mixture (O/kg)
In this method, equation (1) was solved for the temperature using a Newton-
Raphson iteration technique.
1
T
lleth_d B ,_
Qt
Differentiation of eq. (1) uith respect to temperature (T) gives 1 !
r!
t
' NS
' rf.h.
dT i:l I I (2)
dt NS
n.c
i=l 1 Pl 1
where, is the constant-pressure specific heat of species I (J/kmol K). I
cP i
$
' ' I
In this method, the temperature was evaluated by integrating equation 121. 1
t
[
Figure 7 Evaluation of Temperature(for LSODE, EPISODE, and CREDO) , :
i
264 ®,
1984012457-256
._,;,, ° t_.
)-
h
.. OF POOP, QUALITY
05 ,4,'/ "=-- >
$ '3WI1ha::)
_, 265
,i
1984012457-257
ORIGINAL PAG_ |_ !
:* OF POOR QUALITY J
HO (s) CPU (s)
I
10-i 0.786
10-6 0.783
10-_ 0.791
: 10"* 7.91
" 10-! 8.0q
10-l° 0.772
Figure 10 Example of Effect of Initial Steplen_th (HO)
on Mork Required by EPISODE-_ (EPS = 10 "s)
for Test Problem 2
I
I
JLnapproximate expression for AT -- the maxlmumallowable _emperature
chan_e allowed before the reaction rate constants k. and _ are updated -- wasderived by requiring t at the maximumrelative error in resultan_ reaction
rates does not exceed the local relative error tolerance (EPS) required
of the numerlcal solution. The approximation for AT is given by
. AT = EPS.T (3) i
max + Nj ; + N !j RT RT -j ,
p
where, T is the current temperature, the bars] [denote absolute value, and I
the maximumis _aken over all forward and reverse reactions.
L
Figure 11 /Ipproxtmatton for AT
?66
....... T
1984012457-258
OR:C;: :
OF'Po0_" QUAL.j_1_
_lethod Test Problem
1 2
_8_ 0.i3_; z,73
CRE_ID 0 •23 1.O_
LSODE-_ 0•31" 0.$2"
LSODE-B 0.29* 0.SI*
EPISODE- k 0 •75* 0.54*
EPISODE-B 0.70" 0.67
CF_HEQ-_ 6.ql* 13.6"
Cf_HEQ-B 5.69* 12.3"
• meChod incorporated eq. (3)
Figure 12 Minimum CPU Time (in seconds on IBM 370/3033 computer)
Required for the Test Problems
267
1984012457-259


Fast algorithms for combustion kinetics calculations: A comparison

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server