Moving Bed Coal Gasifier Dynamics Using MOC and MOL Techniques

17 Moving Bed Coal Gasifier Dynamics Using MOC and M O L Techniques RICHARD STILLMAN

Downloaded by CORNELL UNIV on October 15, 2016 | http://pubs.acs.org Publication Date: September 21, 1981 | doi: 10.1021/bk-1981-0168.ch017

IBM Scientific Center, P.O. Box 10500, Palo Alto, C A 94304

The method of characteristics and both method of lines techniques, continuous-time discrete-space and continuous-space discrete-time, were used to solve the system of hyperbolic partial differential equations representing the dynamics of a moving bed coal gasifier with countercurrent gas-solid heat transfer. The adiabatic plug flow model considers 17 solids stream components, 10 gas stream components and 17 reactions. The kinetic and thermodynamic parameters were derived for a Wyoming subbituminous coal. The inherent numerical stiffness of the coupled gas-solids equations was handled by assuming that the gas stream achieved steady state values almost instantly. Calculated dynamic responses are shown for step changes in reactor pressure, blast temperature, steam flow rate, and coal moisture. Both steady state convergent and limit cycle responses were obtained.

The chemical i n d u s t r y is beginning to s h i f t away from dependence on n a t u r a l gas and petroleum to the use of c o a l as the b a s i s f o r some of t h e i r hydrocarbon feed s t o c k s . Examples can be found in ammonia manufacture, methanol p r o d u c t i o n , a c e t i c anhydride s y n t h e s i s and s y n t h e t i c g a s o l i n e p r o d u c t i o n . Electric utilities are a l s o l o o k i n g at c o a l f o r use in combined c y c l e and f u e l c e l l based power p l a n t s . A common f i r s t step in many of these p l a n t s is the g a s i f i c a t i o n of c o a l , and s i n c e g a s i f i e r o p e r a t i o n can impact the operation of other u n i t s in the p l a n t , it is u s e f u l to p r e d i c t the dynamic behavior of a c o a l g a s i f i e r r e a c t o r using a s i m u l a t i o n model. The mathematical model can help to provide a b e t t e r understanding of the complex dynamic behavior e x h i b i t e d by the a c t u a l g a s i f i e r when subjected to simple and m u l t i p l e d i s t u r bances. The extent of dynamic t e s t i n g that can be performed on a s i m u l a t i o n model, as w e l l as the type of information which can be 0097-6156/81/0168-0331$09.00/0 © 1981 American Chemical Society Fogler; Chemical Reactors ACS Symposium Series; American Chemical Society: Washington, DC, 1981.


332

CHEMICAL REACTORS

obtained, o f t e n goes much beyond that which would be allowable or even p o s s i b l e in a given commercial i n s t a l l a t i o n . Some work has already been done on the s i m u l a t i o n of t r a n s i e n t behavior of moving bed c o a l g a s i f i e r s . However, the analys i s is not based on the use of a t r u l y dynamic model but i n s t e a d uses a steady s t a t e g a s i f i e r model p l u s a pseudo steady state approximation. For this type of approach, the time response of the g a s i f i e r to r e a c t o r input changes appears as a continuous sequence of new steady s t a t e s . Yoon, Wei and Denn (1, 2^, 3) consider the time response of a g a s i f i e r to small changes in operating c o n d i t i o n s such as might occur during normal operation of the r e a c t o r . They regard the time r e q u i r e d to reach a new steady s t a t e , f o l l o w i n g a step change in operating c o n d i t i o n s , as the most u s e f u l measure of t r a n s i e n t response. For small step changes, they estimate this response time, and the changes in r e a c t o r v a r i a b l e s during a transient, using a psuedo steady s t a t e technique. T h e i r technique i n v o l v e s removal of the time v a r i a b l e from the system of dynamic equations by assuming that the space o r i g i n moves at the same v e l o c i t y as the v e l o c i t y of the thermal wave f o r the maximum bed temperature. They give d e t a i l e d c a l c u l a t e d r e s u l t s f o r small step changes in c o a l feed r a t e s f o r both ash discharge and s l a g g i n g g a s i f i e r s . The g a s i f i e r modeling technique used by Hsieh, Ahner and Quentin (4) is based on the c o n s t r u c t i o n of a data space representing steady s t a t e r e a c t o r c o n d i t i o n s using the U n i v e r s i t y of Delaware steady s t a t e model described in Yoon, Wei and Denn (1). They made an a n a l y s i s of the chemical r e a c t i o n r a t e s and thermal capacitance e f f e c t s to develop the dynamic algorithms used to simulate the dynamic trends between the steady s t a t e points. T h e i r dynamic responses are thus estimated by using the quasi steady s t a t e data bank, a l i n e a r i n t e r p o l a t i o n r o u t i n e and the derived dynamic algorithms. No a c t u a l d e t a i l s of t h e i r model are given. However, they do show a short 9 minute g a s i f i e r t r a n s i e n t response f o r e x i t gas composition and temperature r e s u l t i n g from a ramp decrease in steam, a i r and c o a l feed r a t e s . D a n i e l (5) has used a s i m p l i f i e d method to develop a short time s c a l e t r a n s i e n t model f o r a moving bed g a s i f i e r . Wei (6) presents a very b r i e f d i s c u s s i o n of c o a l g a s i f i c a t i o n r e a c t o r dynamics. He describes the t r a n s i e n t response to a small step change as a s o f t t r a n s i e n c e in which the movement from one steady s t a t e to another one nearby takes place as a wave through a s e r i e s of pseudo steady s t a t e s . He p o i n t s out that the hard t r a n s i e n c e of s t a r t up and major upset in r e a c t o r operation are not w e l l understood. One of the purposes of this paper is to increase this understanding. Although the pseudo steady s t a t e approximation provides a u s e f u l t o o l f o r estimating some aspects of g a s i f i e r dynamics, it does not provide the means to examine the f u l l range of dynamic behavior that one would expect to f i n d f o r a g a s i f i e r . Therefore, a d i f f e r e n t approach has been taken here in that a nonlinear

Fogler; Chemical Reactors ACS Symposium Series; American Chemical Society: Washington, DC, 1981.

17.

STiLLMAN

Moving Bed

Coal Gasifier Dynamics

333

dynamic model r e p r e s e n t i n g moving bed g a s i f i e r dynamics is solved d i r e c t l y so that the g l o b a l aspects of the dynamic behavior can be examined. One f u r t h e r note, the U n i v e r s i t y of Delaware g a s i f i e r model used in the pseudo steady s t a t e approximation assumes that the gas and s o l i d s temperatures are the same w i t h i n the r e a c t o r . That assumption removes an important dynamic feedback e f f e c t between the countercurrent flowing gas and s o l i d s streams. This is p a r t i c u l a r l y important when the burning zone moves up and down w i t h i n the r e a c t o r in an o s c i l l a t o r y manner in response to a step change in operating c o n d i t i o n s .


Process D e s c r i p t i o n A moving bed g a s i f i e r is a v e r t i c a l r e a c t o r with countercurrent flowing gas and s o l i d s streams. Coal enters the top of the r e a c t o r and ash (and/or c l i n k e r or molten slag) is removed from the bottom. A mixture of steam and oxygen (or a i r ) enters the bottom of the r e a c t o r and the raw product gas e x i t s from the top. An a d i a b a t i c steady s t a t e plug flow model has been developed by S t i l l m a n (7, 8) f o r this type of g a s i f i e r . For that model, the f o l l o w i n g sequence of p h y s i c a l and chemical events was assumed to take place in the r e a c t o r . Heat was extracted from the hot e x i t i n g gas by the e n t e r i n g s o l i d s stream so that the c o a l temperature was i n c r e a s e d and the c o a l moisture was evaporated. A f u r t h e r increase in the c o a l temperature caused the c o a l v o l a t i l e matter to be r e l e a s e d , l e a v i n g a char. In the g a s i f i c a t i o n zone, some of the char reacted with the carbon d i o x i d e , water and hydrogen gas stream components. The oxygen in the feed gas burned all or almost all of the remaining char to provide the heat necessary to run this endothermic process. An ash r e s i d u e was l e f t from the combustion r e a c t i o n . Some or all of this ash melted and then e i t h e r s o l i d i f i e d to form c l i n k e r s or e l s e remained in a molten s t a t e , depending on the r e a c t o r o p e r a t i n g c o n d i t i o n s . The water gas s h i f t r e a c t i o n and the methanation r e a c t i o n were a l s o assumed to take place in the gas stream. The model considered 17 components in the s o l i d s stream: water, hydrogen, n i t r o g e n , oxygen, carbon, s u l f u r , ash, s l a g , c l i n k e r , w a t e r ( v s ) , hydrogen(vs), carbon d i o x i d e ( v s ) , carbon monoxide(vs), methane(vs), hydrogen s u l f i d e ( v s ) , ammonia(vs), tar(vs), where (vs) i n d i c a t e s a v o l a t i l e s o l i d component. The gas stream had 10 components: water, hydrogen, n i t r o g e n , oxygen, carbon d i o x i d e , carbon monoxide, methane, hydrogen s u l f i d e , ammonia, t a r . A set of 17 r e a c t i o n s was w r i t t e n to simulate the r e a c t o r events and they included 1 r e a c t i o n f o r d r y i n g , 8 p a r a l l e l r e a c t i o n s f o r d e v o l a t i l i z a t i o n , 5 r e a c t i o n s f o r g a s i f i c a t i o n and 3 r e a c t i o n s f o r combustion. The k i n e t i c and thermodynamic parameters f o r these r e a c t i o n s were derived f o r a Wyoming subbituminous c o a l .


