Computer Applications to Chemical Engineering - American Chemical

11 A n Iterative Approach for the Solution of a Class of

Downloaded by UNIV OF CALIFORNIA SANTA BARBARA on April 3, 2018 | https://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch011

Stiff O D E Models of Reacting Polydispersed Particles FERHAN KAYIHAN Department of Chemical Engineering, Oregon State University, Corvallis, OR 97331 The mathematical models of the reacting polydispersed parti cles usually have stiff ordinary differential equations. Stiff ness arises from the effect of particle sizes on the thermal transients of the particles and from the strong temperature depen dence of the reactions like combustion and devolatilization. The computation time for the numerical solution using commercially available stiff ODE solvers may take excessive time for some sys tems. A model that uses Κ discrete size cuts and Ν gas-solid re actions will have K(N + 1) differential equations. As an alter native to the numerical solution of these equations an iterative finite difference method was developed and tested on the pyrolysis model of polydispersed coal particles in a transport reactor. The resulting 160 differential equations were solved in less than 30 seconds on a CDC Cyber 73. This is compared to more than 10 hours on the same machine using a commercially available stiff solver which is based on Gear's method. In the following sections some background information on stiff ordinary differential equations will be given and the gener al finite difference approximations for particle temperatures will be derived. Later, the technique will be applied to coal pyrolysis in a transport reactor where the difference equations for reaction kinetics will be discussed and the calculation results will be compared with those obtained by the previously established tech niques. BACKGROUND Mathematical models that contain ordinary differential equa tions face an inherent computational d i f f i c u l t y associated with the stiffness of the equations. Stiffness of ordinary differen t i a l equations depends on the relative magnitudes of the response modes or the characteristic time constants of the system being modeled. In solid fuel conversion problems where particles of varying sizes are considered the differential equations for the thermal transients of the particles are usually s t i f f . Estab0-8412-0549-3/80/47-124-217$05.00/0 © 1980 American Chemical Society Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.


218

COMPUTER APPLICATIONS TO CHEMICAL ENGINEERING

lished numerical techniques, like Gear's method, are available for the solution of these problems. However, this becomes verycostly for systems with large number of coupled differential equa tions. The present work illustrates an alternative approach through finite difference approximations which is shown to be ef fective for problems dealing with the reaction of polydispersed particles. After a brief description of the stiffness conditions in ordinary differential equations the approximate solution tech nique w i l l be presented in general terms. Then, a specific appli cation to coal devolatilization w i l l be discussed in detail. Consider a nonlinear model of η differential equations g

= f C y . t ) . y(t )

where y = (y y . . . ίΥ Α ) gives 2

0

- y

o

yn).

(l)

o

The linearization of this model around

0

dy dt

f ( y , t ) + f ( y , t )(y-y ) + f\ ( y , t ) (t-t ) o ο y ο ο ο t o ο ο

(2)

w

dy dt

(3)

g(t)

where J is the Jacobian at ( y , t ) and g(t) is only a function of the independent variable t. The local solution is 0

0

V

y(t) = e

0

J (t-x) e g(T)dx 0

(4)

t For a moment, assume that Eq. (3) is correct for a l l t where Eq. (4) gives the true solution. Then the response modes of the system are given by the eigenvalues ( λ ) of the Jacobian J . For a stable system, the eigenvalues have the property that 0

Re(A ) i

i=l,...,n

< 0

(5)

For the special case of distinct eigenvalues the solution is writ ten in terms of -t/τ.

-t/τ

where

= - Ι / λ for i = l , . . . , n are the system time constants. In a transient problem we are interested in the solution for t *w n b. At • [ra (t ) ( 1 - — - m (t a. 2à. J

)

Eq. (8) can be approximate-

-At /à. y

8

n

)(1

i +

1

J

At -At /a. - ) e ] 2a. (9) n

J

J

J

