Heat Transfer in Packed Bed Heat Recuperators—I. A Precision Numerical Analysis

Heat Transfer in Packed Bed Heat Recuperators—I. A Precision Numerical Analysis

Heat Transfer in Packed Bed Heat Recuperators-I. A Precision Ibmerical Analysis* by STANLEY H. JURY Department of Chemical Albuquerque, New Mexico ...

774KB Sizes 0 Downloads 14 Views

Heat Transfer in Packed Bed Heat Recuperators-I. A Precision Ibmerical Analysis* by

STANLEY

H. JURY

Department of Chemical Albuquerque, New Mexico and

MARIO

Department Mexico

Engingeering,

University of New Mexico,

c. BERBANO

of Engineering,

University of New Mexico,

Albuquerque,

New

An extensive literature survey is reported that verifies the need for a more precise representation of the heat transfer film coefficient in packed beds. A precision numerical analysis scheme is developed for a much more general model than previously reported. Application of the results to the experimental data will be reported in Part II of the paper.

ABSTRA~:

Norution

a, c PG CIT CS 4 DO hR

h hw &, &I &I, Rso t z

El9

PG

PS PT TA

62,

E3,

E4

interstitial area per unit volume of bed (ft2/ft3) gas heat capacity (Btu/lb moleoF) heat capacity of plastic tube (Btu/lb”R) solid heat capacity (Btu/lb”R) inside diameter of tube wall (ft) outside diameter of plastic tube (ft) coefficient of radiation heat transfer from tube wall (Btu/hrft’“R) gas phase heat transfer film coefficient (Btu/hrft“‘R) heat transfer film coefficient from gas phase to tube wall (Btu/hrft2”R) inside radius of plastic tube and aluminum casing, respectively (ft) radius of inner and outer radiation shield, respectively (ft) time (hr) interstitial gas velocity (ft/hr) distance from bed entrance (ft) external porosity emissivity of plastic tube, inner and outer radiation shields, and aluminum casing, respectively gas density (lb moles/ftq solid density (lb/ft”) plastic tube density (lb/ft’) temperature of ambient air (OR)

* This work was supported by the National GK-31078 and Amendment No. 1.

285

Science Foundation

under Grant NSF

S. H. Jury and M. C. Berbano 70 (0 TG

7s * 7, TW

bulk temperature of gas at z and initial bed temperature (“R) solid temperature at r and t (“R)

t (“R)

gas phase temperature at gas-solid interface (“R) tube wall temperature at L and t (“R)

I. Introduction

The so-called simple or Schumann model (1) has been the primary basis of analysis of the heat recuperator. It is assumed that the velocity and thermal properties of the fluid are constant; it also neglects the existence of thermal gradients within the low thermal conductivity solid particles. The fluid velocity varies along the length of the bed. In addition, as the particle size increases radial temperature gradients will increase and intraparticle conduction can henceforth become a very significant factor. It is the ultimate purpose of this investigation to determine the functional relation of the heat transfer film coefficient that is consistent with a solution of the heat recuperation problem which will involve a kinetic mechanism that considers both intraparticle conduction and variable velocity down the bed. Experimental time-temperature breakthrough profiles subsequent to a step change in the inlet fluid temperature will be compared (in Part II) with the theoretically predicted in order to check the validity of the kinetic model. Here we report a precision numerical analysis scheme which enables one to determine by a least square method a functional relation for the heat transfer film coefficient without trial and error bed calculations. Ultimately, it is hoped that this type of analysis for the pure heat recuperation case can be extended to the percolation operations in general. IL Historical Background

of the Heat Recuperator

In 1926, Anzelius (2) and a year later Nusselt (3) reported analytical solutions to the basic recuperator problem. In 1929, Schumann (1) obtained a reduced solution which was first used by Furnas (4) in 1930 in determining the coefficient of heat transfer between a flowing gas and a bed of granular particles. More recently (5) the work of Schumann was extended to show that for longer times the same form of solution applied to internal diffusion as well as gas film diffusion (17). Due to the non-linear nature of the heat recuperator one must resort to numerical solutions in precision work. Several of these have been reported (24,27,14,15,20) and all leave something to be desired. A novel numerical analysis scheme was reported (16) for the solution of systems of equations describing simultaneous heat and mass transfer in a packed bed of sorbent granules. Adiabatic conditions and variability of the gas velocity and the physical properties with temperature and concentration were all considered. The important problem of stability was solved for the solution of mixed partial differential equations for the bed. This scheme of numerical analysis has not previously been applied to a heat recuperator as such. 286

Journal of The Franklin

Institute

Heat Transfer in Packed Bed Heat Recuperators--I Steady state experiments (27,ll) involve generating heat in the bed and feeding cold gas until a steady state prevails. Kudryavtsev et al. (19) found that steady and transient experiments did not lead to the same film coefficient for comparable conditions. Lijf and Hawley (21) emphasized that this should not be so if data were properly analyzed. The present state of affairs is somewhat less than satisfactory and in the present work we hope to improve upon this situation. III. Theoretical Treatment of the Heat Recuperator For the case of a hot dry gas stream passing through an initially cold bed of granular solids, the physical system, in its most general form, is quite complicated. The model

A major objective here is to include among other things the effect of internal heat conduction in the low thermal conductivity particles and variable velocity down the bed. Our simplifying assumptions are (1) constant physical properties of solid; (2) negligible interparticle radiation and conduction effects (3) onedimensional plug flow; and (4) negligible axial conduction in the gas phase compared to heat carried in gas by bulk flow. Recuperation

equations to include tube wall and radiation loss effects

A sketch of the packed bed is shown in Fig. 1. The thermal balance for the gas phase previously obtained (16) for adiabatic operation can be modified to omit adsorption heat and include heat lost to the plastic tube wall as follows:

:zpdstic holder

Absorbent Solids-

2

f AZ

,

1 2

2=0

FIG. 1. Heat and mass balance

Vol

303, No. 3, March 1977

for the bed.

287

S. H. Jury and M. C. Berbano

The h in Eq. (1) can be represented Dittus-Boelter type, which is h = V2(1)$

as a functional

relation

(2)

(NRc)v2(2)(NpJv2(3) P

where NRe respectively, conductivity parameters.

of the

and NPr are the dimensionless Reynolds and Prandtl numbers, Dp is the average particle diameter in feet, k~ is the gas thermal in Btu/hrft”R, and V2(1), V2(2) and V2(3) are dimensionless The wall heat transfer coefficient hw = CWZ

(N&0.83

was taken from the literature (4). The overall gas phase mass balance is aPG

d(upo)._O

at+

a2

(4)

The thermal behavior of the solid spherical particles is described by

(5) r = radial distance from center of sclid (ft) kS = average solid thermal conductivity

(Btu/hrft”R)

The thermal balance at the gas-solid interface is +h($-TG)=O

(6)

where R is the particle radius in feet. If one assumes ideal gas behavior and constant pressure operation P = &-JTG.

(7)

Here, P is total pressure in lb/in*abs, and R is the gas law constant in (lb/in2)(ft3)/lb moleoR. In addition to the foregoing equations one must also specify the initial and boundary conditions such as 7G = 76 = u= 7s =

measured feed gas temperature constant constant constant

7s = 7:

(87&r) = 0 Equation (6) 288

at at at at at at at

z=O, t>O r=O, z>O z=O, t>O t=O, r>O r=R, t>O r=O,

t>O

r= R, t>O.

(8) (9) (10) (11)

(12) (13) (14)

Journal of The Franklin Institute

Heat Transfer in Packed Bed Heat Recuperators-l

In addition to heat storage in the thbe that contains the bed, one should also consider radiation heat loss from the bed via the wall to the vacuum insulation. A heat balance on an elemental section of the tube will consequently yield: arw 4hwDI at=p,~n(&-D',)

4hRDo

(%-7w)--

p&q-(D;-

0:)

(7w - 7A) 7w=70(1),

t=o,

z>o.

(15)

The equation for the radiation coefficient hR is hR = 1.712~ 10-gCR(&+r$)(rw+rA).

(lo)

Here, CR is a constant

(17) From the foregoing mathematical description of the proposed model one finds that there are five primary dependent variables: (1) bulk gas temperature 76, (2) bulk velocity of the gas V, (3) bulk gas density 76, (4) plastic wall temperature TV, and (5) adsorbent solid temperature TV. The first four variables are functions only of time and position along the bed. The last variable Ts, is a function of time and radial position within the solid particle. Secondary variables, which are functions of the primary variables, are the physical and thermal properties of the gas. Difference quotient approximations

The bed equations, i.e. Eqs. (1) and (4), which are first-order hyperbolic-type, can be generally represented thus

+-y+A4$= u(x, t)

and of

(18)

where u is a dependent variable at time t and position x, and M is the variable coefficient. The differencing method can be simply illustrated through the special case, i.e.

au

z =

t-

by using the power series u(x, t)=u(x, in the integration

t1)+

I

(1-c)7+c

(19)

(2 )I2

{4x, f2) - 4x, t1))

(20)

of Eq. (19) with respect to time thus

4% t2)- 4% td At p

Vol. 303. No. 3, March 197,

tl

ub, t)

=)+&;