CHEMICAL REACTORS

334

The same model d e s c r i p t i o n w i l l be used as the b a s i s d e r i v i n g the g a s i f i e r dynamic model. A l l of the k i n e t i c thermodynamic parameters w i l l be taken from S t i l l m a n (7^, 8 ) .

for and

Dynamic Model


The c o n t i n u i t y equations f o r mass and energy w i l l be used to d e r i v e the h y p e r b o l i c p a r t i a l d i f f e r e n t i a l equation model f o r the s i m u l a t i o n of moving bed c o a l g a s i f i e r dynamics. Plug flow (no a x i a l dispersion) and a d i a b a t i c (no radial gradients) o p e r a t i o n w i l l be assumed. The mass balance dynamic equations f o r the s o l i d s stream are given by ^£i +il£l= 3

t

3

2

ri 1

1=1,2....,17 3=1,2

17

where CjS is the c o n c e n t r a t i o n of s o l i d s component j and the molar f l u x of s o l i d s component j defined as FjS

Ξ CjS uS

FjS

is

j=l,2,...,17.

The a i j values are the s t o i c h i o m e t r i c c o e f f i c i e n t s f o r component j in r e a c t i o n i , r i is the r a t e of reaction i , uS is the l o c a l s o l i d s stream v e l o c i t y , t is r e a l time, and ζ is distance measured from the bottom of the r e a c t o r . The corresponding mass balance dynamic equations f o r the countercurrent flowing gas stream are given by fCjG _ 3

t

= £ 3

z

1

a

i

j

r

l

1

β

1

§

2

§

. . . ,

7

j=18,19,...,27

where the molar f l u x of gas component j FjG = CjG uG

1

is

defined as

j=18,19,...,27.

The energy balance dynamic equation f o r the s o l i d s stream 3 < 1 > S

+ ^

8t

9z

is

= hGS AGS(TG - TS) - £ r i ΔΗ1 i i-1,2,...,12,15,16,17

where the energy d e n s i t y and energy are defined by

f l u x f u n c t i o n s f o r the

| S = pS CS TS C


solids

17.

Moving Bed Coal Gasifier Dynamics

STiLLMAN

335

and


= G(t,z) TS(t,z) TG(t,z).

j=l,2,...,17 j=18,19,...,27

t=0

z-0,...,L

The boundary c o n d i t i o n s needed f o r the model are the input molar fluxes and the temperatures of the s o l i d s and gas feed streams, and the i n l e t gas p r e s s u r e . At the top of the r e a c t o r :


336

CHEMICAL

FjS(t.z) TS(t,z).

j-l.2

9 * * * 9

17

t=t

z=L

t=t

z=0

REACTORS

At the bottom of the r e a c t o r :


FjG(t,z) TG(t,z) P(t,z).

j-18,19 9

·

·

·

9

27

T h i s set of h y p e r b o l i c p a r t i a l d i f f e r e n t i a l equations f o r the g a s i f i e r dynamic model represents an open or s p l i t boundary-value problem. S t a r t i n g with the initial c o n d i t i o n s w i t h i n the r e a c t o r , we can use some type of marching procedure to s o l v e the equations d i r e c t l y and to move the s o l u t i o n forward in time based on the s p e c i f i e d boundary c o n d i t i o n s f o r the i n l e t gas and i n l e t s o l i d s streams. However, it is important to note that there is an inherent numerical s t i f f n e s s in the coupled g a s - s o l i d s equations because the gas stream moves through the r e a c t o r much more r a p i d l y than the s o l i d s stream. In a t y p i c a l example, while it only takes the gas about 7 seconds to move through the r e a c t o r , it takes the s o l i d s stream about a 1000 times longer. T y p i c a l r a t i o s of gas v e l o c i t y to s o l i d s v e l o c i t y are about 400, 4200, 1200, at the top of the r e a c t o r , in the burning zone, and at the bottom of the r e a c t o r , r e s p e c t i v e l y . The s o l i d s and gas v e l o c i t i e s represent the two c h a r a c t e r i s t i c d i r e c t i o n s f o r our h y p e r b o l i c system. I f we p l o t these v e l o c i t y curves on a r e a c t o r length versus time graph, the c h a r a c t e r i s t i c curves f o r the gas w i l l be e s s e n t i a l l y h o r i z o n t a l in comparison to the s o l i d s stream c h a r a c t e r i s t i c because of the l a r g e gas to s o l i d s v e l o c i t y r a t i o s . Making the assumption that the gas stream characteristic curves are indeed h o r i z o n t a l f o r all p r a c t i c a l purposes, is equivalent to s e t t i n g the time p a r t i a l d e r i v a t i v e s f o r the concent r a t i o n s and the energy d e n s i t y equal to zero in the original system of p a r t i a l d i f f e r e n t i a l equations f o r the gas. Using this approximation reduces the gas equations to a set of steady s t a t e equations. Thus our f i n a l dynamic model f o r a moving bed c o a l g a s i f i e r c o n s i s t s of a set of h y p e r b o l i c p a r t i a l d i f f e r e n t i a l equations f o r the s o l i d s stream coupled to a set of o r d i n a r y d i f f e r e n t i a l equations f o r the gas stream. Shampine and Gear (9) c a u t i o n t h a t , f o r systems c o n t a i n i n g elements with d i f f e r e n t time scales, removing s t i f f n e s s by changing the model may be r i s k y because it might be d i f f i c u l t to r e l a t e the s o l u t i o n of the modified model to that of the original model. In our case, no such d i f f i c u l t y was found. The steady s t a t e c o n d i t i o n s p r e d i c t e d by the modified dynamic model were monitored by the steady s t a t e model 8)·


17.


STILLMAN

337


Distance Method of L i n e s The distance method of l i n e s (continuous-time discrete-space) technique is a s t r a i g h t f o r w a r d way f o r o b t a i n i n g the numerical s o l u t i o n of time dependent p a r t i a l d i f f e r e n t i a l equations with one spatial variable. The original system of p a r t i a l d i f f e r e n t i a l equations is transformed i n t o a coupled system of time dependent ordinary d i f f e r e n t i a l equations by using s p a t i a l f i n i t e d i f f e r e n c e formulas to r e p l a c e the s p a t i a l differentiation terms for a discrete set of spatial grid points. The number of o r d i n a r y d i f f e r e n t i a l equations produced by this o p e r a t i o n is equal to the original number of p a r t i a l d i f f e r e n t i a l equations m u l t i p l i e d by the number of g r i d p o i n t s used. Thus, although we now have a l a r g e r number of equations to c o n s i d e r , u s u a l l y the augmented system of o r d i n a r y d i f f e r e n t i a l equations is e a s i e r to solve n u m e r i c a l l y than the original smaller system of p a r t i a l d i f f e r e n t i a l equations. First order hyperbolic d i f f e r e n t i a l equations transmit d i s c o n t i n u i t i e s without d i s p e r s i o n or d i s s i p a t i o n . Unfortunately, as Carver (10) and Carver and Hinds (11) p o i n t out, the use of s p a t i a l f i n i t e d i f f e r e n c e formulas introduces unwanted d i s p e r s i o n and spurious o s c i l l a t i o n problems i n t o the numerical s o l u t i o n of the differential equations. They suggest the use of upwind d i f f e r e n c e formulas as a way to d i m i n i s h the o s c i l l a t i o n problem. T h i s follows d i r e c t l y from the concept of domain of i n f l u e n c e . For h y p e r b o l i c systems, the domain of i n f l u e n c e of a given v a r i a b l e is downstream from the p o i n t of r e f e r e n c e , and t h e r e f o r e , a n a t u r a l consequence is to use upstream d i f f e r e n c e formulas to estimate downstream c o n d i t i o n s . When necessary, the unwanted d i s p e r s i o n problem can be reduced by using low order upwind d i f f e r e n c e formulas. The Lagrange i n t e r p o l a t i o n polynomial was used to develop the s p a t i a l f i n i t e d i f f e r e n c e formulas used f o r the distance method of lines calculation. For example, the two point polynomial f o r the s o l i d s f l u x v a r i a b l e F ( t , z ) can be expressed by 2

z

F(t,z) • / " \z(k-l)

(

k

)

W k - 1 ) - z(k)/

Z

Z

(

k

1

}

+ ( " " W k ) \z(k) - z(k-l)/

where k represents the g r i d point index number. The index number increases in value from top to bottom of the r e a c t o r . I f we take the p a r t i a l d e r i v a t i v e of the two point polynomial with respect to ζ at index point k, we o b t a i n 3F(t,z) 9z

m

F(t,k) -

F(t,k-1)

Δζ

which is the two point upwind formula f o r the s o l i d s stream. The l o c a l g r i d spacing is i n d i c a t e d by Δ ζ . In a s i m i l a r f a s h i o n , more accurate higher order formulas can be developed. The four point upwind formula is


338

CHEMICAL REACTORS 3F(t,z)

m

H F ( t , k ) - 18F(t,k-l) + 9F(t,k-2)

-

2F(t,k-3).

6ΔΖ

3z

The four p o i n t upwind biased formula is given by 3F(t,z)

m

2F(t,k+l)

+ 3F(t,k) - 6 F ( t , k - l ) + F(t,k-2) 6ΔΖ

3z

and the four p o i n t downwind biased formula 8F(t,z)


8z

m

-F(t,k+2)

+ 6F(t,k+l)

is

- 3F(t,k) - 2 F ( t , k - l ) . 6ΔΖ

