Numerical solution of stiff boundary valued problems in kinetics and diffusion

Chemical Etwwmn~

Scmce

1979 Vol

34 pp 527-537

Pergamon Press

Pnstcd I” Great Bntmn

NUMERICAL SOLUTION OF STIFF BOUNDARY VALUED PROBLEMS IN KINETICS AND DIFFUSION J 0

L

WENDT,

Department

C

H

MARTINEZ,

of Chermcal Engmeenng, (Received

D

Umverslty

16 January

G

LILLEYt

of Anzona,

1978, accepted

and T Tucson,

L

CORLEYS

AZ 85721, U S A

14 July 1978)

Abstract-Two new numerical methods for the solution of stiff boundary valued or&nary differential equations are presented and compared The specific problem solved IS that of diffusion and reaction m a char pore, where eight species diffuse and react through ten free radical combustion reactions A new vanable mesh finite difference algonthm IS presented, and it IS shown how existing techniques for treating stiff mmal value problems can be adapted for the boundary valued case The numerical techniques presented have general application m problems mvolvmg dlffuslon and stiff chemical reactions, such as occur m combustion, gasdication, and reactor design MTRODUCTION

Many chemical engmeenng problems involve the solution of ordinary differential equations descnbmg the flow of a chemically reacting gas where the reactlon rate coefficients, and therefore charactensttc time constants, differ widely m magnitude Such systems are commonly referred to as “stdf”[l-31, and usually require spectal numerical methods for their solution Homogeneous combuslon reactions are known to be notonously stiff, and for many years It was common practice to involve the pseudo stationary state hypothesls[4], even though this approxunatlon can be mvahd [5] Development of computational methods to solve stiff nutlal value ordinary dlfferentlal equations has reached an advanced level[ i-3] However, the sltuatlon as regards to numerical procedures for solving stll? boundary valued problems IS not nearly as promismg[6,7], even though the problems occur frequently m combustion and reactor engmeenng A simple modlficahon of mttlal value techniques utthztng a shooting procedure loses its attraction when eight or more dependent vanables must be considered Two uutial value techniques deserve special mention, not because they are necessanly supenor for all stiff applications, but because they form the basis for the work reported here Traenor’s@] technique expands the reaction rate (of formation of species I) term, R,, m a Taylor Senes about the backward point, but makes the drastic approxlmatlon that the JacobIan, dRJt3Yj IS diagonally donunant and therefore solves for Yk+’ directly This approxlmatlon 1s very seldom vahd in kinetic systems, and so Traenor’s method, although computationally simple, 1s not among the most effective for stiff uutlal value problems It will be shown below, however, that when combined with successive underrelaxation, Traenor’s approxlmatlon IS useful m sotvmg stiff boundary valued problems tPresent address Department of Mechanical Engmeenng, Concordla University, Montreal, Canada SPresent address Energy and Envuonmental Research, Santa Ana Cahforma, U S A

*

527