= 4x9 td+O--P)U(X, t2)

(21)

o
S. H. Jury and M. C. Berbano Based on a stability analysis it has been shown (16)that unconditional stability prevails. Application of method to discrete heat recuperator model

The gas phase thermal balance could have been written for a finite bed length AZ in which the gas phase is assumed mixed to uniformity. In this case, one obtains a&r+

1, t)7&+

1, t)

at

=z

1

Id& tM4

-u(l+

thiu,

1, t)po(l+

+-${&I+

1, &(I+

1, t)-7&+

--$+r+1, PG

t) 1, t)) 1, t)}

t)-r&+1,

t)

(22)

I

where I refers to the bed increment prior to the current bed increment I+ 1 being computed. The overall gas phase mass balance (Eq. 4) can likewise be written as

If one assumes ideal gas behavior and constant pressure operation P=&o(I+i,

t)r&+i,

so that

t)=constant

(24)

and, for purposes of scaling sets G(r+ 1, t) = u(I+ 1, t)p&+

1, t)

(25)

then for current time t = t2 the thermal balance (Eq. 22) reduces to C(I+ 1,

tdTG(r+

1,

td=

w

h)TG(&

f2)

-y{TG(I+i,

t&&I+i,

t2)}

PG

-E{TG(r+i, PG

t2)-TW(r+i,

t2)}.

(26)

I

Following second-order interpolation and integration, the density equation (Eq. 23) can be multiplied by r&I+ 1, tz). This will yield pG(r+

1,

f2)TG(r+

1,

thPG(r+

1,

thG(r+

=p-$G(& td-G(r+i,

290

1,

tZ)

t1))7G(r+i,

t2)

Journal of The Franklin

Institute

Hear Transfer in Packed Bed Heat Recuperators-I

Equation

(26) is now substituted

PG(0%U)-PoU+

I, fJ%(I+

in this multiplied density equation giving

I, f2)

=pE(RR0)7G(I+l,

tz)+(l-@)${G(I,

-WI,

&c(I+l,

h)-- u{T:(l+

f2bGt6

eC

fz)

1, Q--7&+

1,

tz)}

PO

+E

imu+1, tz)-- nvu+L12H)

(28)

where RR0 = G(I, tl)- G(I+ 1, tl) PGUMI)