D i f f e r e n t combinations of s p a t i a l f i n i t e d i f f e r e n c e formulas were t r i e d to determine the best set f o r our system of equations. The two p o i n t upwind formula was found to be best f o r the s o l i d s component molar f l u x e s . The low order formula was used because most of the g a s i f i e r r e a c t i o n s t u r n o f f a b r u p t l y when a component disappears and this creates sharp d i s c o n t i n u i t i e s . Higher order formulas tend to f l a t t e n out d i s c o n t i n u i t i e s , and in some cases, this causes m a t e r i a l balances to be l o s t which then leads to numerical i n s t a b i l i t y problems. M a i n t a i n i n g component m a t e r i a l balance is an important a i d to p r e s e r v i n g numerical s t a b i l i t y in the c a l c u l a t i o n s . The low order formulas minimized these d i f f i culties. The four point upwind biased formula worked best f o r the s o l i d s stream energy f l u x c a l c u l a t i o n . Some downstream informa t i o n was u s e f u l because of the countercurrent flow of the gas and s o l i d s streams. To keep the same o r d e r , the four point downwind biased formula was used at the top of the r e a c t o r and the four p o i n t upwind formula was used at the bottom. Accuracy and c a l c u l a t i o n time are h i g h l y dependent on the number of s p a t i a l g r i d p o i n t s used. More g r i d p o i n t s give b e t t e r accuracy but c a l c u l a t i o n time increases a c c o r d i n g l y . To r e s o l v e this dilemma, a v a r i a b l e g r i d s t r u c t u r e was used in the c a l c u l a tions. In the top part of the r e a c t o r , where the d r y i n g and d e v o l a t i l i z a t i o n reactions were t a k i n g p l a c e , a coarse g r i d was used. In the bottom part of the r e a c t o r , where the g a s i f i c a t i o n and combustion r e a c t i o n s were o c c u r r i n g , a f i n e r mesh was used. A t o t a l of 82 g r i d p o i n t s was used f o r the method of l i n e s c a l c u l a tions. With the v a r i a b l e g r i d s t r u c t u r e , the top t h i r d of the r e a c t o r had 13 nodes, the middle t h i r d had 21 nodes, and the bottom t h i r d had 48 nodes. G r i d spacing has no e f f e c t on the use of the two p o i n t upwind formula but it does e f f e c t the use of the four p o i n t formula. Therefore, g r i d r e d u c t i o n was done in a p r e s c r i b e d manner. For any one change, the g r i d spacing could only be cut in h a l f and the g r i d change had to remain in e f f e c t f o r at l e a s t 3 node p o i n t s .


17.


STiLLMAN

339

T h i s r e s t r i c t i o n allowed the c o e f f i c i e n t s f o r the four point upwind biased formula to sequence through the f o l l o w i n g values f o r one g r i d spacing change


2.0 3.2 2.25 2.0

3.0 -1.5 2.0 3.0

-6.0 -2.0 -4.5 -6.0

1.0 O.3 O.25 1.0

f u l l step s i z e h a l f step s i z e h a l f step s i z e h a l f step s i z e

which is then repeated f o r each succeeding change. The v a r i a b l e s needed in the dynamic model c a l c u l a t i o n s are the s o l i d s molar and energy f l u x e s . However, in the distance method of l i n e s technique, when we r e p l a c e the s p a t i a l d i f f e r e n t i a l terms by f i n i t e d i f f e r e n c e formulas, the time d e r i v a t i v e in the remaining d i f f e r e n t i a l equation is in terms of e i t h e r compo nent c o n c e n t r a t i o n or energy d e n s i t y , which we do not want. T h e r e f o r e , the f o l l o w i n g change was made in the distance method of l i n e s model. Replacing CjS with i t s equivalent FjS/uS and taking the p a r t i a l d e r i v a t i v e with respect to time gives 3

/

F

J

S

\

at" Vus /

1

F

3 JS

us

at

FjS

3uS

us us at

N e g l e c t i n g the a c c e l e r a t i o n term, we have the approximation acjs

at

_ ι

us

aFjs

at

which gives us the d e s i r e d f l u x v a r i a b l e in the model. A s i m i l a r change was used to convert from energy d e n s i t y to energy f l u x . Time Method of L i n e s The time method of l i n e s (continuous-space discrete-time) technique is a h y b r i d computer method f o r s o l v i n g p a r t i a l d i f f e r e n t i a l equations. However, in i t s standard form, the method gives poor r e s u l t s when c a l c u l a t i n g t r a n s i e n t responses f o r h y p e r b o l i c equations. M o d i f i c a t i o n s to the technique, such as the method of decomposition (12), the method of d i r e c t i o n a l d i f f e r e n c e s (13), and the method of c h a r a c t e r i s t i c s (14) have been used to c o r r e c t this problem on a h y b r i d computer. To make a comparison with the d i s t a n c e method of lines and the method of characteristics r e s u l t s , the technique was used by us in i t s standard form on a d i g i t a l computer. The original system of p a r t i a l d i f f e r e n t i a l equations is transformed i n t o a system of o r d i n a r y d i f f e r e n t i a l equations by r e p l a c i n g the time d i f f e r e n t i a l terms with time f i n i t e d i f f e r e n c e formulas. The number of equations in the new system is the same as the original number of equations. However, it is necessary to s t o r e intermediate r e s u l t s at s p a t i a l nodes f o r both current and previous time increments.


CHEMICAL REACTORS

340

The Lagrange i n t e r p o l a t i o n polynomial was again used to develop the f i n i t e d i f f e r e n c e formulas. To avoid a d d i t i o n a l iterations, only upwind d i f f e r e n c e s were used. The two p o i n t upwind formula f o r the s o l i d s stream c o n c e n t r a t i o n v a r i a b l e at any l o c a t i o n ζ w i t h i n the r e a c t o r f o r time t is given by 8C(t,z)

β

at

C(t,z) - C ( t - l , z )

At

where At is the time increment. 3C(t,z)

m

3C(t,z) - 4 C ( t - l , z ) + C(t-2,z)

at Downloaded by CORNELL UNIV on October 15, 2016 | http://pubs.acs.org Publication Date: September 21, 1981 | doi: 10.1021/bk-1981-0168.ch017

The three p o i n t upwind formula is

2At

and the four p o i n t upwind formula is 8C(t,z)

m

l l C ( t , z ) - 18C(t-l,z) + 9C(t-2,z) - 2C(t-3,z).

8t

6At

The same 82 p o i n t v a r i a b l e g r i d s t r u c t u r e was used in the time method of l i n e s c a l c u l a t i o n s as was used f o r the d i s t a n c e method of l i n e s c a l c u l a t i o n s . A l s o , the three and f o u r p o i n t upwind formulas were found to attenuate the c a l c u l a t e d step responses too much and they were d i s c a r d e d . Method of C h a r a c t e r i s t i c s The method of c h a r a c t e r i s t i c s (15) is a n a t u r a l way for s o l v i n g h y p e r b o l i c p a r t i a l d i f f e r e n t i a l equations. The technique is based on l o c a t i n g the c h a r a c t e r i s t i c propagation paths or d i r e c t i o n s f o r the p a r t i a l d i f f e r e n t i a l equations and i n t e g r a t i n g the r e s u l t i n g o r d i n a r y d i f f e r e n t i a l equations along these d i r e c t i o n s . Thus, as with the method of l i n e s , this technique t r a n s forms our problem from s o l v i n g p a r t i a l d i f f e r e n t i a l equations to s o l v i n g o r d i n a r y d i f f e r e n t i a l equations. For the g a s i f i e r dynamic model, the c h a r a c t e r i s t i c d i r e c t i o n s are given by the s o l i d s and gas stream v e l o c i t i e s : dz = uS

+ direction

= -uG

- direction.

dt dz dt These two f a m i l i e s of curves c r o s s each other at a number of common nodes. However, the assumption was made e a r l i e r that only the steady s t a t e equations would be used f o r the gas stream c a l c u l a t i o n s . Therefore, we only need to consider the a p p l i c a t i o n of the method of c h a r a c t e r i s t i c s to the s o l i d s stream partial d i f f e r e n t i a l equations.


17.


STiLLMAN

341


To use the technique f o r our system o f equations, we f i r s t make the assumption that the s o l i d s stream v e l o c i t y is piecewise constant f o r very small a x i a l s e c t i o n s of the r e a c t o r , i . e . , f o r the l o c a l i n t e g r a t i o n step. The s o l i d s v e l o c i t y s t i l l v a r i e s s i g n i f i c a n t l y w i t h i n the g a s i f i e r , but i t s change is assumed piecewise r a t h e r than continuous. I f this is done, then the s o l i d s stream energy balance dynamic equations can be r e w r i t t e n as

^

W

+ W

uS

at

8z

m

RHS

where RHS is the right-hand s i d e of the original equation. Since s o l i d s v e l o c i t y has been assumed piecewise constant, the t o t a l d i f f e r e n t i a l f o r the energy f l u x is given by

dz

3t

P u t t i n g the above two equations i n t o v e c t o r - m a t r i x n o t a t i o n , —

—

1

1

US

8ψε

RHS

~

dt

dz

di/;S 3z~

we have a system of simultaneous l i n e a r a l g e b r a i c equations in terms o f the f i r s t partial derivatives. The c h a r a c t e r i s t i c s o l u t i o n f o r these equations is obtained when the determinant of the matrix vanishes. From l i n e a r equation theory, the determinant must then a l s o vanish when the column v e c t o r on the right-hand s i d e o f the v e c t o r - m a t r i x equation is s u b s t i t u t e d f o r e i t h e r of the columns in the matrix on the l e f t . S u b s t i t u t i n g f o r the second column and s e t t i n g the determi nant t o zero gives

i_^=RHS. uS dt

Likewise, s u b s t i t u t i n g f o r the f i r s t

column gives