where the quantities a.(t y W

) = [a.(t ) + a.(t

ν

+

b

= f y v = £

W

J]/2

t

ji n- »

/ 2

1

+

W

i

2

are assumed to be constant for the interval A t . The corresponding expression for mj(t ) can be derived ac cording to the particular reaction. In general terms the rate of mass loss, -dmj/dt, is a function of particle radius r j , kinetic rate R(Tj, Tp, Ρ ) , and gas phase compositions y ' s . This can be n

n

Squires and Reklaitis; Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

11.

KAYIHAN

Iterative Approach

221

written as dm. dt

f(r.)

(10)

R (T., T , P) g (y 's) F

G

Eq. (10) can be approximately integrated through analytical techniques similar to Eq. (7) over the time increment A t either by assuming a constant average particle temperature between T j ( t ) and T j ( t i ) or by using a linear temperature profile in this range. For combustion reactions Levenspiel (4) gives the constant temperature integration for reaction and gas or ash diffusion controlled cases. The integration of the pyrolysis kinetics w i l l be demonstrated in the following section. Eq. (9) and the integrated version of Eq. (10) constitute the finite difference approximation for the model. The solution is generated through successive iterations on the complete temperature profiles and other variables. A convenient i n i t i a l guess can be generated by using dmj/dt = 0 which corresponds to the simple heating of the particles. The details of the computation procedure is described in the following example on coal pyrolysis. n

n


n

EXAMPLE:

PYROLYSIS OF POLYDISPERSED COAL PARTICLES

At low pressures the kinetics of thermal decomposition of pulverized coal particles can be represented by an i n f i n i t e l y large number of parallel f i r s t order reactions as shown by Howard, et a l . (5,(>). For those cases where data on species production rates are not available the devolatilization kinetics can be written in terms of one expression which combines the kinetics of a l l possible reactions by using a normally distributed activation energy function. However, i f data are available on species kinetics, then devolatilization can be characterized by a f i n i t e number of f i r s t order reactions. For example, Suuberg (7) gives the kinetic parameters for 15 parallel reactions for eight identifiable v o l a t i l e species in the pyrolysis of Montana l i g n i t e . This particular set of data was used in this study to examine the computational problems that arise in modeling of pulverized coal pyrolysis. For reactor design calculations i t is necessary to know the total devolatilization rate as well as the species production rates. Therefore, one needs to include in the reactor model a l l the reaction rates that are available for the devolatilization of the particular coal. Kayihan and Reklaitis (8) show that the kinetic data provided by Howard, et a l . (5,6) can be easily i n corporated in the design calculations for fluidized beds where the coal residence times are long. However, i f the residence time of pulverized coal in the reactor is short as i t is in entrained bed reactors, then the handling of ordinary differential equations arising from the reaction kinetics require excessive machine computation time. This is due to the stiffness of the differential equations. It is found that the model equations cannot be solved


222



by non-stiff numerical techniques and that the Gear's method re quires excessively long computation times due to the large number of equations in the model. However, using the method described above, the execution time is significantly reduced while main taining the same accuracy as the Gear's method. The Model. The physical system considered is an entrained bed reactor. Pulverized coal, carried by a gas stream, mixes with hot gas at the reactor entrance. As coal particles are carried their temperatures increase and devolatilization takes place. For practical purposes the particle size distribution was approximated by 10 discrete cuts. Since the devolatilization kinetics are high ly temperature dependent and the temperature transient of a p a r t i cle is affected by i t s size and mass separate account must be taken of each of the 15 reactions in each of the 10 different size particles. Without any detail on the derivation of the model, the equations can be summarized as: Mixed gas flowrate F = F + G (1 - e ο o

"

t A

G u

)+m J

v

Σ Σρ.ν.. . M il

(11)

ο .

Mass flowrate of particles with size j m. = m p.(1 - Σ v..) 3 3 i

j = 1, . . . , Κ

(12)

0

Temperature of mixed gas -t/τ T

= [G Cl - e

F

V T

+ (F C o

G

+

m c)T

Ο

ο

- C Σ m.T. - (m Σ Σ p . ν . . ) A h ] / F C j 33 j J !J R

(13)

n

0

Particle temperature dT. dt

1

"

{ T

F - j T

+

9

pr/

* Σ k . (v. οι ι

"V

"HT

- v..)e

-E./RT. m. } / — 1

cpr. ==Jm p. 3k

3

1

iy

±

j = 1,

X

Q

K

(14)

Devolatilization kinetics dv.. * -E./RT. = k.(v. - v..)e dt οι ι ij i = 1, N 1

v

j = 1,

3

J

K


(15)


11.

223

Iterative Approach

KAYiHAN

As shown before the total number of differential equations is K(N + 1 ) . In this study, Κ = 10 and Ν = 15. The choice of 10 is considered to be a reasonable number for the characterization of particle size fractions for design calculations. Therefore, the number of differential equations, 160, was not a r t i f i c i a l l y re duced by taking only a few discrete cuts. However, the p o s s i b i l i t y of representing the particle behavior with one average size was ex plored. Different averages like mean surface, surface mean, v o l ume mean, etc., were t r i e d ; but, none proved to be applicable for this problem where heat transfer and devolatilization occur simul taneously. The stiffness of the system arises partly from the wide range of reaction rates and from the strong temperature dependence of the kinetics given by Eq. (15). The data used in this paper are those obtained by Suuberg (7). As a relative measure of the stiffness of the sysFem the variable time constants of Eq. (15) were calculated at various temperatures between 500°K and 1000°K. The values for -E./RT (k.e 15 1, ) 1

V

01

are given in Table I. There is significant variation between the smallest and largest value evaluated at each temperature where the ratios of the two (the stiffness ratios) range from 10? to 10l6. But with a wide range of particle sizes i t is possible to have a large particle at 500°K while a small particle being already at 1000°K. Therefore, the actual ratio of largest to smallest time constant can be on the order to 2.319 χ 10 /.9831 χ ΙΟ" - Ι Ο which indicates that the eigenvalues of the Jacobian for the sys tem can vary as much as by a factor of 10 . 2ΰ

4

2 4

24

Table I. RXN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

The Wide Range of Time Constants of the Pyrolysis Kinetics (Sec):

800 Κ 600 Κ 700 Κ .8665E+14 .2222E+10 .8005E+06 .5038E+05 .1109E+03 .1126E+01 .6440E+04 .3261E+02 .6190E+00 .1299E+10 .1090E+07 .5374E+04 .2321E+12 .2220E+09 .1207E+07 .9359E+02 .1258E+01 .4965E-01 .3639E+10 .1724E+07 .5536E+04 .2886E+09 .1945E+07 .4573E+05 .2915E+05 .6264E+02 .6252E+00 .2797E+11 .7221E+07 .1471E+05 .6623E+07 .8989E+03 .1130E+01 .1022E+10 .7704E+06 .3504E+04 .1380E+10 .3279E+06 .6274E+03 .4524E+02 .5271E+00 .1869E-01 .8953E+10 .1145E+07 .1377E+04 Stiffness Ratios: .1006E+17 .1915E+13 .4215E+10 .6458E+08

500 Κ .2319E+21 .2647E+09 .1054E+08 .2633E+14 .3916E+16 .3904E+05 .1641E+15 .3165E+12 .1583E+09 .2951E+16 .1719E+13 .2407E+14 .1637E+15 .2305E+05 .2526E+16

900 Κ 1680E+04 .3173E-01 .2835E-01 .8627E+02 .2091E+05 .4019E-02 .6367E+02 .2475E+04 .1737E-01 .1187E+03 .6271E-02 .5284E+02 .4823E+01 .1392E-02 .7379E+01

1000 Κ .1210Ε+02 .1825Ε-02 .2406Ε-02 .3164Ε+01 .8155Ε+03 .5379Ε-03 .1789Ε+01 .2400Ε+03 .9881Ε-03 .2512Ε+01 .9831Ε-04 .1844Ε+01 .9818Ε-01 .1743Ε-03 .1125Ε+00

.1502E+08 .8294Ε+07


224


Another reason for the stiffness of the system is the wide range of thermal transients as indicated before. The time con stant for the thermal transient depends on the square of the p a r t i cle size. Therefore for a system where the ratio of largest to smallest size is 100 the stiffness ratio becomes 10 . For computations the physical properties are assumed to be constant. Particle sizes are discretized by 10 equal-weight cuts from the Rosin-Rammler distribution Downloaded by UNIV OF CALIFORNIA SANTA BARBARA on April 3, 2018 | https://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch011

4

w(x) = exp (-0.0004X )

(16)

1,5

The physical properties and the reactor operating conditions are given in Table II. Table II.

Physical properties and parameters used for calculations

%

(gm/sec)

:

1.0

F

(gm/sec)

:

0.5

(gm/sec)

:

8.5

G

Q

0

T

F o

(°K)

:

400

T

G o

(°K)

:

1200

(

:

0.05

:

0.5

t

G

c

s e c

)

(cal/gm°K)

C

(cal/gm°K)

:

0.7

k

(cal/gm sec °K)

:

1.2 χ 10

(cal/gm)

:

0

(gm/cm )

:

Ah p

R

3

1.4

Solution by a S t i f f Solver. T r i a l solutions with various single size particles representative of the average size deter mined from Eq. (16) indicated that for 80 to 90% completion of de v o l a t i l i z a t i o n and under typical operating conditions the coal residence time in the reactor must be about one second. To check different computation schemes while considering a l l particle sizes only a fraction of this estimated residence time was used. A commercial s t i f f ordinary differential equation solver sub routine, DVOGER, is available in the IMSL Library (3). This sub routine uses Gear's method for the solution of s t i f f ODE s with analytic or numerical Jacobians. The pyrolysis model was solved using DVOGER and the analytical Jacobians of Eqs. (14) and (15). For a residence time of 0.0511 in dimensionless time, defined as t/t where 2 1

-

C p r

AV



11.

Iterative Approach

KAYiHAN

225

and with relative error criterion set at 0.001 the computations took 8.59 minutes on a CDC Cyber 73 computer. With physical prop erties used in this work and with the particle size distribution given by Eq. (16) the characteristic time constant t is 0.1839 seconds. This means that the residence time used was about 1/100 of the expected coal residence time in an actual reactor. Con sidering the stiffness of the kinetic relations over a wide temper ature range and the different temperature time histories of differ ent particle sizes i t is reasonable to expect that the computation times for t/t - 5 or for t - 1.0 sec w i l l be proportionately longer. The main reason for slow computations is the dimensions of the sparse 160 χ 160 Jacobian. Although other routines are available through the work of Hindmarsh (1,2) to handle banded Jacobians or Jacobians with certain structures there are no commercial routines that could be used for the particular Jacobian that arises in the modeling of the low pressure pyrolysis of polydispersed coal p a r t i cles. The Finite Difference Approximations and the Iterative Tech nique. The reaction kinetics are dependent on the rest of the equa tions through particle temperatures. If we have an i n i t i a l guess on the coal particle temperature history during i t s total residence in the reactor then the relations given by Eq. (15) are decoupled from the rest and can be solved separately. Consequently, the v o l a t i l i z a t i o n estimates, v j ( t ) s , can be used to update particle temperatures with Eqs. (11) through (14). F i r s t , consider the approximate analytical integration of Eq. (15). At any time t v o l a t i l i z a t i o n in particle size j due to reac tion i is given by * t -E./RT. f

V

ii

(

t

)

=

V

ii

[ 1

"

" oi

e x p (

k

f

Θ

1

3

d

t

)

]

(

1

8

)

ο Now consider the time increments of At and the discrete points o = 0, t j = At, t2 = 2At, etc. Then, Eq. (18) can be written as * η t -E./RT. v.,(t) = ν [1 - e x p ( - k . Σ / e dt)] (19) ij i j οr m=l t mFor convenience define t -E./RT. J

J

t

m

Φ..(t

) = /

m

e

1

J

1

J

dt

(20)

m-1 = T . ( t ), then the integral is simply " -E./RT.(t ) Φ.. (t ) = At e (21) IJ m However, for any other case where the particle is not iso thermal during a time increment an analytical solution is not available. But, i f i t is assumed that At is small enough so that the temperature profiles within each time increment are straight lines, then the following approximate relation for Φ-jj (t ) can be It T . ( t 3

m

J 1

3

m

1

3

m

m


226


obtained RT.(t ) ι j

J

m

J

0

- ΤA t where Downloaded by UNIV OF CALIFORNIA SANTA BARBARA on April 3, 2018 | https://pubs.acs.org Publication Date: May 30, 1980 | doi: 10.1021/bk-1980-0124.ch011

J(l

s

ι RT.(t J -E./RT.(t ) - 2 J , )e " ] ι 1

J

-E./RT.(t )

1

m

(22)

Τ (t ) - Τ (t ) . ) = -2-3 LJUzL ;j m At (

t

v

T φT r .. (tt )) Φ τ .(t J . j nr j m-1' Now Eq. (19) becomes and

(

j..(t

n

) * v. n' ι

IJ

η Σ Φ . . ( t ))] ij w

[1 - exp ( - k . oi

L

r

v

(23)

J

Assume that it a l l T j ( t ) s and v j ( t ) 's are available from the previous iteration (14) can be :ion or from the i n i t i a l guess. Eq. (14) integrated between t i and t i f the following coefficients are assumed to be constant within the time interval as defined below: [m.(t ) + m.(t )]/2 c p r . a.(t ) = " i-S J (24) 3 n m p. 3k f

n

n

n

n

2

J

pr.

W

•

n

1

0

2

- à -

"

T (t n

Eq. (14)

-

T.(t >

C

J

) + T.(t ) 2

N

J

"V

N

(

2

6

)

gives -At/a.(t )

y

( 2 5 )

A

) + Τ (t ) 2

F

T (t ) = p

C

J

ν

" y v i j

η

)

e

3

ι

"ι V V l Σ

n

+

-At/a.(t ) w

3

J

1

1

+

2ΓΤΓ7 3 η

)β

e

J

n

)

π -At/a (t ) 1

( 2 ? )

With these finite difference approximations the solution can not be achieved in one step but the complete temperature and de v o l a t i l i z a t i o n profiles must be calculated repetitively. As indi cated in Eqs. (22), (23) and (27) the right hand sides have varia bles which must be evaluated at the same time mode as the left hand side. Instead of a simultaneous solution an iterative approach was taken where the right hand sides can be evaluated by the values obtained in the most recent iteration. The complete calculation procedure can be summarized as follows:


11.

Iterative Approach

KAYIHAN

Step 1.

227

Divide total time into NP equal increments, guess a temperature profile for each particle size, and set v.. =0.


A good i n i t i a l guess is the temperature profile that would occur i f there were no devolatilization reactions, but only heat transfer effects on p a r t i cles. Step 2.

Evaluate Eqs. (11), crements.

(12) and (13)

for a l l time i n -

Step 3.

Evaluate Eq. (23) for a l l particle sizes, reactions and time increments.

Step 4.

Update temperature profiles through Eq. (27).

Step 5.

Stop i f convergence is achieved; otherwise, go to Step 2.

This procedure had converged in 4 or 5 iterations to four significant figures for a l l cases tried in this study. The accuracy of the calculations depends on the time increment At because the finite difference approximations become more accurate as At gets smaller. A summary of some iteration results and a comparison between this technique and the numerical integration with Gear's method w i l l be presented after the following discussion on the s t a b i l i t y of the temperature equation. Stability of the Temperature Equation. Consider Eq. (27) where T j ( t ) is dependent on T j ( t i ) and T j ( t ) computed by the previous iteration through the f i r s t term on the right and through bj(t ). The latter can cause i n s t a b i l i t y in the explicit iterative calculations for some choices of At. F i r s t , consider an idealized case where c = C and Ah = 0. Then b j ( t ) = 0 and Eq. (27) takes the form -At/a.(t ) T.(t ) - T.(t J e -At/a.(t ) T (t )(l - e ) (28) n

n

n

n

R

3

3

n

3

+

n

p

n

n

X

3

n

n

For this ideal case the iterations for T j ( t ) w i l l be stable for a l l choices of At. Therefore, this approximation is a good start for exploratory reactor design calculations because of its simp l i c i t y and s t a b i l i t y . Usually, however, the above assumptions w i l l not hold. To establish a safe At for stable calculations observe that b j ( t ) in Eq. (27) is on the order of pr. h [c - C)T.(t - ihl 3k j n-r R n

n

2

v


228


Substituting this into Eq. (27)

T.(t ) * T . ( t . ) j

η

3

n-1 2


3 - Σ v..

{e-

η

(t

p

j V

A t / a

J(l

(

J

ι

* T (t )(l - e

gives

j η -At/a (t ) 1>

J

- - Γ Τ

+

-At/a

n

J

(t )

pr ) - -gjL- Α Η [ Σ ν . . ( κ

-At/a -Σν..(ΐ

ij' n-l

i

JK

2a

L—)

j(V

e

ν (

1-

(t ) 3

n

(29)

A safe choice for At w i l l be At * 2a.(t

3

η

)

(30)

This would be restricted by the smallest particle size, j = 1, and by the extent of devolatilization. Usually, the volatile matter of coal is not more than 50% by weight. Therefore, choose cpr 2

With physical properties given in Table I, this becomes At = 1944 Γ and for r

χ

2

= 10 pm

At = 0.0019 sec or NP = 514 peints for 1 second integration time. While, for r

= 5 ym -4

At = 4.9 χ 10

sec

and NP = 2058 points. Because i t is necessary to adjust the time increment according to the smallest particle size the number of discrete points, NP, may be unnecessarily high. This would increase the computation time without significantly improving the accuracy. In this study, instead of adopting the above s t a b i l i t y criterion, i t was ob served that the smaller particles w i l l respond to temperature


Iterative Approach

11. KAYiHAN

229

changes in their environment much faster than the large particles. Therefore, i t was assumed that for an arbitrary At i f a j ( t ) , the i n i t i a l temperature time constant of the particle, is smaller than At, then i t w i l l respond to changes in Tp instantaneously, i . e . , T . ( t ) = T ( t ) i f a j ( t ) * At. This allows At to be adjusted un t i l desired accuracy is achieved in the converged iterative results. The procedure described here has provided a stable and an accurate iterative technique without unnecessarily large number of time i n crements. 0


n

F

n

0

Comparison of the Computation Results. As indicated above, Gear's method was used to solve the model equations only for a fraction of the total residence time in the reactor which took 8.59 minutes of machine computation time. The same set of equa tions was solved by the approximate iterative technique for the same time interval in 5.8 seconds of computer time. As a com parison of the accuracy overall devolatilization v = Σ Σ v-jj as predicted by the two techniques are plotted on a dimensionless scale in Figure 1. The definitions for the dimensionless quanti ties used are: T

Dimensionless time = t/t Dimensionless devolatilization = v / v * where t is the time constant for the temperature transient of nondevolatilizing average size particle and v * is the mass fraction of total volatile matter in coal. Typical reactor calculations would take about five times the average particle time constant t. For these cases the approximate relations were used to predict some of the reactor conditions as a function of time. To give an idea of the dependence of the com putation time and accuracy on the time increment, At, three dif ferent cases are compared in Table III. These results show that the proposed technique provides a fast and reliable method for the solution of s t i f f ODE models of reacting polydispersed particles. Recently, Turton (9) applied this method successfully to the modeling of wood char combustion in a transport reactor. T

CONCLUSIONS Models for the reacting polydispersed particles contain s t i f f ordinary differential equations. The stiffness is due partly to the wide range of thermal time constants of the particles and part ly to the high temperature dependence of reactions like combustion and devolatilization. As an alternative to the established solu tion techniques based on Gear's method an iterative approach is developed which uses the f i n i t e difference representations of the differential equations. The finite differences are obtained by



Figure 1.

0.00

— Ω? -

.08

".01

Δ

Δ

DIFFERENCE

'.02 ".03 DIMENSIONLESS

Δ

FINITE

COAL

TIME

'.04

SOLUTION

POLYDISPERED

(IMSL)

OF

Δ

D

'.05

•Δ

PARTICLES

'.06

Comparison of the iterative finite difference solution of this work to Gear's method

Δ

ITERATIVE

Δ

METHOD

GEAR'S

•

DEVOLATILIZATION


11. KAYiHAN

Iterative Approach

231

Table III. Sample calculations for particle temperatures and for total devolatilization showing the effect of At on accuracy and execution time. Number of Points =100


Number of Points = 50 Time 0.0000 0. 5000 1. 0000 1. 5000 2. 0000 2. 5000 3. 0000 3. 5000 4. 0000 4. 5000 5. 0000

T4 .3333 .8113 .9076 .9105 .9077 .9056 .9041 .9032 .9024 .9018 .9014

T7 .3333 .5563 .7063 .8044 .8579 .8834 .8947 .8994 .9011 .9016 .9015

VT .0000 .3173 .4808 .5832 .6600 .7096 .7536 .7794 .8060 .8369 .8564

11.879 CP seconds execution time

Time 0. 0000 0. 5000 1. 0000 1. 5000 2. 0000 2. 5000 3. 0000 3. 5000 4. 0000 4. 5000 5. 0000

T4 .3333 .8252 .9085 .9104 .9075 .9055 .9041 .9031 .9024 .9018 .9014

T7 .3333 .5660 .7123 .8083 .8598 .8843 .8951 .8996 .9012 .9016 .9015


Number of Points = 150 Time 0.0000 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000 4.0000 4.5000 5.0000

T4 .3333 .8292 .9087 .9103 .9075 .9055 .9041 .9031 .9024 .9018 .9013

T7 .3333 .5692 .7143 .8095 .8604 .8845 .8952 .8996 .9012 .9016 .9014

VT .0000 .3156 .4814 .5838 .6604 .7116 .7541 .7800 .8072 .8382 .8566

VT .0000 .3173 .4816 .5838 .6605 .7112 .7541 .7799 .8069 .8380 .8566

Time = t/t T4 T7 VT

V Go T

=

V Go T

VV




232


approximate analytical integrations through small time increments. This technique allows successive iterations on the complete solu tion through sequential evaluations of the equations rather than a simultaneous approach. Application to the modeling of polydis persed coal pyrolysis at low pressures gives encouraging results. For this case the direct solution of the model equations with a s t i f f solver like DVOGER in the IMSL Library is expected to take more than 10 hours on a CDC Cyber 73. Other s t i f f ODE solvers that can handle particular Jacobian structures are not applicable to this problem. The method proposed here reduces the computation time to less than 30 seconds with three digit accuracy in the com puted variables. This iterative technique would be practical for reactor design and optimization studies. NOMENCLATURE 2

a. S3 c

surface area of particle, cm

C

specific heat of gas, cal/gm°K

E

activation energy of i t h reaction, cal/gmol

F

mass flowrate of mixed gas, gm/sec

v

specific heat of coal, cal/gm°K

F

Q

i n i t i a l mass flowrate of coal carrier gas, gm/sec

G

Q

i n i t i a l mass flowrate of hot gas, gm/sec

k k

thermal conductivity of gas, cal/gm sec °K Q

Κ

pre-exponential factor of ith reaction,

1/sec

number of discrete particle size cuts

m

coal mass flow rate at reactor entrance, gm/sec

ΐϊκ

mass flowrate of char particles in the jth size cut, gm/sec

Ν

number of gas-solid reactions

p. J r r

mass fraction of m that is in size cut j ο particle radius of size cut j , cm root-mean-square average radius of particle size distribution, cm

Q

F

v

R

gas constant, cal/gmol °K

t t

time, sec time constant for the temperature transient of size r , fined in Eq. (17), sec

Tp

temperature of gas stream F , °K

Tp ο

temperature of gas stream F , °K


de

11.

KAYIHAN

T

Iterative Approach

temperature of gas stream G , °K

b


233

O

T\ v.*

temperature of particles of size j , °K maximum v o l a t i l i z a t i o n due to reaction i , gm of volatiles/gm of original coal

vj

volatiles produced in size j through reaction i , gm/gm

w

Rosin-Rammler mass fraction size distribution function

χ

particle diameter in w(x), cm

Ah

R

heat of reaction, cal/gm

At

time increment for discrete solution, sec

p

density of coal, gm/cm

Tç

characteristic mixing time of hot and cold gas streams, sec

3

ACKNOWLEDGEMENT The computer time was provided by the Oregon State University Milne Computer Center. LITERATURE CITED 1.

Hindmarsh, A.C., "GEARB: Solution of Ordinary Differential Equations Having Banded Jacobians," Lawrence Livermore Laboratory Report UCID-30059, June 1977, Rev. 2.

2.

Hindmarsh, A.C., "Preliminary Documentation of GEARBI: Solution of ODE Systems with Block Iterative Treatment of the Jacobian," Lawrence Livermore Laboratory Report UCID-30149, December 1976.

3.

The IMSL Library, International Mathematical and Libraries, Inc., July 1977, 6th Ed.

4.

Levenspiel, O., "Chemical Reaction Engineering," 1972, Second Edition, Wiley.

5.

Anthony, D.B., Howard, J.B.,

6.

Suuberg, E.M., Peters, W.A., Howard, J.B., and Development, 1978, 17, 37.

7.

Suuberg, E.M., Sc.D. Thesis, Dept. of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA, 1977.

8.

Kayihan, F., Reklaitis, G.V., "Modelling of Staged Fluidized Bed Coal Pyrolysis Reactors," Paper presented at 85th National AIChE Meeting, June 1978.

9.

Turton, R., "Combustion of Wood Char in a Transport Reactor," M.S. Thesis, Dept. of Chemical Engineering, Oregon State University, Corvallis, OR, 1979.

AIChE J.,

Statistical

1976, 22, 625. I&EC Process Design

Squires and Reklaitis; November 5, 1979.Computer Applications to Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1980.

RECEIVED

Computer Applications to Chemical Engineering - American Chemical

Recommend Documents