= p&I+ 1,

t2bGu+

1,

t*)

(gas law)

where ~~(1) and ~~(1) are the initial bed density (lb/ft3) and initial bed temperature (OR), respectively. Before obtaining an explicit expression for T&I+ 1, t2), the wall temperature T&I+ 1, f2) must first be determined. If one also applies second-order interpolation followed by integration to Eq. (15), the result can be solved for T&I+ 1, tz) to give T&+

1,

t2) =

W(1+1,

tl)+To++(I-p)T’%(I+I,

tz)+YRQ

(29)

1+(1-p)(y+?R)

where 4Drhw A\t ‘=fiC&D:-D:) 4DohR At ‘YR= pTC,(D; - D’,) Yo=p{Y%(I+I,

tl)-(Y+YR)%‘(l+I,

tl)).

(30)

One can set p=o and define the normalized

(31)

temperatures

TG(I+

1, t2) =

TW(I+

1, f2) =

b)-%(1)

%(~-dI) %(I+ T&f)

TA _ -

?+(~+l,

1, fd-R&) - 76 (1)

TA-TG(I) 7~

@f)

Vol. 303. No. 3. March

1977

(33) (34)

-

TG (1)

r Ts(b, K)- V?(I) UT( tz, K) = z To(I)

Finally, the ~~(I-tl,

(32)

(35)

t2) in Eq. (29) is eliminated by means of Eq. (33). The 291

S. H. Jury and M. C. Berbano

resulting equation produces a result solvable for TG(I+ 1, tz) via the simplifying relations (31-35) as follows: TG(I+

1, tz) = QA - UT(t2, l)+ QB

(36)

where 1

QA=

(37)