^=RHS. dz Both of these equations are o r d i n a r y d i f f e r e n t i a l equations f o r the energy f l u x along the s o l i d s v e l o c i t y c h a r a c t e r i s t i c curve. E i t h e r one can be used, but we chose t o implement the second



342

CHEMICAL REACTORS

equation f o r our method of c h a r a c t e r i s t i c s dynamic model. A s i m i l a r d e r i v a t i o n can be done to develop the corresponding s o l i d s stream mass balance equations. The i n t e g r a t i o n of these equations is c a r r i e d out f o r f i x e d time s l i c e s along the v a r i o u s c h a r a c t e r i s t i c curves w i t h i n the r e a c t o r . Since the s o l i d s v e l o c i t y v a r i e s down the r e a c t o r , this f i x e d time s l i c i n g gives the e f f e c t of a v a r i a b l e g r i d s t r u c t u r e with a v a r i a b l e number of nodes. The l o c a t i o n of the l a s t node is determined by the bottom of the r e a c t o r r a t h e r than by the time s l i c e . At any given time during a t r a n s i e n t response c a l c u l a t i o n , nodes can be e i t h e r added or subtracted to handle the changing c o n d i t i o n s w i t h i n the r e a c t o r . A time s l i c e of 1.5 minutes was used f o r most of the method of c h a r a c t e r i s t i c s c a l c u l a t i o n s . T h i s time s l i c e gave a t o t a l of 79 nodes f o r the base case initial c o n d i t i o n steady s t a t e . With the v a r i a b l e s o l i d s v e l o c i t y , the top t h i r d of the r e a c t o r had 13 nodes, the middle t h i r d had 15 nodes, and the bottom t h i r d had 51 nodes. T h i s is very s i m i l a r to the v a r i a b l e g r i d s t r u c t u r e used in the method of l i n e s c a l c u l a t i o n s . Numerical Considerations As the gas and s o l i d s streams move through the gasifier, d i f f e r e n t r e a c t i o n s are slowly s t a r t i n g and a b r u p t l y stopping as components disappear at d i f f e r e n t l o c a t i o n s in the v a r i o u s zones. While some r e a c t i o n s are proceeding v i g o r o u s l y , other r e a c t i o n s are j u s t s t a r t i n g at very low r a t e s . Extremely steep axial temperature and molar f l u x g r a d i e n t s are present in the burning zone area. A l s o , the zone l o c a t i o n s are c o n t i n u a l l y s h i f t i n g up and down the r e a c t o r during a t r a n s i e n t p e r i o d . T h u s , we would expect the g a s i f i e r d i f f e r e n t i a l equations to e x h i b i t a high degree of numerical s t i f f n e s s . To examine the extent of this problem, an eigenvalue a n a l y s i s was done f o r many d i f f e r e n t r e a c t o r t r a n s i e n t and steady s t a t e time p r o f i l e s . The range of s t i f f n e s s r a t i o s (absolute value of r a t i o of l a r g e s t to smallest eigenvalue, r e a l p a r t s only) observed f o r the d i f f e r e n t r e a c t o r zones was as f o l l o w s : d r y i n g zone d e v o l a t i l i z a t i o n zone g a s i f i c a t i o n zone burning zone ash zone

10E5 10E2 10E1 10E2 10E0

-

10E10 10E8 10E4 10E9 10E1.

T h i s i n d i c a t e s that the equations are extremely s t i f f with a c o n s t a n t l y changing s t i f f n e s s r a t i o throughout the r e a c t o r . Thus, the equations might be very s t i f f f o r a time, then moderately s t i f f f o r a while and then m i l d l y s t i f f at v a r i o u s other l o c a t i o n s w i t h i n the r e a c t o r . The eigenvalue mix was found to be very s i m i l a r f o r all of



17.

STiLLMAN


343

the p r o f i l e s . On the average, 36.54% of the eigenvalues were negative r e a l , 3.26% were p o s i t i v e r e a l , 1.39% were negative complex, O.16% were p o s i t i v e complex and the remaining 58.65% were zero. A l s o , the full Jacobian was s i n g u l a r at all locations w i t h i n the r e a c t o r . With the presence of p o s i t i v e r e a l eigenvalues and p o s i t i v e and negative complex eigenvalues, we would hope to f i n d a few o s c i l l a t o r y t r a n s i e n t responses, and indeed, some l i m i t c y c l e responses were obtained. The r a t i o of the imaginary to the r e a l part of the p o s i t i v e complex eigenvalue seemed to be the d e c i d i n g f a c t o r between steady s t a t e convergent or l i m i t c y c l e responses. T h i s r a t i o was never observed to go above about 5 f o r the steady s t a t e r e s u l t s , but it went as high as 22 f o r the l i m i t c y c l e runs. S t i f f n e s s is not a problem f o r the t r a n s i e n t part of the c a l c u l a t i o n s since in the t r a n s i e n t r e g i o n i n t e g r a t i o n step s i z e is l i m i t e d by accuracy r a t h e r than by s t a b i l i t y (9). Nonstiff i n t e g r a t i o n codes would be expected to perform b e t t e r than s t i f f codes in this case. For the time method of l i n e s and the method of c h a r a c t e r i s t i c s , the d i f f e r e n t i a l equations are n u m e r i c a l l y in a t r a n s i e n t s t a t e f o r the i n t e g r a t i o n code even though the r e a c t o r is in a steady s t a t e c o n d i t i o n . T h i s is due to the r e a c t i o n s t u r n i n g on and o f f at v a r i o u s l o c a t i o n s as the integration proceeds along the r e a c t o r . However, f o r the distance method of l i n e s technique where i n t e g r a t i o n is done at f i x e d distance nodes, s t i f f n e s s could be a problem as we approach steady s t a t e c o n d i t i o n s . Jacobian based s t i f f codes can not handle the i n t e g r a t i o n in i t s present form because the Jacobian is always s i n g u l a r . In many s e c t i o n s of the r e a c t o r , the number of a c t i v e d i f f e r e n t i a l equations is c o n t i n u a l ly i n c r e a s i n g or decreasing. Codes based on v a r i a b l e order multistep methods are very i n e f f i c i e n t under these conditions because they must be r e s t a r t e d whenever a change o c c u r s . This means that we would be using higher order code f o r some p a r t s of the r e a c t o r and lower order code in the more r a p i d l y changing sections. A f i x e d higher order method would probably be b e t t e r . A l s o , f o r m u l t i s t e p methods, the i n t e g r a t o r work- space has to be saved f o r each of the nodes. Another p o s s i b i l i t y would be to use a m u l t i r a t e code which i n t e g r a t e s each d i f f e r e n t i a l equation i n d i v i d u a l l y using d i f f e r e n t step s i z e s . O r a i l o g l u (16) and Gear (17) d i s c u s s this approach but t h e i r procedure uses m u l t i s t e p Jacobian methods which are not e f f i c i e n t f o r our system of equations. What we need is a f i x e d order s i n g l e step m u l t i r a t e method. We f i n a l l y decided to solve the problem by using a simple s i f t i n g procedure s i m i l a r to the one used by Emanuel and Vale (18). Before the i n t e g r a t o r was c a l l e d , the d e r i v a t i v e s were determined f o r the s o l i d s stream molar f l u x equations. Any f l u x equation that had a d e r i v a t i v e below the value of the s i f t paramet e r was considered to be i n a c t i v e at that time. The energy f l u x equation was always considered a c t i v e . Only the a c t i v e equations



344

CHEMICAL REACTORS

were then i n t e g r a t e d f o r the given time s t e p . T h i s procedure was repeated f o r each node at each i n t e g r a t o r time s t e p . In this way the s t i f f n e s s r a t i o could be reduced to l e s s than 100 in the d r y i n g zone, l e s s than 10 in the d e v o l a t i l i z a t i o n zone and to about 1 in the g a s i f i c a t i o n , b u r n i n g , and ash zones, while s t i l l maintaining c a l c u l a t e d r e s u l t s very c l o s e to those obtained f o r u n s i f t e d runs. Various i n t e g r a t i o n methods were t e s t e d on the dynamic model equations. They i n c l u d e d an i m p l i c i t i t e r a t i v e m u l t i s t e p method, an i m p l i c i t E u l e r / m o d i f i e d E u l e r method, an i m p l i c i t midpoint averaging method, and a modified d i v i d e d d i f f e r e n c e form of the v a r i a b l e - o r d e r / v a r i a b l e - s t e p Adams PECE formulas with local extrapolation. However, the best i n t e g r a t o r f o r our system of equations turned out to be the variable-step fifth-order Runge-Kutta-Fehlberg method. This e x p l i c i t method was used f o r all of the c a l c u l a t i o n s presented here. The initial c o n d i t i o n s w i t h i n the r e a c t o r needed to s t a r t the dynamic model c a l c u l a t i o n s were e s t a b l i s h e d by u s i n g the steady s t a t e model Ç 7 , 8) to c a l c u l a t e a f i r s t estimate. T h i s estimate was adjusted f o r use with the d i s t a n c e method of l i n e s , time method of l i n e s , and method of c h a r a c t e r i s t i c s programs by running the i n d i v i d u a l programs to a steady s t a t e c o n d i t i o n without changing input c o n d i t i o n s . The dynamic model c a l c u l a t i o n s are done in a two phase process. In the f i r s t phase, the s o l i d s stream dynamic equations are i n t e g r a t e d f o r the r e a c t o r while keeping gas stream c o n d i t i o n s constant. When necessary, intermediate gas stream values are obtained by i n t e r p o l a t i o n between storage nodes. These c a l c u l a t i o n s proceed down the r e a c t o r f o r the distance and time method of l i n e s , and back up the r e a c t o r f o r the method of c h a r a c teristics. In the second phase, the gas stream steady state equations are i n t e g r a t e d from the bottom to the top of the r e a c t o r while keeping the s o l i d s stream c o n d i t i o n s constant. Intermediate s o l i d s stream values are obtained by i n t e r p o l a t i o n between storage nodes. The numerical s t a b i l i t y requirement f o r the c o u p l i n g of the g a s - s o l i d s c a l c u l a t i o n s in the distance method of l i n e s model was estimated to be uS At/Az 1.67 which sets the time g r i d s i z e spacing.

