these will lead to a lower rate of deactivation. By the time the catalyst at the reactor entrance needs reactivation the catalyst near the exit may still have a long life expectation before reactivation. A better operation is to reverse the direction of flow of the reaction mixture after a certain time period. I n this way, the catalyst in the whole bed will reach the time for reactivation nearly simultaneously, and the catalyst service life before each reactivation \vi11 be prolonged. The effects of h. { B , and on p are relatively low and are not shown here.
rc
ra
= rate of main reaction = rate of adsorption
rd
=
rs
= rate of fouling reaction = actual rate of main reaction = rate of reaction at 0 = 0 if interior of the catalyst
r
r1 rl
SI
= = = =
T t x
rate of desorption
pellet were equally available for reaction adsorbed fouling material S absolute temperature time spatial coordinate
GREEKLETTERS Acknowledgment
The author expresses his great appreciation to UCLA Computing Facility for the use of the computer.
a, CYf = coefficient of 0 in Equations 38 and 37, respectively = DA.~A~/DB~PB, CB = ksDAepao/kDt,pco {c = effectiveness factor tl 170
= initial effectiveness factor
Nomenclature
tlR
=
8
+,
= = = = =
av
= time average. [XI,,, =
A AI B
BI C
= = = = =
= = = = =
=
= = = = = = = =
reactant A adsorbed ,4 product B adsorbed B component C concentration of adsorbed component i adsorbed C effective diffusivity of component i in the internal pores of the catalyst pellet local vacancy factor, the fraction of adsorption sites that remain active at a certain point average vacancy factor gas modified Thiele modulus adsorption equilibrium constant of component i
K,p,,
reaction velocity constant of main reaction adsorption velocity constant of component i desorption velocity constant of component i reaction velocity constant of the fouling reaction total concentration of adsorption sites. active or inactive = active center = partial pressure of component z = partial pressure of component i at external surface = radius of spherical pellet = gas constant
~0.10
E D
relative effectiveness, q / q o dimensionless time, k,t 0 at which q R = 0.10 dimensionless spatial coordinate, x / R relative rate of deactivation dimensionless partial pressure of component i, pi/pio
fJt
x
fJ
0
dt =
x de
Literature Cited
Archibald, R. C., Greensfelder, B. S., Znd. Eng. Chem. 37, 356 (1945). Herington, E. F. G., Rideal, E. K., Proc. Roy. SOC.(London) 184A, 414 1194;) \ - -
’ - / ’
Hogan, J. P., Banks, R . L., Lanning, W. C., Clark, A , , Znd. Eng. Chem. 47. 752 (1955). Masamune; S., Smith, J . M., A.Z.Ch.E. J . 12, 384, 1966. Peters, E. F., Zletz, A . , Evering, B. L., Znd. Eng. Chem. 49, 1879 11957’1.
Pozzi, A’, L., Rase, H. F., Zng. Eng. Chem. 50, 1075 (1958). Voorhies, A , , Znd. Eng. Chem. 37, 318 (1945). Wheeler, A , , Adban. Catd. 3, 250 (1950).
CHIEH C H U Unioersity uf Calqornia Lus Angeles, Calif. RECEIVED for review April 24, 1967 ACCEPTED May 2, 1968
NUMERICAL SOLUTION OF COUPLED, ORDINARY DIFFERENTIAL EQUATIONS A wide variety of problems involving ordinary differential equations can be linearized about a trial solution and then put into finite difference form. The resulting coupled, tridiagonal matrices can b e solved readily on a digital computer, and the nonlinear problem can then be solved by iteration. problems in the physical sciences can be reduced to ordinary differential equations. The availability of high speed digital computers and a generalized method of solution allows such problems to be treated without the drastic approximations frequently needed to obtain analytic solutions. The original problems are often nonlinear and involve several dependent variables, but by a proper linearization of such problems a convergent iteration scheme usually results, although convergence cannot generally be assured. Consequently, it is appropriate to discuss first the solution of coupled, linear differential equations, followed by illustrations of the linearization me tliod. ANY
514
I&EC FUNDAMENTALS
Since other, very different techniques work well with initialvalue problems, attention is restricted here to boundary-value problems-that is, with boundary conditions at x = 0 and at x = L or x = m . Such problems of interest to chemical engineers arise, for example, in the following situations: Mass transfer into a semi-infinite, stagnant medium ((‘penetration” model). Mass transfer in a stagnant film or a porous solid, as encountered rvith heterogeneous catalysis or porous electrodes. Mass transfer in boundary layers possessing profiles similar in the distance along a surface. This can include both free and forced convection, and for large Schmidt numbers the similarity of the hydrodynamics ceases to be esential.
Velocity distributions in similar boundary layers. Distribution of charge and mass in diffuse. electric double layers (perhaps of less interest to chemical engineers). The concentration distributions of several species may be coupled among themselves or with the velocity and temperature fields for a number of reasons. The diffusion coefficients: viscosity, and other physical properties depend upon the composition, temperature, and pressure. The interfacial velocity at the surface is related to the rate of mass transfer. The species may be charged and interact with each other through the electric potential. The components may be involved in heterogeneous or homogeneous reactions described by equilibrium or rate expressions. For free convection, the fluid motion results from density differences created by nonuniform composition or temperature.
j = l
I x =-h
2
3
4
5
I
1
I
I
I \
0
h
2h
3h
L
I
ZIrnoge point Figure 1 . Image point for treatment of boundary conditions with derivatives
If we use the finite difference form of Equation 3, we can expedite matters by the introduction of an “image point” at x = -h, outside the domain of interest (see Figure 1). Then the coefficients in Equation 10 become Z f , k
= - - B t , k ( l ) = ‘/z
Di,k(I) = hQ,k,
disk>
Gi(1) = hfi The references illustrate some problems which have been or could be treated by the present method (Grens and Tobias, 1964; Sewman, 1963, 1966; Tiewman and Hsueh, 1967; Newman and Tobias, 1962). Coupled, Linear, Difference Equations
A set of n coupled, linear, second-order differential equations is represented as
where y k ( x ) are the n unknown functions. Subscript i denotes the equation number, and each of the equations can involve all of the u n k n o w ~ ~ys k, , through the sum. The equations are linear-that is, coefficients ai,k,bi,k:and ct,k are independent of the unknowns. For central difference approximations of the derivatives with a mesh distance h,
d2yk,dx2
=
bi;(xj + h) f
dyJdx =
b k ( ~ f j
~
h)
~ jh ) - 2 y k ( s j ) ] / h 2 + O(hZ) - yi;(xj - h)]/2h f O(h2)
k
(
one obtains the following difference equations n
A l , k ( j ) Y k ( j- 1) k=l
+ B i , k ( j )Y k ( j ) + =
Zi,eYk(j- 2 )
D t , k ( j ) Y t ( jf 1) = G A j )
(4)
Pk(Xj)
(5)
Bi,,(j)
-
= -2ai,k(.Yj)
Df,k(j)=
ai,ic(xj)
(h/2)bt,k(,Yj)
(6)
+
(7)
h2Ci,k(Xj)
+ (h/’2)bi,tjxj)
G t ( j ) = h*gi(xj)
=
Gi(j) (13)
where coefficients Zi,k again allow the introduction of complex boundary conditions at the upper limit (x = L ) of the domain of interest, in the same way that the coefficients xi,kdo at x = 0.
\Ve turn now to the method of solution of the difference Equations4: 10, and 13. F o r j = 1, let Y k ( j )take the form n =
Ai,k(j) = Qt,k(XI)
+ Ai,b(j)Yk(j- 1) f B200 storage locations for the calculation of 20 components and 100 stages. compared to 10,000 locations for the largest array in the earlier method. The calculation time was perceptibly lonqer for the present method. The earlier method takes advantage of the fact that one set of variables, the temperatures on the stages, controls all the compositions if the flow rates are fixed. The subsequent adjustment of the flow rates to satisfy the enthalpy balances is less certain, whereas the present method simultaneously treats the enthalpy balances in the linearization. The present method would make it easier to introduce composition-dependent activity coefficients for the liquid phase without upsetting the convergence characteristics. The present method provides greater flexibility in the specification of column operation, allowing the specification of the bottom-product flow rate and the reflux flow rate to be replaced by conditions such as a given mole fraction of a component in the top or bottom product, a given product amount of a component, a fixed reboiler or condenser temperature? or a fixed reboiler or condenser heat load. Discuss ion and Conclusions
The procedure outlined here for solving coupled, nonlinear, difference equations by linearization and subsequent iteration is general and flexible and has proved useful in a number of problems. For special problems it may be possible to devise more efficient methods, but with a loss of generality and an expense of effort.
Two other methods might occur to one faced with a problem of the type treated here. One is to linearize and decouple the equations by taking the coefficients of the derivatives to be given by a trial solution-for example, approximate yldyz/dx by y1Odyz/dx. Then the decoupled equations are solved one after another in a cyclic process producing new functions to be used as a trial solution. I n general, the convergence behavior is poorer than for the present method, although there are special problems where the coupling is not strong and the method works. A second method would be to treat the problem as an initial value problem and to fabricate the needed initial conditions. This method requires little storage space, but the adjustment of the added initial conditions so as to satisfy the boundary conditions a t x = L can be tricky or impossible. This method is similar to plate-to-plate distillation calculations, which are very sensitive to, say: the value used for the mole fraction of a light component in the bottom product. The errors in the present method arise from three sources: convergence errors for the nonlinear problem (which can be made negligibly small here), errors in the difference approximation to the differential equations (which decrease with the mesh distance, h ) , and round-off errors in the computer (which increase as the mesh distance is decreased). Convergence may not be possible if there are sharp variations of the unknowns in some region of x; in such a case a singular-perturbation method may be appropriate. A remark is in order for first- and third-order differential equations. For the purpose of computation, a third-order equation involving, say, d 3 ~ t / d x 3could be replaced by two second-order equations (with w = dyJdx and y2 = dre/dx) or a first- and second-order equation (with w = y l and j z = dai/dx). I n this way the finite difference forms still involve only the points at j - 1, j , and j 1. For a first-order equation it is
+
probably better to use a backward difference rather than a central difference. The order of the approximation will still be h2 if the coefficient takes on its average value: K d j i l d x = '/z
[K(j)
+ K ( j - l)][Yi(j) Y l ( j - l)l/h +O(h2)
(25)
Availability of Computer Programs
A general program for solving coupled, linear, difference equations as outlined here has been developed and is available from the author. This can be used as a subroutine for programs which set up the coefficients A i , k , Bi,k, etc., and control the iteration procedure. This report (UCRL-17739) also contains two distillation programs, the earlier one (Newman, 1963) with improvements in the calculation of compositions and the adjustment of flow rates to satisfy the enthalpy balances and one which uses the more general linearization procedure. literature Cited
Grens, E. A , , 11, Tobias, C. W., Ber. Bunsenges. Physik. Chem. 68, 236-49 (1964). Newman, John, Hydrocarbon Process. Petrol. Refiner 42, No. 4, 141-4 (196'3). \----,.
Newman, John, IND. ENG. CHEM.FUNDAMENTALS 5 . 525-9 (1966). Newman, John, Hsueh, Limin, Electrochim. Acta 12, 417-27 (1967). Newman, John, Tobias, C. W., J . Electrochem. Sod. 109, 1183-91 (1962).
JOHN XEWMAN Uniuersitj of California Berkeley, Calif. 94720
RECEIVED for review August 21, 1967 ACCEPTEDMarch 4, 1968 Work supported by the U. S. Atomic Energy Commission,
FILM BOILING FROM A SPHERE An analysis of film boiling from a sphere to a saturated liquid in forced convection is presented. The heat transfer coefficient from a sphere is similar to that obtained by Bromley for forced convection film boiling from cylindrical tubes.
c
onsider a sphere in a floxving saturated liquid with film boiling occurring a t the sphere surface. Kobayashi (1965) performed an analysis of film boiling for this case in a manner similar to that of Bromley et al. (1953) for a cylinder. Kobayashi's analysis was refined by Hesson and Witte (1966). The analysis of Kobayashi resulted in a set of equations which were amenable to a numerical solution only. This paper presents an analysis which results in a simple expression for the heat transfer coefficient from a sphere to a flowing, saturated liquid.
-c
Figure 1. Model for film boiling around a sphere
Figure 1 shows the model representing the film boiling phenomenon. A solution for the heat transfer coefficient for this case can be achieved under the following assumptions. The
vapor and liquid are both pure and incompressible. The sphere surface temperature is uniform. Thermophysical properties are constant with temperature. The flow is laminar. T h e vapor layer is thin compared with the sphere radius and the liquid-vapor interface is smooth. The velocity and temperature profiles in the film are linear. Radiative heat transfer may be neglected. The liquid velocity a t the interface can be calculated from potential flow theory. An energy-mass balance a t the liquid-vapor interface for a small differential element of vapor adjacent to the sphere surface gives kATdA ~6 - hf,'dw
(11
where k is the thermal conductivity, AT is the sphere wall temperature minus the saturation temperature, dw is the rate of vapor produced in the element, and hfo' is the latent heat of vaporization as derived by Rohsenow (1956) hiv' = h,, VOL.
+ 0.68 C,AT
7 NO. 3
AUGUST 1 9 6 8
(2) 517