4(1+yR&W

QD+

+

1

&DI(~+~+YR)

{TW(I+

1,

tl> +

yRTA)

,,=76&fhG(I)

TG(I)

.

With TG(I+l, tz) at hand one can compute p&+1, f2) from the gas law. The rest of the variables can be subsequently determined. For instance, Eq. (32) can be written for T&+1, tz) thus, TG(I+l,

b)=~G(fl+TG(I+l,

(38)

b){TG(~-~G(~)).

The terms T~(I+ 1, t2) and r&I+ 1, t2) in Eq. (28) are eliminated via Eqs. (29) and (38), respectively. The result can be explicitly solved for G(I+ 1, f2). And with G(I+ 1, tz) then at hand the velocity u(l+ 1, t2) can be computed from Eq. (25). Solution of the solid temperature UT(t2, K) is by means of the Thomas algorithm (18). Finite difference approximations of diffusion equation

Using Eq. (35) it is possible to transform the diffusion equation (Eq. (5)) to the form aUT(J+l,K)=

at

(Y

a2UT(J+1,K) ar2



t,=J+l.

(39)

If we define at

A=a(Ar)'

(40)

then all the methods of finite differencing the diffusion equation can be written into a general finite difference form, thus CJTfJ + 1. K) - t 77-f .1. Kj =A6{UT(J+l,K+l)-2UT(J+l,K)+UT(J+l,K-1)) + A(1 - O){UT(J, K + 1) - 2 UT@, K) + UT(J, K - 1)); 292

lourn.

0<8<1. of The

Franklm

(41)

lnst~rute

Heal Transfer in Packed Bed Heat Recuperators-I It can be shown that this expression reduces to the Schmidt equation for the special case of 8 = 0. Also if 8 = 1 is substituted in Eq. (41) the result is that of von Neumann. Subsequent researches produced two notable results. One is the Crank-Nicolson method (6) which is stable for all A for the case of 8 =f. The Douglas method (7) is likewise stable for all A when 0 =f- 1/12A. Aside from being unconditionally stable the method by Douglas has also a high order of accuracy since it has the least truncation error. In the present work the Douglas method was used in the data analysis. The Thomas algorithm The general finite difference expression for the diffusion equation (Eq. 41) can be rearranged to separate all the terms at current time level .I+ 1 from the terms at time level J so that, for convenience -28UT(K+1)+~(1+2bB)UT(K)-28UT(K-l)=DE(K)

(42)