based on the s p e c i f i e d

storage node



17.

STiLLMAN


345

An IBM 370/158 computer was used f o r all of the c a l c u l a t i o n s . As would be expected, the c a l c u l a t i o n speed was d i f f e r e n t f o r each of the three dynamic models. The distance method of l i n e s model ran at a speed of O.23 times r e a l time (4.3 times slower than r e a l time) which was very slow. For short time t r a n s i e n t c a l c u l a t i o n s , the speed could be increased to r e a l time speed, but long term numerical s t a b i l i t y r e q u i r e d the slower speed. The time method of l i n e s model ran 3.92 times f a s t e r than r e a l time. However, the c a l c u l a t e d t r a n s i e n t responses were i n c o r r e c t and the model could not be used f o r that purpose. For a time s l i c e of 1.5 minutes, the method of c h a r a c t e r i s t i c s model ran 1.56 times f a s t e r than r e a l time. I f the time s l i c e was increased to 5.0 minutes (fewer nodes), the speed increased to 4.75 times r e a l time but the gas stream accuracy was reduced. T h e r e f o r e , the 1.5 minute time s l i c e was used f o r the c a l c u l a t i o n s shown here. Simulation Results Based on the dynamic model presented in the previous s e c t i o n s , three computer programs were w r i t t e n to simulate moving bed g a s i f i e r dynamics using the method of c h a r a c t e r i s t i c s and both method of l i n e s techniques. Table I l i s t s the ash (Lurgi) g a s i f i er operating data f o r the base case initial c o n d i t i o n s . Except for the m u l t i p l e steady state runs, all of the step change response c a l c u l a t i o n s were made using these initial condition values. The proximate, u l t i m a t e and s i m u l a t i o n model a n a l y s i s of the Roland seam subbituminous c o a l used in the c a l c u l a t i o n s are given in S t i l l m a n (7). Figure 1 shows the e x i t gas temperature time response to a step change in c o a l moisture from 34.67 to 27.00 wt % f o r the three dynamic models. T h i s is a l a r g e step change i n v o l v i n g a s o l i d s m a t e r i a l wave moving through the r e a c t o r and it was i n t e n d ed to provide a severe t e s t f o r the three methods. The upper curve was c a l c u l a t e d by the method of c h a r a c t e r i s tics program and it e x h i b i t s a true l i m i t c y c l e or sustained o s c i l l a t i o n response (19). The middle curve was c a l c u l a t e d by the d i s t a n c e method of l i n e s program. The response is attenuated and s t r e t c h e d out. The f i n a l long term o s c i l l a t i o n s had random unequal periods and they were out of phase with the MOC r e s u l t s . The lower curve was c a l c u l a t e d by the time method of lines program. The initial p a r t of the response is s i m i l a r to the DMOL results but then the temperature i n c o r r e c t l y l e v e l s out to a steady s t a t e c o n d i t i o n . Thus, it was evident that the distance and time method of l i n e s techniques were not as accurate as the method of c h a r a c t e r i s t i c s procedure f o r c a l c u l a t i n g the g a s i f i e r step responses and they were d i s c a r d e d . All of the remaining c a l c u l a t i o n s were done by the MOC program. The s o l i d s temperature at the O.3 m (1 f t ) l e v e l w i l l be used to d i s p l a y the step change responses. That l o c a t i o n was chosen because it is i n i t i a l l y s l i g h t l y above the burning zone


346

CHEMICAL REACTORS


Table I.

Base Operating Data f o r Ash Discharge

Reactor Bed Height Reactor I n t e r n a l Diameter E x i t Gas Pressure Dry Coal Feed Rate Coal Moisture Content Dry Coal/Oxygen Ratio Steam/Oxygen R a t i o I n l e t Coal Temperature I n l e t Gas Temperature E x i t Gas Temperature E x i t S o l i d s Temperature S o l i d s Τ at O.3 m ( 1 f t )

2.74

(Lurgi) Reactor

m (9.00 ft)

3.70

m

2.84

MPa ( 2 8 . 0 atm)

(12.14 f t )

2 0 6 7 . 7 kg/h-m-m ( 4 2 3 . 5 l b / h - f t - f t ) 3 4 . 6 7 wt %

2 . 8 0 wt/wt 8 . 2 0 mol/mol 78.0 360.6

oc ( 1 7 2 . 4 °c

oF)

( 6 8 1 . 0 F)

275.6

oc ( 5 2 8 . 1 F)

372.8

oc ( 7 0 3 . 0 F)

912.6

oc ( 1 6 7 4 . 6 F)

Dry E x i t Gas mol % Hydrogen

Carbon Dioxide Carbon Monoxide Methane Nitrogen, Ammonia, Hydrogen S u l f i d e , Tar

40.83

31.47 14.53 10.88 2.29

300 MINUTES AFTER STEP CHANGE

600

Figure 1. Exit gas temperature response, coal moisture reduced to 27.00 wt %


17.

STiLLMAN



where l a r g e temperature f l u c t u a t i o n s can occur when zone s h i f t s back and f o r t h in the r e a c t o r .

347 the burning

I n l e t Gas Temperature. A few step response runs were made f o r changes in i n l e t gas ( b l a s t ) temperature. Only steady s t a t e convergent responses were observed and the f i n a l r e s u l t s a r e summarized in Table I I . F i g u r e 2 shows the s o l i d s temperature response when the i n l e t gas temperature is reduced from 360.6 oC (681.0 oF) to 305.0 oC (581.0 o F ) . This change puts l e s s heat i n t o the r e a c t o r which slows down the r e a c t i o n s in the g a s i f i e r above the burning zone. T h i s causes the burning zone t o s h i f t down the r e a c t o r a short d i s t a n c e producing the r a p i d temperature drop e x h i b i t e d in F i g u r e 2. Figure 3 shows the response when the i n l e t gas temperature is increased to 416.1 oC (781.0 oF). Since more heat is a v a i l a b l e in the g a s i f i e r , the r e a c t i o n r a t e s are increased, the burning zone s h i f t s slowly up the r e a c t o r a short d i s t a n c e , and the s o l i d s temperature a t the O.3 m ( 1 f t ) l e v e l g r a d u a l l y increases as the burning zone moves a l i t t l e c l o s e r . Coal Moisture. The r e s u l t s f o r changes in c o a l moisture content are given in Tables I I I , IV and V, and F i g u r e s 4 and 5. In a d d i t i o n to steady s t a t e convergent r e s u l t s , we now have underdamped (decaying o s c i l l a t i o n ) and l i m i t c y c l e responses. F i g u r e 4 shows the s o l i d s temperature response when the i n l e t c o a l moisture is reduced from 34.67 to 32.91 wt %. There is a 66 minute delay in the response which represents the transport time f o r the s o l i d s stream wave to reach theO.3m ( 1 f t ) l e v e l . As the burning zone begins to s h i f t up the r e a c t o r , the s o l i d s temperature peaks and s t a r t s back down during the next 50 minutes which is the time it takes the s o l i d s wave to move out o f the g a s i f i e r . The s o l i d s temperature o s c i l l a t e s a few more times as the l o c a t i o n o f the burning zone s h i f t s up and down before s e t t l i n g out to a steady s t a t e p o s i t i o n s l i g h t l y above i t s initial location. The underdamped steady s t a t e r e s u l t s f o r reducing the c o a l moisture content t o 31.05 wt % are given in the t h i r d data column of Table I I I . The s o l i d s temperature response s t a r t s out s i m i l a r to the F i g u r e 4 r e s u l t s except the temperature r i s e s to a higher peak value o f 1133 c (2072 © F ) . S t a r t i n g with an amplitude of about 5 0 °C ( 9 0 F), the temperature response then o s c i l l a t e s in a decaying manner f o r the next 2300 minutes u n t i l it reaches a f i n a l steady s t a t e value. The f i n a l l o c a t i o n of the maximum solids temperature is s l i g h t l y below the O.3 m ( 1 f t ) l e v e l . The limit c y c l e response caused by decreasing the c o a l moisture content to 27.00 wt % is shown in F i g u r e 5. Again there is a steep s o l i d s temperature r i s e as the s o l i d s m a t e r i a l wave reaches the O.3 m ( 1 f t ) l e v e l , but in this case, the decaying temperature o s c i l l a t i o n s only l a s t about 2 0 0 minutes before the g a s i f i e r s e t t l e s i n t o a sustained o s c i l l a t i o n mode with an e s t i -

American Chemical Society Library 1155 16th St. N. W. Fogler; Chemical Washington, D. Reactors C. 20038

ACS Symposium Series; American Chemical Society: Washington, DC, 1981.

348

CHEMICAL REACTORS


Table I I .

I n l e t Gas Temperature, F i n a l Steady State Results

I n l e t Gas Temperature

oC oF

305.0 581.0

388.3 731.0

416.1 781.0

E x i t Gas Temperature

oC °F

269.4 517.0

276.7 530.0

277.6 531.7

E x i t S o l i d s Temperature

oC oF

330.9 627.7

398.1 748.6

424.3 795.7