A very efficient algonthm for solving stiff kinetic initial value problems 1s due to Tyson[9] and has since found widespread apphcatlon m theoretical chemical kinetic pre&ctlons[lO, 111 Tyson’s method evaluates the Jacobian iU?,/aY, at the backward step, k, and solves for Y:+’ through a full matnx mverslon at each step Tyson’s method 1s therefore not fully imphclt m the numerical sense, since no iteration IS mvolved, but IS especially powerful when combmed with a capablhty of rapid step size changes In this paper we focus on the stiff boundary valued problem of diffusion and reaction of I species m a smgle char particle pore This problem exhlblts many of the stiff charactenstlcs to be found m other problems such as a one dimensional lammar flame or reactor with axial diffusion We propose to solve this using an unequally spaced (contracting) gnd, with N gnd points, where the smallest gnd spaces occur where species gradients are steepest For the pore model, this LS at the pore mouth, for a lammar flame it would be at the flame-holder The finite difference scheme IS then solved using first, mversion of I tndlagonal matrices each of rank N (based on Traenor’s approxlmatlon) and second, mverslon of an N x N block tndlagonal matnx, where each block IS an I x Z matnx (based on Tyson’s approxlmatlon) Conare times computer vergence effectiveness and compared Special precautions are noted and guidelines for treatment of progressively sttffer systems are presented It should be noted that we do not claim to present an algorithm suitable for all stiff two point boundary valued problems, but rather only for those problems where stiffness arises through widely separated magmtudes of reactton rate coefficients Computational techniques for solution of other stiff situations, such as boundary layer type problems, have been presented elsewhere[6,12] Furthermore, we approach the problem from a pragmatic, applications onented vlewpomt, and therefore leave the mathematical convergence cntena to others it IS hoped, that our approach and Nonetheless, expenences will prove usefuI to other chemical engineers Identical) when faced with slmdar (but not necessarily

J 0

528 stiff boundary reaction

valued

problems

STATEMENT

mvolving

diffusion

L and

OF THE PROBLEM

We consider here a system of ordinary differential equations resulting from diffusion and reaction m a (char) pore -D,,

$i=

, Yz)

RAY,, Yzr l’s

For simplicity we have assumed, pseudo at constant temperature The boundary Fig 1) are

c.$*

binary diffusion conditions (see

x=0

x=L

y,=yp

(3)

The reaction rate of formatton of species I, R,, composed of both heterogeneous and homogeneous (combustion) reaction expressions, is generally a nonlinear function of all the concentrations Table 1 shows the kinetic scheme and reaction parameters utilized for this problem Reactions l-7 are the conventional, elementary reactions for the combustion of moist carbon monoxide[13] At temperatures greater than 1400°K reactions (2)-(S) are very fast, and indeed are sometimes partially equihbnated These are the chain branching reactions Reactions (5) and (7), which are the dissociation and recombination reactions are slower, as are the surface reactions[l4] @j(lO) The system of eight nonlinear coupled differential equations can therefore be considered to be stiff Previous experience with the shooting technique for only three non-stiff equations describing a similar problem[l5] indicated that the shooting technique would not work Therefore a finite difference method was formulated, leading to a set of algebraic equations that can be solved by matrix mechanics, provided the nonlinear terms are treated as known source terms Since the R, are very nonlinear in the Y,, a simple approach is to calculate these terms explicitly in terms of known Zero gradients ////l/l

///f////J

WEND-F ef al

values of the Y,, and then to iterate to convergence This approach, predictably, did not work because of instabilities arising from the stiff kinetics Following Tyson[9], a reasonable approach IS to expand R, about the backward iteration, and use

R:+,$$Y;“- Y,“)

R,*+’ =

where k denotes the kth iteration The Jacobian aR,/aY, can be calculated analytically for each iteration, and at each mesh point FINiTE DlFFERE NCE

Consider the unequally spaced or collapsing grid structure shown in Fig 1 In the case examined here It was found essential to have unequally spaced grid points, since gradients were not uniform over the region of Interest If n denotes the gnd point in question, AX, and AXE are the intervals to the west and east of n, where AXE = f AX, and f < 1 For example, with 9 mtervals and f = 0 8 the last interval is only 17% of the first interval m size Choice of the gnd adJustment factor is somewhat arbitrary, but should be such that sufficiently small intervals are used in regions where gradients are steepest The finite difference form of eqn (1) along this gnd 1s given by