where DE(K)=2(1-B)UT(I,K+l)+;{l-ZA(l-B)}UT(J,K) +2(1-@UT(J,K-1). All the terms on the left-hand side of Eq. (42) are understood level J+ 1. By defining

(43)

to be at the time

BETA = f

(44)

UT(K) = BE(K) * UT(K + 1) + GE(K)

(45)

and and by using the Douglas criterion 1

e=?-=

1

(46)

one can show (I-7) (47) BE(K)=(2+~BETA)-(l-~)BE(K-I) DE(K)+(l-y)GE(K-I) GE(K) =

(48) (Z+zBETA)-(l-y)BE(K-1)’

Before engaging the Thomas algorithm the limiting values for BE(K) and GE(K) at K = 1 must first be determined. In addition to considering the stability and accuracy of finite difference approximations of the diffusion Vol

3U3.

No.

3. March

1977

293

S. H. Jury and M. C. Berbano equation and the bed equations, one must also take into account these same effects when finite differencing the boundary equation for convective heat transfer. This is discussed in the work of Price and Slack (25). Consequently, the finite difference approximation of the boundary equation at the gas-solid interface (i.e. Eq. (6)) can be obtained to produce the result &

UT&, 2) = TG(I+ 1, tJ-

UT(r2, 1)

(49)

where N is the number of radial increments within the solid particle. This expression is substituted in the equation for the bulk gas temperature (Eq. 36) to eliminate TG(I+ 1, tz). The result is solved for UT(t*, 1) thus Nk, UT(2,l)

=

(N- l)h Ar

QB

UT(2 2)+

ks

,

-+l-QA h Ar

(50)

ks

-+1-QA h Ar

Hence after setting K= 1, one obtains

BE(l)=

GE(l)=

Nks (N- 1)h Ar k L+l-QA h Ar

(51)

(52)

k QB S-+l-QA h Ar

where QA and QB are defined at Eq. (36). Thus, one can carry out the computations for the solid temperatures as follows: (1) Assume functionality for film coefficient and calculate initial values and input parameters to solve for BE(l) and GE(l) via Eqs. (51) and (52). (2) Calculate set of values for BE(K), DE(K) and GE(K) using recurrence formulas given in Eqs. (43,47,48); K = 2, . . . , N- 1. (3) Next calculate DE(N) and GE(N). Computed GE(N) corresponds to the solid temperature at first increment from center of granule, UT(2, N), i.e. UT(2, N) = GE(N). (4) With set of BE(K), GE(K) calculated, compute the solid temperatures toward the surface of the solid via Eq. (43). That is, solve for UT(2, K); K=N-l,..., N-13. Determination

of increment sizes in the three-dimensional

gridwork

In this work the numerical analysis was carried out to determine the sizes of At, Ar and AZ which would give sufficiently small error compared to the limiting result for extremely small values of At, Ar and AZ, 294

Journal

of The Franklin

Institute

Heat Transfer in Packed Bed Heat Recuperators-I

The following are the results. At = 0.0005556 hr (2 set) Ar = 0.0001773 ft (0.00213 in.) AZ = 0.003237 ft (0.0388 in.)

maximum error: maximum error: maximum error:

An overall treatment of the error involved in the experimental tional measurement are reported in the original work (24).

1.24% 1.50% 1.17%. and computa-

Computer programming

The numerical computational scheme that was used in determining the parameters via the least-squares method is shown in the computer flow diagrams given in Figs. 2-4. Only the flow charts for the computer program written for the Drierite breakthrough data analysis (to be considered in Part II) are presented here since the concept of solution and most parts of the program are essentially the same for all the samples. Two experimental recuperator breakthrough runs having different flow rates were simultaneously analyzed. The calculated temperature breakthrough profile was obtained via the Thomas algorithm. Most of the computer time was spent in the subroutine “Thomas” in which the solid temperature distribution is calculated. This is a distinguishing feature in that it permits direct calculation of the bulk gas temperature, and not via the trial and error matching, for example, of the heat flux on the fluid side of the interface with the flux on the solid side of the interface. In the subroutine “Search” the interpolation methods employed in determining the parameters were carried to the point that the relative error became less than 1 x 10-4. Three of the parameters that were determined were the V2(1), V2(2) and V2(3) in the functionality for the film coefficient, i.e. Eq. (2). Another was CW, the parameter for the wall coefficient defined in Eq. (3). The fifth and last parameter was the shape factor for each adsorbent sample, I,!J.The author found that with the specifications for At, Ar and AZ increment sizes as presented in the preceding section, the computer time necessary for complete execution of the computer program to determine the parameters was a little more than 2 hr. Debugging work and the initial phase of preliminary estimation of the initial parameters, i.e. preliminary curve-fitting, were done with the departmental PDP-15 computer. References (1) T. E. W. Schumann, “Heat transfer: A liquid flowing through a porous prism,” J. Franklin Inst., Vol. 208, pp. 405-416, 1929. (2) A. Anzelius, %ber Erwarmung Vermittels Durch Stromender Meiden,” Zeit, Angew. Math. Me&., Vol. 6, p. 291, 1926. (3) W. Nusselt, “Die Theorie des Winderhitzers,” Zeit. ver. deut. kg., Vol. 71, p. 85, 1927. (4) C. C. Furnas, “Heat transfer from a gas stream to a bed of broken solids,” Trans. AIChE, Vol. 24, pp. 142-193, 1930.

Vol.

303. No. 3. March

1977

295

S. H. Jury and M. C. Berbano

I

FIG. 2. Flow chart of main computer program.

296

Heat Transfer in Packed Bed Heat Recuperators-I Main Program Input data -- RETA, NN, N; initial solid temp. distribution UT(l,J)=O J;t.N; boundary vatues for BE(K) GE(K) at K=t

r K=K+l

I Compute BE(K). DE(K) and GE(K) using the Thomas recurrence formulas

I

NO

Calculate DE(K) ond GE(K) at K=N Computed GE(N) 1s solid temp. at first increment from center 1

I L’LiI

L=2

1

-4 K-N-L*1

I

Compute UT(2,K) via recurrence formula

I

FIG.

2.

Flow

chart of subroutine

THOMAS.

(5) S. H. Jury and W. Licht, Jr., “Drying of gases-adsorption wave in desiccant beds,” Chem. Engng. Prog., Vol. 48, pp. 102-109, 1952. (6) J. Crank and P. Nicholson, “A practical method for numerical evaluation of partial differential equations of heat conduction type,” Proc. Cambridge Phil. Sot. Math. Phys. Sci., Vol. 43, pp. 50-67, 1957. (7) J. Douglas, Jr., “The solution of the diffusion equation by a high order correct difference equation,” J. Math. Phys., Vol. 35, pp. 145-151, 1956. (8) R. B. Bird, W. E. Stewart and E. N. Lightfoot, “Transport phenomena,” Wiley, New York, 1960. (9) B. W. Gamson, G. Thodos and 0. A. Hougen, “Heat mass and momentum transfer in the flow of gases through granular solids,” AIChE J., Vol. 39, pp. l-35, 1943. (10)W. H. Gauvin and S. Katta, “Momentum transfer through packed beds of various particles in a turbulent flow region,” AIChE J., Vol. 19, pp. 775-783, 1973. (11)B. M. Gillespie, E. D. Crandall and J. J. Carberry, “Local and average interphase transfer coefficients in a randomly packed bed of sphere,” AIChE J., Vol. 14, 3, pp. 483-490, 1968.

Vol.

303, No. 3, March

1977

297

S. H. Jury and M. C. Berbano

FIG. 4. Flow chart of subroutine

298

SEARCH.

Journal of The Franklin Institute

Heat Transfer

in Packed

Bed Heat Recuperators--I

(12)H. Green, “The effect of non-uniformity and particular shape on average particle size,” J. Franklin Inst., Vol. 205, pp. 713-729, 1927. and heat transfer mechanism in (13) D. Handley and P. J. Heggs, “Momentum regular-shaped packings,” Trans. Inst. Chem. Engnr., Vol. 46, pp. T251-T264, 1968. (14) D. Handley and P. J. Heggs, “The effect of thermal conductivity of the packing material on transient heat transfer in a fixed bed,” Znt. .J. Heat & Mass Trans., Vol. 12, pp. 549-570, 1969. of breakthrough curves in packed beds: I(15) C. P. Jefferson, “Prediction Applicability of single parameter models,” AZChE., Vol. 18, 2, pp. 409-420, 1972. (16) S. H. Jury, “Numerical solution of systems of mixed partial differential equations associated with transport phenomena with exchange,” J. Franklin Inst., Vol. 292, 3, pp. 203-223, 1971. (17) S. H. Jury, “An improved version of the rate equation for molecular diffusion in a dispersed phase,” AZChE J., Vol. 16, 1, pp. 1124-1126, 1967. (18) G. H. Bruce, D. W. Peaceman, H. H. Rachford and J. D. Rice, “Calculations of unsteady-state gas flow through porous media,” Trans. Am. Inst. Mining & Metallurgical Engnr., Petroleum Branch Vol. 198, pp. 79-92, 1953. (19) Y. V. Kudryavtsev ef al., “Unsteady State Heat Transfer,” Illiffe, London, 1966. (20) P. K. Leung and D. Quon, “A computer model of the regenerative bed,” Can. J. Chem. Engng., Vol. 44, p. 26, 1966. heat transfer between air and (21) G. 0. G. LSf and R. W. Hawley, “Unsteady-state loose solids,” Znd. Engng. Chem., Vol. 40, (lo), pp. 1061-1070, 1948. (22) D. F. Moline and J. 0. Hougen, “Thermal conductivity of granular solids through which gases are flowing,” Chem. Engng. Prog., Vol. 48, (3), pp. 147-149, 1952. (23) J. R. Mondt, “International Development in Heat Transfer,” Boulder Meeting, p. 614, 1961. (24) M. C. Berbano, “The Kinetics of Heat Recuperation,” Ph.D. Thesis, The University of Tennessee, Knoxville, 1974. (25) P. H. Price and M. R. Slack, “Stability and accuracy of numerical solutions of the heat flow equation,” Br. J. appl. Phys., Vol. 3, pp. 379-384, 1952. (26) 0. A. Saunders and H. S. Ford, “Heat transfer in the flow of gas through a bed of solid particles,” J. Iron & Steel Inst., Vol. 141, (l), pp. 291-316, 1940. (27) J. R. Arthur and J. W. Linnette, “The interchange of heat between a gas stream and solids particles-Part I,” J. Chem. Sot. (Part I), pp. 415-424, 1947. (28) J. M. Smith, “Chemical Engineering Kinetics,” 2nd Ed., McGraw-Hill, New York, 1970.

Vol. 303, No. 3, March 1977

299