S o l i d s Τ at O.3 m (1 f t ) °C oF

778.7 1433.6

951.4 1744.6

986.4 1807.5

40.92 32.08 13.83 10.88 2.29

40.82 31.45 14.56 10.88 2.29

40.82 31.43 14.59 10.87 2.29

Dry E x i t Gas Hydrogen Carbon Dioxide Carbon Monoxide Methane N i t r o g e n , Ammonia, Hydrogen S u l f i d e , Tar

mol %

T,°C 900

800

700 0

Figure 2.


600

Solids temperature atO.3m, inlet gas temperature reduced to 305.0°C


17.


STiLLMAN

349

T, C

1900 \-


1000

INLET GAS Τ 416.1

C 600


Figure 3.

Solids temperature atO.3m, inlet gas temperature increased to 416.1 °C

Table I I I .

Coal Moisture Change, F i n a l Steady State R e s u l t s

Coal Moisture

wt %

33.80

32.91

31.05


oC oF

276.0 528.8

281.2 538.1

290.4 554.8


oC oF

375.1 707.1

372.1 701.7

368.1 694.6

S o l i d s Τ a t O.3 m (1 f t ) oC oF

838.6 1541.5

928.5 1703.3

1112.9 2035.2

41.06 31.91 13.85 10.90 2.28

41.08 31.95 13.79 10.90 2.28

41.11 32.02 13.69 10.90 2.28

Dry E x i t Gas Hydrogen Carbon Dioxide Carbon Monoxide Methane Nitrogen, Ammonia, Hydrogen S u l f i d e , Tar

mol %


CHEMICAL REACTORS

350

Table IV.

Coal Moisture


30.08

wt %, L i m i t Cycle R e s u l t s

oC °F


°C OF

S o l i d s Τ a t O.3 m ( 1 f t )oC


°F


Max

Min

297.4 567.3

291.7 557.1

367.3 693.1

366.4 691.6

1126.1 2058.9

1070.4 1958.8

41.00 31.78 14.06 10.88 2.28

41.22 32.27 13.31 10.92 2.28

mol %

Estimated L i m i t Cycle P e r i o d 81.4 minutes

Table V.

Coal Moisture

27.00

wt %, L i m i t Cycle Results Max

Min


c oF

306.9 584.4

303.3 577.9


c oF

367.2 693.0

366.6 691.9

(1 f t ) °c oF

1123.6 2054.5

1095.8 2004.5

S o l i d s Τ a tO.3m


mol % 41.25

32.30 13.27 10.91 2.27

41.36

32.57 12.86 10.94 2.27




STiLLMAN

τ, F

1800


1600

h


Figure 4. Solids temperature atO.3m, coal moisture reduced to 32.91 wt %

T,°F 1100

2000

1000

1700 300 MINUTES AFTER STEP CHANGE

Figure 5. Solids temperature atO.3m, coal moisture reduced to 27.00 wt %


352

CHEMICAL REACTORS


mated p e r i o d of 82.8 minutes. The l o c a t i o n of the s o l i d s peak burning zone temperature moves up and down in rhythm with the o s c i l l a t i o n s but always remains s l i g h t l y above the O.3 m (1 f t ) level. The decaying type response means that the l i m i t cycle s u r f a c e , which encloses an e q u i l i b r i u m or steady s t a t e p o i n t , was approached from the outside d i r e c t i o n . Steam Flow. The s team/oxygen feed r a t i o is an important v a r i a b l e f o r c o n t r o l l i n g the gas and s o l i d s temperatures in the burning zone and changing this r a t i o has a s i g n i f i c a n t e f f e c t on g a s i f i e r operation. The s i m u l a t i o n r e s u l t s f o r steam flow r a t e step changes are given in Tables VI-X and Figures 6, 7 and 8. Besides having steady s t a t e convergent, underdamped steady s t a t e , and l i m i t c y c l e responses, we a l s o have an underdamped l i m i t c y c l e response. Figure 6 shows the steady s t a t e convergent results for reducing the steam flow r a t e by 1.875%. Since the oxygen flow r a t e remains constant at its initial v a l u e , this step change reduces the steam/oxygen molar r a t i o from 8.20 to 8.05. The response shows a s m a l l almost l i n e a r temperature r i s e f o r about 40 minutes, and then a h e a v i l y damped o s c i l l a t o r y drop to the f i n a l steady s t a t e c o n d i t i o n . Even though reducing the steam flow r a t e r a i s e s the gas peak temperature, the i n c r e a s e in this run is not enough to overcome the e f f e c t of the reduced gas stream flow which puts l e s s heat i n t o the upper p a r t of the r e a c t o r . This s h i f t s the burning zone down the r e a c t o r and causes the f i n a l solids temperature at the O.3 m (1 f t ) l e v e l to drop by about 41 ©C (74 oF). The underdamped steady s t a t e r e s u l t s f o r a 2.25% decrease in steam feed r a t e are given in the t h i r d data column of Table V I . T h i s step change reduces the steam/oxygen molar r a t i o to 8.02. The underdamped response s t a r t s out very s i m i l a r to the F i g u r e 6 r e s u l t s and appears to be headed f o r a steady state condition. However, beginning at about the 305 minute p o i n t , the burning zone begins to move back up the r e a c t o r . Thus, during the p e r i o d of 305 to 365 minutes, the s o l i d s temperature at the O.3 m (1 f t ) level r i s e s from 868 c (1594 F) to 962 c (1763 F). The response then s e t t l e s back i n t o a decaying o s c i l l a t i o n mode u n t i l steady s t a t e c o n d i t i o n s are reached 445 minutes l a t e r . The underdamped l i m i t c y c l e r e s u l t s f o r a 2.50% decrease in steam feed r a t e (8.00 steam/oxygen molar r a t i o ) are given in Table VII. Once a g a i n , the response s t a r t s out l i k e the F i g u r e 6 response but then goes i n t o a decaying l i m i t c y c l e mode. However a g a i n , beginning at about the 340 minute p o i n t , the burning zone shifts back up the r e a c t o r . During the next 60 minutes, the s o l i d s temperature at the O.3 m (1 f t ) l e v e l goes from 860 °C (1580 ° F ) to 956 ©c (1752 F). The response then r e t u r n s to the decaying l i m i t c y c l e mode f o r the next 935 minutes u n t i l a very small s t a b l e l i m i t c y c l e is reached. The l i m i t c y c l e responses f o r a 5.0% decrease and a 5.0%


17.

STiLLMAN


Table VI.

353


Steam Feed Change, F i n a l Steady State Results

Steam Feed Rate Change

%

-1.25

-1.875

-2.25


oC oF

273.2 523.8

271.9 521.5

271.3 520.3

Exit Solids

°C oF

374.1 705.3

374.8 706.7

372.4 702.4

S o l i d s Τ a t O.3 m (1 f t ) oC oF

886.9 1628.5

872.3 1602.1

952.3 1746.2

40.66 31.36 14.74 10.94 2.30

40.58 31.31 14.84 10.97 2.30

40.53 31.28 14.90 10.99 2.30

Temperature


Table V I I .

mol %

Steam Feed Rate Change -2.5%, L i m i t

Cycle R e s u l t s


©c oF

Max 271.0 519.9

Min 270.8 519.5


c oF

372.8 703.0

372.7 702.8

S o l i d s Τ a t O.3 m (1 f t ) c oF

948.8 1739.8

945.6 1734.1

40.50 31.24 14.96 11.00 2.30

40.51 31.27 14.92 11.00 2.30


Estimated L i m i t

mol %

Cycle P e r i o d 82.0 minutes


354

CHEMICAL

Table V I I I .

Steam Feed Rate Change -3.75%, L i m i t ι

°F

Max 273.5 524.3

Min 267.7 513.8

oC op

375.2 707.3

372.6 702.6

S o l i d s Τ a t O.3 m (1 f t ) °C

977.8 1792.1

889.4 1632.9

40.15 30.76 15.74 11.03 2.32

40.53 31.62 14.47 11.08 2.30



°C

OF


REACTORS


mol %


Table IX.

Steam Feed Rate Change -5.0%, L i m i t Cycle Results


oC OF


oC °F

S o l i d s Τ a t O.3 m (1 f t ) oC OF


Max 271.7 521.1

Min 265.2 509.3

376.8 710.2

373.6 704.4

958.7 1757.6

860.1 1580.2

39.92 30.52 16.15 11.08 2.33

40.37 31.53 14.64 11.15 2.31

mol %




STILLMAN

Table X.

Steam Feed Rate Change +5.0%, L i m i t Cycle R e s u l t s


oC oF


°C OF

S o l i d s Τ at O.3 m (1 f t ) °C OF


355


Max 279.8 535.7

Min 275.3 527.6

374.7 706.5

372.9 703.3

844.5 1552.1

802.6 1476.7

41.46 31.93 13.68 10.66 2.27

41.66 32.34 13.05 10.69 2.26

mol %


1800

k

1600


Figure 6.

600

Solids temperature atO.3m, steam feed rate reduced 1.875 %


356

CHEMICAL REACTORS

T,°F

1800 \-


1600


Figure 7. Solids temperature atO.3m, steam feed rate reduced 5.0 %

1600

Y

1400


Figure 8. Solids temperature atO.3m, steam feed rate increased 5.0 %


17.

STiLLMAN


357