-24, yk+* I”+, + [ AX,(AX, + AX,) 1 =&+,&I’

Heterogeneous

React Ion Concentrotlons, 1’1, I

I

I X=L lntervols

(Y;,“-

-2L.L A * n = AX,(AX, + AXE)

//f//////////////I/Speclfled Axaal Dlff uston ond Homogeneous React eon 111/f///// /,,,,,,,,,

N-2

n

Y;,,

(5)

where values at the k + 1 Iteration, which are unknown, are to be calculated from values at the kth iteration, at all mesh points, n - 1, n, n + 1 (n = 3 through N - 1) and for all species I = 1 through I Defining

x=0

w

FORMULATION

----I

ax&-+--dxE Rg I Pore model Physical and numerical schemes

Yy

(54

Numencai

solution of sttif boundary valued problems in klnetlcs and drffusron Table 1 Parameters

Chemical

Species

02

co

co2

H20 (4)

(1) (2) (3) Chemical

Reactions

Forward k

used for pore model OH

H 0 Rz (5) (6) (7)(S)

Rate Coefficient

Reverse

(gm/mole.cm.sec)

Rate Coefficzent

k (gm-mole,cm,sec)

1) CO + OH - CO2 + H

5 6x1011Toexp(-l.080/RT)

7 29x1013Toexp(-23.410/RT)

2)H+02=OHtO

2 24x1014Toexp(-16.800/RT)

1 3x1013To

3) 0 + H2 = OH + H

1 7~10~~T~exp(-9.450/RT)

7 3~10~~T~exp(-7300/RT)

4) 0 + II20 = OH + OH

5 57x1013Toexp(-18,000/RT)

5 75x1012Toexp(-760/RT)

5) HP + OH = Hz0 + H

2 19x101'Toexp(-5150/RT)

8 41x1013Toexp(-20.100/RT)

6) O2 + M = 0 + 0 + M

51 19 101aT-1exp(-118.700/RT)

25 63~10~~T-~exp(-34O/RT)

7) Ii20 + M = H+Oli+X

3 0x1015Toexp(-105.000/RT)

15 756~10~~T~exp(14628/RT)

8) Cg + 402 -f CO

376 1 T1 exp(-40,00O/RT)

0

9) cs + co* + 2co

3 6x10' To exp(-85,30O/RT)

0

10) cs + Hz0 + n2+co

1 2x10' To exp(-80,00O/RT)

0

Reactions Reactions

1) through

5)

6) and 7) are bimolecular

Reaction

8) has forward

expression

Reaction

9) has forward

expression

forward - 21r

and ternmlecular

I(1 + asY1) 3 L r = pore radius a8 = 16 4T - 2/r

RgfY3/(1

r. - pore radius.

10) has forward

expression

2 = ;

r = pore

2DrPn l&n = AXWAX, -24, c’ n = AX,(AX, +

(Sb)

+ a 9 Y 2 + b9Y3) 1 a9 = 3 28x10 -4Texp(40300/RT) = 2 59xT exp(-6100/RT)

I k10fY4/(1 f a10Y5 + bloYo)7 L J b9

r-

Reaction

reverse

KSfYl

[

radius.

al0 = 8 2x10 -l"T exp(66700/RT) -14 b = 1 74x10 Texp(82500/RT) 10

a”d

w

AXEI

B, 2 =

A,*=0

(74

2Dtm/@~~)*

G’b)

C, z = -~L&,/(AXE)~