i n c r e a s e in steam feed flow r a t e are given in F i g u r e s 7 and 8 · For these r u n s , the steam/oxygen molar r a t i o s are 7.79 and 8.61, respectively. Both of the curves show a c l a s s i c decaying l i m i t c y c l e form which i n d i c a t e s that the l i m i t c y c l e surface was approached from the o u t s i d e . E x i t Gas P r e s s u r e . G a s i f i e r pressure has an important e f f e c t on the composition of the e x i t raw gas and on the o p e r a t i o n of the reactor. At normal o p e r a t i n g p r e s s u r e , increases in reactor pressure i n c r e a s e the methane and carbon d i o x i d e content of the raw gas and decrease the hydrogen and carbon monoxide content. Decreases in pressure have the opposite e f f e c t . The r e a c t o r pressure step change r e s u l t s are given in F i g u r e s 9 and 10, and Tables X I , XII and X I I I . Figure 9 shows the steady s t a t e convergent results for reducing the e x i t gas pressure from 2.84 MPa (28 atm) to 2.53 MPa (25 atm). The pressure decrease lowers the gas d e n s i t y but increases the gas v e l o c i t y s i n c e the molar feed r a t e of the i n l e t gas remains constant. T h i s puts l e s s heat i n t o the r e a c t o r which lowers the e x i t gas temperature and causes the burning zone to s h i f t down the r e a c t o r a short d i s t a n c e . This is r e f l e c t e d in the F i g u r e 9 s o l i d s temperature response which s t a r t s out with a short f a l s e temperature r i s e but then e x h i b i t s a r a p i d temperature drop to a f i n a l steady s t a t e c o n d i t i o n (non-minimum phase response). Figure 10 shows the l i m i t c y c l e response produced when the e x i t gas pressure is increased to 3.45 MPa (34 atm). T h i s change decreases the gas v e l o c i t y and s h i f t s the burning zone up the r e a c t o r a short distance where, in this case, it o s c i l l a t e s up and down with a p e r i o d of about 82.6 minutes. In F i g u r e 10, the s o l i d s temperature response at the O.3 m (1 f t ) l e v e l begins with a short f a l s e temperature drop and then q u i c k l y r i s e s and goes i n t o an expanding o s c i l l a t i o n mode before s e t t l i n g into a stable l i m i t c y c l e response. The expanding o s c i l l a t i o n s i n d i c a t e that the l i m i t c y c l e surface was approached from the i n s i d e . The Table XII c a l c u l a t i o n had a more pronounced expanding o s c i l l a t i o n phase as it took 846 minutes to reach a s t a b l e l i m i t c y c l e response. No experimental g a s i f i e r data was found to v e r i f y any of the s i m u l a t i o n model r e s u l t s . However, l i m i t c y c l e responses have been experimentally observed f o r the pressure in combustion chambers and b o i l e r s (20). This lends credence to our c a l c u l a t e d r e s u l t s s i n c e c o a l combustion is an important f a c t o r in g a s i f i e r operation. M u l t i p l e Steady S t a t e s . Not only can the n o n l i n e a r i t y of our dynamic model equations produce l i m i t c y c l e responses, it can give r i s e to m u l t i p l e steady s t a t e c o n d i t i o n s in which the same set of o p e r a t i n g parameters can produce d i f f e r e n t r e a c t o r p r o f i l e s . By accident, three of these m u l t i p l e steady s t a t e responses were obtained and they are summarized in Table XIV. To check the v a l i d i t y of the pressure step change r e s u l t s ,



CHEMICAL REACTORS

1400

h


600

Figure 9. Solids temperature atO.3m, exit gas pressure reduced to 2.53 MPa


17.

STiLLMAN

Table X I .


E x i t Pressure Change, F i n a l Steady State Results

E x i t Gas Pressure


MPa atm

2.53 25.0

3.14 31.0

oC

268.2 514.7

278.2 532.7

384.1 723.3

370.1 698.2

785.2 1445.4

992.8 1819.0

41.23 31.00 15.00 10.50 2.27

40.44 32.30 13.58 11.37 2.31

OF

Exit Solids

Temperature

oC OF

S o l i d s Τ at O.3 m (1 f t ) oC Downloaded by CORNELL UNIV on October 15, 2016 | http://pubs.acs.org Publication Date: September 21, 1981 | doi: 10.1021/bk-1981-0168.ch017

°F


Table X I I .

359

E x i t Pressure

mol %

3.24 MPa, L i m i t Cycle Results


°C °F

Max 280.7 537.2

Min 275.6 528.1

Exit Solids

°C oF

370.3 698.6

368.9 696.0

S o l i d s Τ at O.3 m (1 f t ) °C oF

1067.7 1953.9

978.9 1794.0

40.10 32.16 13.90 11.52 2.32

40.35 32.80 12.95 11.59 2.31

Temperature


mol %



360

CHEMICAL REACTORS


Table X I I I .

E x i t Pressure 3.45 MPa, L i m i t Cycle R e s u l t s


oC oF

Max 283.6 542.4

Min 276.4 529.6

Exit Solids

oC oF

369.2 696.5

367.3 693.1

S o l i d s Τ at O.3 m (1 f t ) oC oF

1126.8 2060.3

1026.6 1879.9

39.66 32.51 13.61 11.88 2.34

39.98 33.46 12.21 12.03 2.32

Temperature


mol %


Table XIV.

E x i t Pressure,

M u l t i p l e Steady State R e s u l t s

E x i t Gas Pressure

MPa atm


°C OF

Exit Solids

Temperature

oC OF

S o l i d s Τ at O.3 m (1 f t ) oC OF


2.84 28.0

2.84 28.0

275.6 528.1

275.7 528.2

270.9 519.7

372.0 701.6

370.4 698.8

379.0 714.2

918.3 1684.9

1002.6 1836.6

782.7 1440.8

40.83 31.47 14.53 10.88 2.29

40.83 31.47 14.53 10.88 2.29

41.03 31.86 13.93 10.90 2.28

2.84* 28.0

mol %

*unstable steady s t a t e r e s u l t s



17.

STILLMAN


361

initial conditions d i f f e r e n t from the base case values were used f o r some o f the t r a n s i e n t c a l c u l a t i o n s . For example, i f the f i n a l steady s t a t e c o n d i t i o n s produced by i n c r e a s i n g the e x i t gas pressure from 2.84 MPa (28 atm) to 3.14 MPa (31 atm) a r e used as the initial c o n d i t i o n s , and the e x i t gas pressure is stepped back to i t s original value o f 2.84 MPa (28 atm), we o b t a i n a steady s t a t e convergent response which r e t u r n s to the original initial c o n d i t i o n . Likewise, i f we step the e x i t gas pressure up to 3.45 MPa (34 atm), we o b t a i n the same l i m i t c y c l e response that we obtained in stepping from 2.84 MPa (28 atm) to 3.45 MPa (34 atm). Next we used the 3.45 MPa (34 atm) l i m i t c y c l e r e s u l t s a t the 810 minute p o i n t f o r our initial c o n d i t i o n s . I f we step the e x i t gas pressure down to 3.14 MPa (31 atm), we o b t a i n the same underdamped steady s t a t e r e s u l t that we obtained by stepping from 2.84 MPa (28 atm) to 3.14 MPa (31 atm). However, i f we step the e x i t gas pressure down to 2.84 MPa (28 atm), we do not r e t u r n to the original initial condition state. Instead, we o b t a i n the m u l t i p l e steady s t a t e r e s u l t given in the t h i r d data column of Table XIV. The e x i t gas temperature, the s o l i d s temperature at the O.3 m (1 f t ) l e v e l , and the raw gas carbon monoxide content are lower than the original v a l u e s . The e x i t s o l i d s temperature, and the hydrogen and carbon d i o x i d e content o f the raw gas a r e higher. This c o n d i t i o n is caused by the burning zone not being able to move back up the r e a c t o r a f t e r i t s initial downward s h i f t . We next t r i e d stepping the e x i t gas pressure down to 3.14 MPa (31 atm), and then a f t e r 270 minutes, stepping on down to 2.84 MPa (28 atm). These r e s u l t s are given in the f i r s t two data columns of Table XIV. A f t e r about 370 minutes, the r e a c t o r reached the spurious steady s t a t e c o n d i t i o n s given in column one, which are very c l o s e to the original initial c o n d i t i o n v a l u e s . However, the r e a c t o r only remained there f o r about 20 minutes, and then spent the next 330 minutes moving to the m u l t i p l e steady s t a t e r e s u l t given in column two. The spurious steady s t a t e r e s u l t s r e p r e s e n t ed an unstable m u l t i p l e steady s t a t e c o n d i t i o n . The column two r e s u l t s show that even though the s o l i d s temperature p r o f i l e down the r e a c t o r is s l i g h t l y d i f f e r e n t from the original p r o f i l e , the e x i t gas composition remains the same (to two decimal p l a c e s ) . F i n a l l y , the e x i t gas pressure was stepped down from 3.45 MPa (34 atm) to 3.24 MPa (32 atm), 3.04 MPa (30 atm), and 2.84 MPa (28 atm) in 270 minute increments. T h i s procedure gave a steady s t a t e convergent response which corresponded to the original initial condition. B i f u r c a t i o n . B i f u r c a t i o n r e f e r s to the switching or branching from one type of response behavior to another as a parameter passes through a critical v a l u e . F o r us, when the parameter is below the critical v a l u e , our step change responses are in the steady s t a t e convergent r e g i o n . When they are above the critical value, our responses are in the l i m i t c y c l e r e g i o n . As the parameter values approach the critical p o i n t , we enter a t r a n s i -



362

CHEMICAL REACTORS

t i o n r e g i o n . Just below the critical p o i n t , even though we s t i l l have steady s t a t e convergent responses, the steady s t a t e is only reached a f t e r a long p e r i o d of decaying o s c i l l a t i o n . The c l o s e r we approach the critical p o i n t , the longer the decaying o s c i l l a t i o n time becomes. Just above the critical p o i n t , we have long periods of e i t h e r decaying or expanding o s c i l l a t i o n s which f i n a l l y r e s u l t in s t a b l e l i m i t c y c l e responses. Thus, the approach to the critical point from e i t h e r above or below is asymptotic. Table XV provides a summary of the c o a l moisture, steam feed r a t e and e x i t gas pressure t r a n s i e n t response runs showing the time r e q u i r e d to reach the given c o n d i t i o n . I t provides a rough estimate f o r the values of the b i f u r c a t i o n p o i n t s f o r these runs. Thus, the b i f u r c a t i o n p o i n t f o r the c o a l moisture step change runs l i e s between 30.08 and 31.05 wt % moisture. For the steam feed r a t e changes, it l i e s between -2.25% and -2.50%, and f o r the e x i t gas pressure, it is bracketed by the 3.14 MPa (31 atm) and 3.24 MPa (32 atm) v a l u e s . Summary The c o n t i n u i t y equations f o r mass and energy were used to d e r i v e an a d i a b a t i c dynamic plug flow s i m u l a t i o n model f o r a moving bed c o a l g a s i f i e r . The r e s u l t i n g set of h y p e r b o l i c p a r t i a l d i f f e r e n t i a l equations represented a s p l i t boundary-value problem. The inherent numerical s t i f f n e s s of the coupled g a s - s o l i d s equat i o n s was handled by removing the time d e r i v a t i v e from the gas stream equations. This converted the dynamic model to a set of p a r t i a l d i f f e r e n t i a l equations f o r the s o l i d s stream coupled to a set of o r d i n a r y d i f f e r e n t i a l equations f o r the gas stream. The method of c h a r a c t e r i s t i c s , the d i s t a n c e method of l i n e s (continuous-time d i s c r e t e - s p a c e ) , and the time method of lines (continuous-space d i s c r e t e - t i m e ) were used to s o l v e the solids stream p a r t i a l d i f f e r e n t i a l equations. Numerical s t i f f n e s s was not considered a problem f o r the method of c h a r a c t e r i s t i c s and time method of l i n e s c a l c u l a t i o n s . For the d i s t a n c e method of l i n e s , a p o s s i b l e numerical s t i f f n e s s problem was solved by using a simple s i f t i n g procedure. A variable-step fifth-order Runge-Kutta-Fehlberg method was used to i n t e g r a t e the d i f f e r e n t i a l equations f o r both the s o l i d s and the gas streams. Step change dynamic response runs revealed that the d i s t a n c e and time method of l i n e s techniques were not as accurate as the method of c h a r a c t e r i s t i c s procedure for calculating gasifier t r a n s i e n t s . Therefore, these two techniques were discarded and the remaining c a l c u l a t i o n s were all done u s i n g the method of characteristics. The n o n l i n e a r dynamic g a s i f i e r model produced a wide v a r i e t y of t r a n s i e n t response types when subjected to step changes in operating c o n d i t i o n s . However, stepping the i n l e t gas (blast) temperature up and down by as much as 55.6 °C (100 °F) only gave steady s t a t e convergent responses. Varying the feed c o a l moisture



17.

STILLMAN

Table XV.


Summary of B i f u r c a t i o n Response Time Data Minutes 375 480 2430 324 315

Type Steady State Steady State Decay Steady State L i m i t Cycle L i m i t Cycle

Steam Feed % Change -1.25 -1.875 -2.25 -2.50 -3.75 -5.00 +5.00

270 375 810 1335 345 300 210

Steady State Steady State Decay Steady State Decay L i m i t Cycle L i m i t Cycle L i m i t Cycle L i m i t Cycle

Pressure MPa(atm) 2.53(25) 3.14(31) 3.24(32) 3.45(34)

270 1890 846 351

Steady State Decay Steady State Expand L i m i t Cycle L i m i t Cycle

Coal Moisture wt % 33.80 32.91 31.05 30.08 27.00



364

CHEMICAL REACTORS

content in the range of 34.67 to 27.00 wt % gave steady s t a t e convergent, underdamped (decaying o s c i l l a t i o n ) , and l i m i t cycle (sustained o s c i l l a t i o n ) responses. I n c r e a s i n g and decreasing the i n l e t steam flow r a t e by up to 5 % produced steady s t a t e conver gent, underdamped steady s t a t e , l i m i t c y c l e , and decaying o s c i l l a t i o n phase l i m i t c y c l e responses. P e r t u r b i n g the e x i t gas p r e s sure in the range of 2.53 MPa (25 atm) to 3.45 MPa (34 atm) gave non-minimum phase steady s t a t e convergent, underdamped steady s t a t e , l i m i t c y c l e , and expanding o c i l l a t i o n phase l i m i t cycle responses. Thus, in these c a l c u l a t i o n s the l i m i t c y c l e surface was approached from both the i n s i d e and the outside d i r e c t i o n s . Changing the e x i t gas pressure a l s o gave three m u l t i p l e steady s t a t e responses in which the same set of operating parame ters produced d i f f e r e n t reactor p r o f i l e s . Finally, a rough estimate f o r the l o c a t i o n of the b i f u r c a t i o n p o i n t s was given f o r the c o a l moisture, steam feed r a t e , and e x i t gas pressure t r a n s i e n t response r u n s .

Nomenclature aij AGS cpG cpS C CG CS CjG CjS F FjG FjS hGS ΔΗί i j k L Ρ ri t t At TG TS uG uS ζ Δζ

s t o i c h i o m e t r i c c o e f f i c i e n t component j in r e a c t i o n l o c a l g a s - s o l i d s heat t r a n s f e r area/volume r a t i o gas molar heat c a p a c i t y s o l i d s molar heat c a p a c i t y concentration variable t o t a l gas c o n c e n t r a t i o n t o t a l solids concentration component j c o n c e n t r a t i o n in gas component j c o n c e n t r a t i o n in s o l i d s molar f l u x v a r i a b l e gas molar f l u x f o r component j s o l i d s molar f l u x f o r component j l o c a l g a s - s o l i d s heat t r a n s f e r c o e f f i c i e n t heat of r e a c t i o n i g a s i f i e r r e a c t i o n index number g a s / s o l i d s component index number d i s t a n c e node index number top of r e a c t o r node index number absolute pressure r a t e of r e a c t i o n i r e a l time time node index number time g r i d spacing gas absolute temperature s o l i d s absolute temperature l o c a l gas v e l o c i t y local solids velocity d i s t a n c e from r e a c t o r bottom l o c a l d i s t a n c e g r i d spacing


i

17.

STILLMAN


365

Greek L e t t e r s |G φε i|>G

gas energy d e n s i t y s o l i d s energy d e n s i t y gas energy f l u x s o l i d s energy f l u x

Literature

Cited


1.

Yoon, H.; Wei, J.; Denn, M. "Modeling and Analysis of Moving Bed Coal Gasifiers"; EPRI Report No. AF-590, Vol 2, 1978. 2. ___________. "Transient Behavior of Moving Bed Coal Gasifica tion Reactors"; 71st Annual AIChE Meeting, Miami, 1978. 3. _____________. AIChE J. 1979, 25, (3), 429. 4. Hsieh, B. C. B.; Ahner, D. J.; Quentin, G. H. in Vogt, W. G.; Mickle, Μ. Η., Eds.; "Modeling and Simulation, Vol 10, Pt 3, Energy and Environment"; ISA: Pittsburgh, 1979; p 825. 5. Daniel, K. J. "Transient Model of a Moving-Bed Coal Gasifi er"; 88th National AIChE Meeting, Philadelphia, 1980. 6. Wei, J. "The Dynamics and Control of Coal Gasification Reac tors"; Proceedings JACC, Vol II; ISA: Pittsburgh, 1978; p 39. 7. Stillman, R. IBM J. Res. Develop. 1979, 23, (3), 240. 8. _____________. "Simulation of a Moving Bed Gasifier for a Western Coal Part 2: Numerical Data"; IBM Palo Alto Scientific Center Report No. G320-3382, 1979. 9. Shampine, L. F.; Gear, C. W. SIAM Review 1979, 21, (1), 1. 10. Carver, M. B. J. Comp. Physics 1980, 35, (1), 57. 11. Carver, M. B.; Hinds, H. W. Simulation 1978, 31, (2), 59. 12. Vichnevetsky, R. Simulation 1971, 16, (4), 168. 13. Buis, J. P. Simulation 1975, 25, (1), 1. 14. McAvoy, T. J. Simulation 1972, 18, (3), 91. 15. Abbott, M. B. "An Introduction to the Method of Characteris tics"; American Elsevier: New York, 1966. 16. Orailoglu, A. "A Multirate Ordinary Differential Equation In tegrator"; University of Illinois Department of Computer Science Report No. UIUCDCS-R-79-959, 1979. 17. Gear, C. W. "Automatic Multirate Methods for Ordinary Differ ential Equations"; University of Illinois Department of Computer Science Report No. UIUCDCS-R-80-1000, 1980. 18. Emanuel, G.; Vale, H. J. in Bahn, G. S., Ed.; "The Perfor mance of High Temperature Systems, Vol 2"; Gordon and Breach: New York, 1969; p 497. 19. Bailey, J. E. in Lapidus, L.; Amundson, N. R., Eds.; "Chemi cal Reactor Theory: A Review"; Prentice-Hall: Englewood Cliffs, 1977; p 758. 20. Friedly, J. "Dynamic Behavior of Processes"; Prentice-Hall: Englewood Cliffs, 1972; p 516. RECEIVED June 3, 1981.


Moving Bed Coal Gasifier Dynamics Using MOC and MOL Techniques

Recommend Documents