(54

and rearrangmg,

529

(7c)

The other boundary condltlon, eqn (3) IS mtroduced settmg Y, N = Y,” and denvmg for n = N - 1

eqn (5) becomes

by

A In Y*+! rn I +B ,#IY:;‘+C,,Y:;:,. =R:,-j&%Y:.

(6)

Equation (6) holds for all values of I between 1 and I and for all values of II between 3 and N - 1 When n = 2 (X = 0). we apply the zero gradlent boundary condltton, eqn (2), to yield Y, 1 = Y, 3 and AX, = AX_= such that eqn

Y:;‘+

(8)

METHODS OF SOLUTION

Method I Dagonal

Jacobran-Truenor’s

approxrmatlon

Foliowmg Traenor, one n&t make the approxlmatlon that the Jacobian & 1s diagonally dommant, 1 e

(6) becomes (3)

=-C,N-,YP-R:N-I-,~,B:,N-,Y:N--I

(;Z:)

Y::‘-$$Y:: = R:,-,$,

&=O i%~Y:z

(7)

Substltutmg

eqns

for

r#~

(9) and (6) one obtams

(9) for the k + 1

530

J 0

a set of I x (N - 2) uncoupled

lteratlon form

equations

A,,y~,‘l,+S,,Y:=‘+C,,Y:::,=F;, 1= 1,1

L

WENDT et al

of the

(10)

k+L

y

I _______-

where

Y

B;,=B:,-J3:n I?

n

=

(104

R!‘, - #La,Y:,

These form a set of I tndlagonal form

:+:,

(lob)

matnx equations

of the

where the off diagonal block elements are each diagonal matrices and the diagonal block elements are each a full matnx given by 1=

1,I

& =&-/3nk

(11)

which can be solved usmg a standard tndlagonal matnx algonthm (TDMA)[16] Clearly an mltlal guess of the profiles of each species IS reqmred, and we found It adequate to assume mltlal constant concentrations with a fairly high value of 0 atoms to allow combustion to proceed Traenor’s approxlmatlon (eqn 9) accomphshes two thmgs first, it attempts to “predict” (albert poorly) the stltT reaction rate term, R,, at the (k + 1)th iteration from known values at the kth lteratlon. second, smce the & are always negative, It makes the diagonal terms m the tndlagonal matnx more posltlve, thus mcreasmg convergence stablhty This second facet IS not present m mltlal value problems It IS therefore hkely that Traenor’s approxlmatlon IS more useful for two point boundary valued problems than for mltlal value problems

fi, IS an I element

(14)

vector given by &=Rnk-&YI

(15)

Direct solution of this block tndlagonal matnx equation ts possible usmg a matnx analog of the TDMA This would entail, however, mverslon of a full I x I matnx N - 2 times per Iteration, and consequently a substantial amount of computer time To overcome this problem we used a different recursIon formula, also based on Gaussian Ehmmatlon, but where YN-~, YN--3, ,Yz were expressed m terms of YN_, thus Yn=Qn+,+&‘n+,Yt.-,

n=N-2.

where Q.. and en can be evaluated relatlonshlp

,2

(16)

from the recurslon

QN=O Method

II

Full Jacobran-Tyson’s

approxlmatton

Equation (6) can be wntten m the form of a single matnx equation, whrch can be partItIoned to give the followmg form t

QN--I =

A&‘-I@N---

Qn = A.-‘[&

~N-IYN~]

- CnQn+z-

(17)

&Qn+,J

PN =I -&kYnk+'

&--I

= R,* -&‘Ynk

= -AN’_&._,

P* = -A”-‘[B,P,+,

(18) - GP”+*l

(12) where A,,, & and G are I x I diagonal matrices,, Y, and R, are vectors, and _@ IS the full Jacoblan Further rearrangement produces a block tndlagonal matnx equation of the form

tMatrlces are denoted

are denoted by a single underhne by a bold character thus Y

thus

4

Vectors

Note that evaluation of Qn and & reqmred only the mverslon of diagonal matrices and simple matnx multiphcatlon The full matnx 8 needs to be mverted (usmg Gauss Jordan Reduction) only once to evaluate Y N-,

Machme preclslon equivalent figures IS reqmred Furthermore,

to slxteen slgmficant m place of the recurslon

Numerical

solution of stuff boundary valued problems m krnetxs and dlffuslon

relatlonshlps (16) Gausstan Ehmmatlon can also be applied, row by row, startmg from the bottom and usmg eqn (19) to calculate YN-_I We obtamed the best results (for our problem) by utlhzmg Gaussian Ehmmatton from the bottom row up, rather than from the top row down, that IS, workmg from the pore mouth where gradtents are steep to the pore end Furthermore, tncreased where gradients are zero stiffness requires more slgmficant figures to be camed, thus making this method somewhat machme dependent Best results for very stiff cases were obtamed using double preclslon on a CDC Cyber 175, Le 32 significant figures

In this section we first present numerical soluttons descnbmg species profiles for the pore model after convergence was obtamed Then comparison of methods to obtam these converged solutions IS presented Figure 2 shows the converged solution for species profiles at 1500°K Note that the free radical concentrattons are very much larger than equlhbnum, even though oxygen has penetrated to the pore end, thus showing both reactton and diffusion effects to be Important Table 2 shows the Jacoblan &, calculated at the pore end for the converged solution The largest real 1s negative etgenvaiue 1s -0 207 x lo’, the smallest -0 1520 X lo3 mdlcatmg a high measure of stdfness Furthermore, the off-dtagonal elements m the Jacobian are important and tt would seem to be a poor approxtmatton that the pti 1s diagonally dommant Figure 3 shows the converged solutton for species profiles at 2000°K Here, gradients are very steep at the pore mouth, mdlcatmg the need for a fine mesh in this region Furthermore here the largest real negative elgenvalue IS -0 332 x 10’ while the smallest IS -0 289 indicating that the system has even stiffer charactenstlcs than at 1500°K Table 2 shows, agam, that the diagonal elements do not dommate Discusston of methods to obtam the results shown in

531

IO -

-

/Ii20 x I02

00

004

002

006

0 0

008

Distance cm

Fig

2 Species profiles at UOO”K-zonverged solution radtus = 0 01 cm, pore length = 0 1 cm)

(pore

Figs 2 and 3 1s therefore pertinent to discussion of methods of solution for progressively stiffer systems Both solutions shown m Figs 2 and 3 were obtamed from uutud guesses assuming constant profiles for Hz0 and O2 equal to the free stream value and by setting the n-utlal mole fraction of 0 atoms to 0 01 everywhere The latter was necessary to allow for rapid lgmtion of the homogeneous reacttons in the pore volume Smce the concentration of any species at the pore end (mstde) was not known, and could vary over many orders of magnitude dunng iteratton to iteration, we report hereafter the 0 atom concentration m gmmoles/cm3 at that location as a function of iteration number as bemg representative of the stablhty and speed of convergence of the methods used In all cases exammed, all other species profiles had converged to a solution when the 0 atom concentration at the pore end had converged

Table 2 Jacobeans, &. at 1500°K and 2000°K at the pore end Jacobian

at

X=(CM)

1500% 0 -7 1 0 0 0 7 -7 7

Jacobian X =

at (CM)

8301D+O3 0 56501)+02-l 0364D+04 1 0364D+04-2 0 0 7518D+03-1 0364D+04 7518D+03 1 0364D+04-2 7518D+O3 0

0 2

0 0 5 3014~~01-2 6086n-03-1 0 0 0 -5 1762D+O3 1 0346D+O5 0 9 62483+02-l L611~+05-2 2 7468D+02 9 3894D+03-8 0819D+04-9 7468D+02-9 6195D+92 1 2611D+05 0 -4 2137D+03-2 2645D+O4-1 7522~+02 7495D+O2

4 1239D+05-1 1860D+04 1 1860D+04-9 2 6585D+05-7 1105J3+04 3316D+05 4 4367D+05-1 6932D+O5

1424D+O6 3 9 1176D+03 0 1176D+03 0 1322D+04-1 9 2694D+O4-4 1 2061D+06-1 2488D+06 3 1 1592D+06-4

4643D+O5

0173D+05 9939D+O3 372OD+O5 5064D+05 53939+05

2000°K 0

-2 6 0

9562D+05 0 3488D+O3-2 0604D+O4 2 0604D+04-1

0 0

2 9245D+05-2 -2 9245D+05 2 9245D+O5

0 0

0604D+04 2 0604D+04-1 0

0 0 7 4 3488D+O2-4 2858~-03-4 0 0 4 0 -8 5114D+O4 7 8937D+05 7 0 4 529RD+04-3 8596D+O5-2 1 8040D+04 1 2149D+05-1 9279D+O5-2 8040~+04-4 7863D+04 3 8596D+05 1 0 -3 6816D+04-9 6586n+O4-2

0 1

8728D+O4 5384D+O4

9623D+O5-4 3147D+O5 6 2852D+05 0097D+05 2 lOOOD+05 0 0097D+05-2 lOOOD+05 0 6985D+05-1 6363D+05-1 8589D+05 0382D+05 2 2163D+05-8 lOOOD+04 5297D+06 7 4894D+05-1 7486D+O5 3976D+06-8 6493D+05 7 0864D+05 336OD+O5 4 8763D+O5-8 9629D+05

J 0 L

532

01stancecm Fig

profiles at 2OOO”K--converged solution radms = 0 01 cm, pore length= 0 I cm)

3 Species

(pore

of methods for pore model at 1500°K At 1500°K the kinetics were moderately sm, as shown m Table 2 Solution of this problem usmg Traenor’s approxlmatlon of a diagonal Jacoblan was not possible, and the method was unstable after the 8th lteratlon as shown m Fig 4 This was because the errors m assuming a diagonal Jacobian outweighed convergence benefits acquired through the increased weight imposed on the diagonal elements of the tndzagonal matnx m equation (11) Successive underrelaxation, S = 0 5, provided sufficient damping such that the instability was removed, but with the penalty that convergence was achieved very slowly at the 853rd iteration as shown m Fig 5 These computer runs were completed with 10 mesh spaces and a gnd adjustment factor, f = 0 8 Furthermore, m all runs reported here, concentrations were prohlblted from being negative Method I did not require more than 8 significant figures to be camed dunng computation It was thought that Method II, based on Tyson’s approxlmatlon, should m pnncrple be stable and converge more rapidly smce It mvolves a better approximation for the R. at the k + 1 Iteration However, m our case, when the uutlal guess IS very far from the true solutron the off diagonal terms of the Jacoblan & evaluated at the backward iteration, change radically m magnitude between each Iteration This appeared to cause numerical mstablhtles within the first three iterations We solved this problem by usmg only the diagonal elements of j3,i for the first three mteratlons and the full Jacobian Companson

WENDT ef al

thereafter Therefore, success was obtamed by using a combmatlon of Traenor’s and Tyson’s approxlmatlons such that the latter was employed at the end of the 3rd iteration, that 1s. only after the off dmgonal elements of @,, had acqmred “reasonable” values This combination of methods allowed convergence to be achieved m 9 iterations as agamst 853 iterations for underrelaxed use of Traenor’s approxlmatlon alone Comparison of the path to convergence of this combmatlon and of the underrelaxed diagonal Jacobian alone 1s shown m Fig 5 and m detail on Fig 6 These results are also for 10 mesh spaces and a gnd adlustment factor f= 0 8 Sixteen significant figures must be carned to avoid round off error mstablhtles (ansmg m the block tndlagonal matnx Inversion) This problem was also solved sufficiently accurately with half the number of mesh points, and wtth factor a constant mesh size that IS with a gnd adjustment of 1 In each case convergence was achieved m 9 lteratrons The lack of Influence of mesh spacing IS probably due to the fact that, as shown m Fig 2, species gradients were not too steep In the case above, the Jacobian was calculated at each mesh point for each iteration This can sometimes be cumbersome so we investigated what would happen if it were assumed that the Jacobian at the center mesh point was representative of all the &,‘s at all the other mesh points This problem converged m 23 Iterations (as shown on Fig 7) mth the first three, again, only utilizing the diagonal elements Table 3 shows a comparison of Method I, modified through use of underrelaxation and Method II, modified for the first three iterations Even though Method II requves much more CPU time per iteration than Method I, it IS still more efficient smce the number of iterations required 1s slgmficantly less Method II, however, requues a much more complicated computer program

of methods for pore model at 2000°K 2000°K the pore model exhibits a high degree of stiffness as discussed above Under these conditions our modified (underrelaxed) Method I, did appear to converge slowly at approximately 892 iterations as shown in Fig 8 Method II, however, was more temperamental and requued (1), at least eighteen mesh intervals (as agamst tune for the previous case) and a gnd adjustment factor of 0 8 (Equally spaced mesh pomts led to mstablhtles), (2), diagonal dominance in the Jacobian and underrelaxatlon for the first three iterations only, to allow convergence When properly guided m this way for the first few iterations, Method II then proceeded to converge within 17 Iterations (Fig 9) and to a better tolerance than did Method I m 892 Iterations However, Method II now required 32 dlmts to be camed dunng the computation, that IS, utlhzmg double precision on a CDC Cyber 175 machme Table 4 shows a comparison of the two methods Comparison At

CONCLUSIONS AND RECOMMENDATIONS

If computer computational

time is not a factor, or if only a few experiments are to be made, we recom-

I

2

lterohon

3

#

4

5

I

6

I

I 7

Converged _* solutlcil

0 lterotlon 7 Y,=l3XU5 mok froc : I 6

I

Rg 4 Oxygenatomconcentration at poreend-methodI, unrelaxed (T = 1500°K)

<

: 40-

)

1

; 50-

L 1 i

j 60-

,

70-

80 -

90

loo

24 -

26 -

28 -

30-

E

0

0”“’

2-

4-

6-

8-

IO-

IZ-

“b x

14-

‘6-

IE-

2 iI 6

z al

f 225 6 zo-

*E

*-Full

32 ..

34,-

367

160

240

’ “1

320

Iteration

400

#

480

1 ““““I

560

Dloganol Jacoblan (Relaxed)

640

720

800

’ ” Fg 5 Comparison of methodI, relaxedandmethodII, modtfied (T = 1500°K)

80

Jocoblan

1 880

534

J 0

L

WENDT el al

-0,

-w

dr-

-Lo

-u-l

--o

--N

_-

O0

Numerical

solution

of strff boundary

p3/

aw

6

valued

pus ad

problems

m kmetlcs

and ddTuslon

535

+o B_o~ x [o]

i g~3ja~o~

6 pus

aJcd 0 *_Ol X [Ol

0

J 0

536

L WENDT el al

Table 3 Companson of methods-DECIO

computer (pore model, T = 15OO”K)

Method I (Relaxed)

Number

Dee

of 10

Average Core

lteratlons

CPU

tune

CPU

to

to

convergence

convergence,

tune/lteratlon,

Requirements,

853

sets

sets

15 0

words

4

11

025

0

of Cyber

M1nlmum

lteratlons 175

No

of

CPU Mesh

to time

39

Cyber 175 computer (pore model, T = 2000°K)

convergence to

convergence

NOTATION

with elements

A, ,,,

Method II Dominant and 3 lteratlOn.5~

to

17 11

2

9

Acknowledgements-This work was supported III part by the U S Department of Energy under Contract EF-77-C-01-2650 and m part by the NatIonal Science Foundation under NSF-RANN Grant No AER 75-03964

(Sa) (I X I) matrices C, n respectively (Sb) (lOa)

(Dugonally Qelaxed

892 4

Intervals

mend using Method I, that IS Traenor’s approximation with underrelaxation This method IS easy to program and can usually be forced into a stable region that will eventually converge by utlhzmg the correct under relaxation factor Method II, modified as described above, however, IS much more efficient computattonally, and so, for multiple computationai expenments, It 1s to be preferred As problems in kmetics and diffusion become progressively stiffer the followmg precautions should be noted (l), the block tndiagonal matnx should be mverted using (block) Gaussian Ehmmatlon, and starting with the top or bottom row depending on which side represents the largest species gradients, or If thts cannot be determmed both mversion algonthms should be tned, (2), double or even triple precision may be required, with increasing number of digits required for increasing stiffness, and (3). increasing the number of mesh points may be requued for mcreasmgly steep species gradients, although this problem can be alleviated by using a collapsmg gnd as was done m this work Although we present results only for the specific chemical engmeenng problem of drffuslon and reaction m a cyhndncal pore we believe the computational methods described above, and the guidelines and precautions presented, are applicable to many other stiff two point boundary valued problems that occur in chemical engmeenng problems of kinetics and diffusion

see eqn diagonal i% “, see eqn see eqn

3

38~

rethod I (Relaxed)

Number

to

9 94

8K

Table 4 Companson of methods
CDC

Method II Dominant Iterations)

(D~aqonally

91

18

& C Dim t”

f

R *PI n

I

k

kf, kr L

IY P _n

Q" R, Rak S ;I Y”L

see eqn (14)-a full matnx rank I see eqn (SC) pseudo bmary dlffuslon coefficient of species I in nitrogen gnd adjustment factor see eqn (lob) see eqn (15) number of species iteration index forward and reverse reaction rate coefficient pore length mesh mdex number of mesh points an I X I matnx defined in recursion formula (18) an Z element vector defined m recursion formula (17) rate of formation of species I, gm moles/cm’ set rate of formation vector evaluated at mesh pomt R at the kth iteration underrelaxation factor distance coordmate concentration of species I, gm moles/cm” concentration vector evaluated at mesh point n, at the kth lteratlon

Greek symbols Jacobian ,?RJaY, evaluated at mesh point II, at the kth iteration Bnk Jacobian matnx evaluated at mesh point n at the kth iteration AXE interval to the right of mesh point n AX, interval to the left of mesh pomt n

S:,”

REFERENCES

Alken R C and Lapldus L , A I Ch E J 1974 20 368 t:; Scmfeld J H , Lapldus L and Hwang M , I&EC (Fundis) 19709266

Numerical

solution

of stiff boundary

[3] Gear C W , Numerical ZmfJai Value Problems m Ordinary Drfferentral Equatrons, Chap 9, II Prentxe Hall, Englewood Chffs, New Jersey 1971 [4] Curtlss C F and Huschfelder J 0, Proc Nat Acad Scr 1952 38 235 [5] Frazier G C and Wendt J 0 L , Chem Engng SCI 1%9 24 95 [6] Flaherty J E and O’Malley R E , Math Comp 1977 3166 [7] Wdloughby R (Ed ), St@ D$erentral Systems Plenum Press, New York 1976 [8] Traenor C E , Math Comp 1966 20 39 [9] Tyson T J , TRW Rep No 9840 6002 RUOO Sept 1964 [IO] Engleman V S , Envrron Proteclion Tech Ser Rep EPA600/2-76-003 Jan 1976

CES Vol

34 No

4-G

valued [Ill

problems

in kmetrcs

537

and &ffuslon

Heap M P , Tyson T J , Clchanowlcz J E , Gershman R , Kau C J , Martm G B and Lamer W S , 16th Symp (Znt ) Comb, p 535, The Combustion Institute, Pittsburgh, PA

[I21 dzson C E I Math Phys 1968 47 351 [13] Baulch D L , brysdale D D , Home D G and Evaluated Kmerrc Data for High Temperature Vols l-3 Butterworths, London 1972, ‘73, ‘77 [14] South I W and Tyler R J , Fuel 1972 51 312 [15] Wendt J 0 L and Schulze 0 E . A Z Ch 102 [16] Camahan B , Luther H A and Wdkes J Numerical Methods, p 4.41 Wdey, New York

Lloyd A C Reacttons,

E 3

,

1976 22

0, Apphed 1969