NONLINEAR EQUILIBRIUM AND AXIAL MIXING EFFECTS I N INTRAPARTICLE DIFFUSION-CONTROLLED SORPTION B Y ION EXCHANGE RESIN BEDS Computer Analysis CHARLES J . COLWELLI A N D J O S H U A S. D R A N O F F Department of Chemical Engineering, Sorthwestern University,Euanslon, Ill. 60101 The effects of nonlinear equilibrium and axial mixing on sorption processes in packed ion exchange resin beds have been studied by digital computer simulation of column performance. Results show clearly that the presence of these effects leads to false estimates of apparent intraparticle diffusivities when a linear model is used in analyzing data. A technique for estimation of true diffusivities from apparent values is presented.
S O R P T I O N of nonionic solutes by packed beds of ion exchange resins has been of interest for many years, both as a potential separation technique and as a model for other liquid-phase sorption processes. Recent work in this laboratory (Colwell, 1965; Colwell and Dranoff, 1966) with dilute aqueous solut'ioiis and packed beds of Dowex-50 (H+) resin demonstrated that intrapart,icle diffusion was the rate-limiting step in such processes. 111addition, it was shown t'hat experimental breakthrough curve data could be represented excellently by the linear equilibrium spherical part,icle intraparticle diffusion model of Rosen (1952, 1954) after proper choice of the diffusivity. However, the resultant diffusivities were different for colunin saturation and elution experiments and apparently affected by flow rates in the packed bed. I t was suggested that slight iionlinearities in the equilibria and axial mixing due to density gradients in the liquid streams were responaihle, respectively, for these variations in apparent diffusivities. The goal of the work described here was to determine the validity of these proposals theoretically by means of digital computer siniulation of appropriate packed-bed models. The computey techniques used and the results obtained are detailed below. This study verified the significance of nonlinear equilihria and axial mixing on apparent intraparticle diffusivitiea anti resulted in a method for estimating true values froiii the apparent ones under certain conditions. The general approach and conclusions discussed are not restricted to specific systems. Their application to some real experimental data will he considered in a subsequent paper.
causes, may be characterized by a diffusional model, the basic differential material balance for solute may be written as in Equation 1,
which relates the bulk or average concentrations in the liquid and solid phases. The latter is given in terms of the local concentration within the resin bead by Equation 2.
The variation of local solid-phase concentration is described by the usual transient material balance for a sphere under the assumption that a diffusion coefficient D will characterize the rate-limiting intraparticle diffusion mechanism.
Appropriate boundary and initial coridit ions for these balances are indicated in Equations 4 to 9. c(2,
0) = c i
(5) c(c0, ' J
Mathematical Model
The process to be simulated is the transient response of a packed bed of uniform spherical particles to a step change in liquid feed composition. The concentrations of solute in the liquid and solid phases are initially constant over the bed and in equilibrium. Liquid flows through the column with a constant velocity which is uniform over the cross section and operation is assumed to be isothermal. With the assuniption that axial mixing, due to whatever Present address, Esso Research and Engineering Co., Florham Park, N. J.
(4
t ) = cz
(6)
O ) = pi
(7)
q (0, z , t ) is finite
(8
p ( E , 2 J t , = K D c (2, t )
(9
Equation 5 is the usual flux boundary condition for the inlet, while Equation 6 corresponds to an open column (Levenspiel and Smith, 1957; Wehner and Rilhelm, 1956) at the downstream end. The last boundary condition expresses the assumption of equilibrium at the particle-solution interface. The choice of a distribution coefficient ( K O ) form of equilibrium relation is made for convenience in linear systems for which the coefficient is constant rather than a function of concentration. Finally, the initial concentrations, ct and pi, also satisfy this equilibrium relation. VOL.
8
NO.
2
MAY
1969
193
1.0
ditions by means of Duhamel’s theorem and other suitable manipulations. This technique has been used by others (Antonson, 1968; Rosen, 1952, 1954; Tien and Thodos, 1959) and is presented in detail elsewhere (Colwell, 1967). Negligible Axial Diffusion. In the case of negligible axial diffusion these manipulations lead to the following equation :
0.8
06
6 0.4
0.2
exp [-n2t,b (r - r ’ ) ] dr’= 0
9 Figure 1.
(20)
where
Nonlinear equilibrium isotherm
This system of equations may be rewritten in a dimensionless form equally applicable to saturation and elution experiments, as follows:
ap - + ap - + fae _ _ _ _ a2p -=o ax ar e a x f f L8x2
(10)
and
Further simplification follows by the change of variables : ?=r-x
(23 1
Z=X
(24 )
with the result that the system equation is now
%+@I= 0 ax Here I is a convenient shorthand notation for the infinite series-integral term. (14) cp(“3,T)
=
S ( P , X,O) =
0
(15)
0
(16)
6 (0, x, r ) is finite
eU,X,
7 ) = K D (x, ~ r
(17)
)
(18)
The type of equilibrium relation considered in this work is indicated in Equation 19, which expresses the distribution coefficient as a linear function of solute concentration.
+
KD = K D ~ [ ~G(1 - cp)] = Ko0f(p)
l&EC
FUNDAMENTALS
p ( f , i I 0) = 0
(28)
The technique of solution is described in the Appendix. Significant Axial Dispersion. In the case of significant axial dispersion, the reduction of the basic equation leads to
(19)
KD, represents the value of the distribution coefficient when the solid is saturated with respect to the feed solution concentration, q = 1. Choice of the parameter G determines the “amount of nonlinearity” in the equilibrium isotherm. Typical behavior of this type of function is illustrated in Figure 1 for KO, = 1.0 and various values of G. G > 0 corresponds to the usual favorable equilibrium case, while G < 0 represents an unfavorable isotherm. This particular form was selected, since it is a simple relation that is representative of the type of data of concern of this work. Equations 10 to 19 describe mathematically the general problem of interest in this work. If a~ and G are both zero, the equations represent the linear model of Rosen (1952, 1954)) while with CYL alone equal to zero, the model corresponds to nonlinear equilibrium with negligible axial dispersion. Computer solutions were obtained for each of these cases on a CDC 3400 digital computer, using finite difference programs devised for this study. The cases of significant and negligible axial dispersion were handled separately, since somewhat different approaches were utilized. For each problem the equations were reduced to a single integro-differential equation with suitable boundary con194
The equation for p(Z, i ) must be solved with the boundary conditions: (27) cp (0, i )= 1
The definition of I is again as in Equation 26, except that r replaces i . I n that connection, there is no point in transforming variables as done previously, since simplification of the equation does not result because of the presence of the second derivative term. The boundary conditions remaining to be satisfied with Equation 29 are:
CPb,0 ) = 0
(30 1
( P ( w , r )= 0
(32 1
The numerical scheme used in solving this problem is also described in the Appendix. Analysis and Discussion of Computed Results
The approach taken in analyzing the results of computer calculations with the above-mentioned programs was to treat the data in the way experimental data might be handled.
Breakthrough curves-that is, curves of cp us. r for fixed values of z-were generated using numerical values for the several physical parameters corresponding to the ranges encountered in the related experimental work. The magnitudes of these parameters are: E
= 0.389
f = 0.395
L= A = u= D=
90.8 em. 5.07 sq. cm. 0.0833 to 1 cc./sec. 0.7 to 2 X
sq. cm./sec.
DL = 0.002 to 0.045 sq. cm./sec. R = 0.01 to 0.05 em. KD, = 0.5 to 2.0 The computed breakthrough curves were then analyzed by means of Rosen's model through the process of matching computed and model curves. The matching was done visually, as with experimental results, and led to a characteristic value of the parameter X (Rosen, 1952, 1954) for each curve. Since
l o w range negligible axial dispersion
this value determines the apparent intraparticle diffusivity, De,,, for that calculated curve, The results are reported below not in terms of D,,,, but rather for more generality as the ratio of the apparent value of to its known value for the curve in question. Thus, this ratio
is a direct measure of the error which arises in the estimation of the intraparticle diffusivity using the linear model when the true equilibrium is nonlinear, all other parameters held constant. Over the range of variables considered, it was possible to fit the computed breakthrough curves well with those of the Rosen model. However, such behavior is found only in a region where the nonlinearities are not too severe. As a guide in correlating the results of these comparisons, the form of the general dimensionless equation and boundary conditions (Equations 29 to 32) was considered. This indicated that the breakthrough curve and hence the values of paPpand ultimately B will be unique functions of four parameters, assuming the column and particle void fractions to be constant. Thus, it is expected that (35 ) Nonlinear Equilibrium with Negligible Axial Dispersion. With D L = 0 and G and KO, fixed, B should depend on (Dt,/R2). Results obtained for various values of the two equilibrium parameters are thus shown in Figures 2 and 3. Figure 2 corresponds t o (Dto/R2)less than 0.9, representing small columns and/or high flow rates. The curves drawn are for KO, = 1.27 and several values of G. Also shown are some individual data points for other values of KD,. Figure 3 presents results only for KO, = 1.27 and G = f 0 . 2 over a larger range of (Dto/R2). More data were not obtained in
OWR' Figure 3.
Diffusivity ratio vs. Dfo/R2
High range negligible axial dispersion
this region, since it corresponds physically to a regime in which axial dispersion becomes significant, thus invalidating the basic model in any case. These two figures point up several very significant facts. First of all, the existence of nonlinear equilibrium causes the apparent intraparticle diffusivity to differ from the true or actual value. Furthermore, when equilibrium is unfavorable (convex upward or G > 0), B is greater than 1.0, making the apparent diffusivity higher than the true value. If equilibrium is unfavorable (convex downward or G < 0), the converse is true. The effect is much greater for the favorable equilibrium case. Secondly, the increases in (Dto/R2)make the deviation of B from 1.0 greater in every case. Thus, as column length is increased or flow rate is decreased, the apparent diffusivity rises if G > 0 or falls if G < 0. This leads to some of the apparent flow rate effects which have been observed previously and falsely ascribed to the phenomenon of film transfer resistance. Several other points are also evident from these data. Figure 2 shows that the curves have significant flat portions for G 5 0.2 within which the apparent diffusivity remains constant. In addition, the value of B is relatively insensitive to KO, when (Dt,/R) < 0.2. The results of Figure 3 for VOL.
8
NO.
2
MAY
1969
195
G = +0.2 correspond to the well known constant pattern breakthrough region. It should be now clear from these data that serious errors may result in the estimation of D in this region if nonlinear equilibria exist. The data of Figures 2 and 3 are restricted to values of G between $0.4 and -0.4. Outside this range it was difficult to fit the computed breakthrough curves well, using the theoretical curves of Rosen. This then defines the zone (for the assumed form of equilibrium isotherm) in which the linear model may be applied in the estimation of intraparticle diffusivities. Within this region, curves such as Figure 2 or 3 may be used to correct apparent diffusivities, as long as axial dispersion is not important. Outside this region, data may be analyzed directly in terms of calculated breakthrough curves for the particular equilibrium isotherm of interest. Generalization of such curves may be impossible, however, because of the many different types of equilibrium models which may be considered. Intraparticle Diffusion with Axial Dispersion. The effects of axial dispersion on the estimation of intraparticle diffusivities were considered next. The principal difficulty encountered in this work was the long computer time required to generate significant data. As a result, fewer data points were obtained than in the previous case. However, sufficient computations were done to make clear some of the general trends to be expected when axial dispersion is important. The measure of the effect of axial dispersion on the estimation of the intraparticle diffusivity has been defined so as to separate the nonlinear equilibrium effects. The ratio used in this case is defined as B* =
Paw
(36)
(Pam )D ~ - 0
If the equilibrium is linear, G = 0 and B* = B. Otherwise, the diffusivity ratio may be found from:
(37) The nature of the results obtained with the computer simulation of this model is illustrated in Figures 4 and 5. Figure 4 presents data for the case of G = 0-that is, linear equilibrium. As indicated in Equation 35, if all other parameters are constant, B* should depend uniquely on DL. Figure 4 shows results for three values of volumetric flow rate and physically likely values of DL. The effect of axial dispersion is most significant a t low flows, becoming almost
0*
0 Figure 4.
0.2
0.4 Q6 0.8 1.0 O, x IO2, cm*/sec.
Effect of axial dispersion on B* linear equilibrium
196
I&EC
FUNDAMENTALS
0*
o
0.2
Figure 5.
0.4
X
0.6
ae
0.1
Effect of column length on B*
KD, = 1.27 D = 1.1 X 10-6 rq. cm./rec.
negligible as the flow rate is increased. Data for G values other than zero show similar behavior but with different curves in every case. Figure 5, a composite of several plots, shows clearly the importance of column length on the estimated diffusivities and demonstrates that the length effect is more significant for favorable than for nonlinear equilibrium. For G = 0 and G = $0.2 the curves are apparently peaked, while for G = -0.2 this may or may not be so. Higher values of DL were required for the latter case in order to produce effects somewhat comparable to those for favorable or linear equilibrium. This difference reflects interplay of axial dispersion and the equilibrium. In the favorable case the concentration front moving through a column tends naturally to become sharper, leading to steeper concentration gradients and in turn enhancing axial dispersion. On the other hand, unfavorable equilibria tend naturally to spread the breakthrough curve, leading to diffuse fronts and consequently less significant dispersion effects. Similarly, the apparently large effect a t small column length in every case is due to the step change in concentration a t the column inlet and the resulting steep concentration gradients near that point. The axial dispersion mechanism tends to propagate these changes down the column rather rapidly a t first, thus yielding a false (low) impression of the intraparticle diffusivity. As the initial step is diminished by absorption in the solid, dispersion may become less significant. The net result is that breakthrough curves from a system with axial dispersion are slightly above the Rosen curves a t the beginning but fit well thereafter. All of these results show clearly that the effect of any axial mixing is always to yield a low estimate of the true intraparticle diffusivity, the error being most significant in cases of favorable equilibrium. Furthermore, the importance of axial dispersion increases as the flow rate of liquid in the packed bed decreases. It has not been possible with the
small amount of data generated here to uncover a general technique for estimating true diffusivities when axial dispersion is important. Nevertheless, if the axial diffusion coefficient DL were known, general breakthrough curves could be determined for specific nonlinear equilibria. However, the details of such curves will not be apparent without considerable computation, which is best left to the specific case under study. These computer studies have demonstrated without question the significance of nonlinear equilibrium and axial dispersion in the estimation of intraparticle diffusivities by curve-fitting experimental breakthrough data with a linear model. For cases in which axial dispersion is not important, it has been possible to generalize these results over a fixed range of nonlinearity of the isotherm, and so to produce a method by which apparent diffusivities may be corrected to yield better estimates of true values. Where axial diffusion is also important, these studies have demonstrated the nature of the effects which may arise, showing that favorable equilibrium cases are affected more significantly and that column length has a definite effect on the estimated values. I n the view of these facts, it is imperative that future work aimed a t estimating diffusivities by the fixed-bed technique make use of such computer analyses in treatment of experimental data. Application of these ideas will be illustrated in a subsequent paper.
mation order, and replacement of the r derivative as shown in Equation A-4, the expression for I is given by A-5.
(A-4 1
Cexp
{-nV(f -
(K
exp
+ 1)l))-
{-nV(t - K l ) } ]
(A-5)
In the present case, [f(cp)cp] is found from Equation 19 as:
Cf(cp)cpl=
(A-6 1
C(1+ G ( 1 - cp)1cpl
I n general, however, any suitable form may be used for f(cp). The actual evaluation of I ( Z , i ) by direct application of Equation A-5 was not desirable, since the number of terms in the first summation increases with each step in the 7 direction, requiring excessively large computer storage. Instead, a convenient recursive relation was deduced after study of the way in which Z evolves. This relation permits the evaluation of I ( 2 , i I ) in terms of I ( Z , 7 ) . Thus,
+
N
I@,?)
=
c Sj(Z,?)
(-4-7 )
j-1
and
Appendix
Si(*, i )= $ j ( Z , i The numerical methods used to generate the foregoing results are sketched here as a guide to others who may wish to apply them in other situations. The case of negligible axial dispersion is considered first. The basic equation is Equation 25, repeated here.
9az+ P I = O An explicit finite difference method was used for this case. The derivative was represented by a standard-first-order approximation :
The heart of the method lies in the evaluation of I . The integral was first broken into a finite number of integrals as shown in Equation A-2.
- 1 ) exp (-jVZ) +
with Si($, 0 ) = 0 j = 1, N
(A-9 )
Here N is chosen suitably large for convergence. The validity of these relations may be best demonstrated by expansion of the series for several i increments (Colwell, 1967). The storage problem is not entirely eliminated by this approach, since N storage locations must be provided a t each Z increment. However, this number is much smaller than the number of i increments required in the evaluation of Equation A-5 (typically 500 to 1000). A further improvement may be made by noting that the term [exp (-j?Z)] will become small as j increases, leading to a modification of Equation A-7 as follows:
(A-10) exp {
-nv (? - 7’))d d
(A-2)
Then the law of the mean was applied to the integral, with the result:
[y’exp
{ - n v (f - 7’)) d7’
(A-3)
I n this form the remaining integral may be evaluated directly. Following this evaluation, interchange of sum-
Here again J is selected suitably large. Equations A-1, A-10, A-8, A-9, and A-4 thus constitute the framework of the numerical technique. The actual computer program wm designed to calculate the values of cp throughout the upper quadrant of the t - Z plane, moving successively along lines of constant Z for sufficiently large i values. The program contained many specific features designed to speed computation and maintain accuracy, all of which are described in detail elsewhere (Colwell, 1967). The program was tested by comparison with the linear equilibrium model results of Rosen (1952, 1954) and by variation of interval sizes. Computation times were of the order of 2 minutes (for Z = 0.2) on the CDC 3400. No stability problems were encountered in these calculations. The structure of the calculations for the case of significant VOL.
8
NO. 2
MAY
1969
197
axial dispersion will be considered next. The equation to be solved is Equation 29.
Determination of a suitable finite difference scheme for this equation and its related boundary conditions required considerable experimentation. It soon became apparent that an implicit scheme would be necessary if stability problems were to be avoided and reasonable intervals employed. I n addition, however, it was also necessary to develop a particular “diagonal” difference approximation to the derivative terms in order to obtain answers which agreed with the limiting case of no axial dispersion. The usually employed difference representations caused the concentration fronts to propagate down the bed faster than the fluid velocity when cy1 was set equal to zero. The basic scheme finally selected, which avoided this problem, resulted from the approximations shown below :
- + -a(a = ax ar
a(a
$(a(z,r)- 2 y ( z - h, 7 - I ) 1
+
i(a(2
- 2h, 7 - 21)
_ -- (a (z + h, 7 )- 2(a ( z , r )+
a2(a
(a (z
h2
ar2
(A-1 1 )
- h, 7 ) (A-12)
The resulting equation is then:
Nomenclature
A = column cross-section area, sq. cm. B = D,,,/D or /3Rpp//3, dimensionless B* = &p/ (Papp)DL=O, dimensionless c = solute concentration in solution, g./ml. ci = initial solute concentration, g./ml. co = feed stream solute concentration, g./ml. D = true intraparticle diffusion coefficient, sq. cm./sec. D,,, = apparent intraparticle diffusioncoefficient, sq. cm./sec. D L = axial diffusivity, sq. cm./sec. f = intraparticle void fraction, dimensionless f ( p ) = equilibrium function, see Equation 19 G = nonlinear equilibrium parameter, dimensionless h = z-increment for computer program, dimensionless I = infinite series integral, see Equation 26 J = number of stored terms in series for I , dimensionless KD = equilibrium distribution coefficient, dimensionless, see Equation 9 KD, = value of KD when cp = 1, dimensionless, see Equation 19
1
L N
= r-increment for computer program, dimensionless = column length, cm. = number of terms in series for I
q
= solute Concentration in particles, g./mL
qi
= initial solute concentration, g./mL
Q
= average solute concentration, g./ml.
r
= radial coordinate in solid particle, cm.
R SJ
= particle radius, cm. = term in series for I , see Equation A-7
1
= time coordinate, sec.
to
= column holding time, =
u
X
= volumetric flow rate, cc./sec. = dimensionless column length coordinate, = z / L = transformed 2 coordinate, 2 = z = dimensionless parameter of Rosen model, see Equation
z
= axial position coordinate, em.
2
f
LAe/u, sec.
33
[b]
(a (1:
- 2h, r - 21)
+ /3Z= 0
GREEKLETTERS(all dimensionless) (A-13)
This five-point formula required modification for use in the first r row and in the last z column because of the limited data available a t those places. I n addition, the boundary condition at z = 0 was handled with a specific equation for y (1, K ) . The development of these approximations is discussed in detail by Colwell (1967). The implicit scheme embodied in Equation A-13 involves the solution of a set of linear equations for (a as a function of z at each time step. Equations were solved by use of Thomas’ method (Lapidus, 1962). The scheme is an iterative one, since the y values at each step must be assumed a t first in order to permit a n initial estimate of the I term. However, the number of iterations required at each step was reduced to 1 or 2 after the first few rows by means of extrapolations based on Lagrange’s formula (Kunz, 1957) to estimate functions in any new row. The I term was evaluated exactly as described for the previous program. The methods to reduce memory requirements discussed previously were very critical in this case, since it was necessary to carry information for all z equations in each of several rows, including the SJ terms at each z-increment necessary for evaluation of 1. This and the time required for solving the linear equations restrict the use of the program somewhat. The program was tested by comparison with known solutions for the limiting cases of cy^ = 0 and cy, = 0. Computation times were approximately 5 minutes for z = 0.2 and 16 minutes for z = 0.8. 198
l&EC
FUNDAMENTALS
D L A (uL)
(YL
= (axial Peclet number)-’ =
cy,
= (particle Peclet nurnber)-l = Dto/R2 = model parameter, see Equation 22
/3
Pap, A E
6 0
= apparent value of
/3
= see Equation A-4 = interparticle void fraction = dimensionless solid concentration of solute = (q - q d / ( c o - C i )
= average value of 0
T/R
p
= radial coordinate =
r
= time coordinate = t / t ,
i (a
9
= transformed coordinate, T - z = solute concentration in solutions, = ( c = Dto+/R2
ci)/(c,
- ci)
literature Cited
Antonson, C. R., Ph.D. dissertation, Northwestern University, Evanston, Ill., 1968. Colwell, C. J., M.S. thesis, Northwestern University, Evanston, Ill., 1965. Colwell, C. J. Ph.D. dissertation, Northwestern University, Evanston, Id., 1967. Colwell, C. J., Dranoff, J. S., A.I.Ch.E. J . 12, 304 (1966). Kunz, K. S., “Numerical Analysis,” McGraw-Hill, New York, 1957. Lapidus, L., “Digital Computation for Chemical Engineers,” McGraw-Hill, New York, 1962. Levenspiel, O., Smith, W. K., Chem. Eng. Sci. 6, 227 (1957). Rosen, J. B., 2nd. Eng. Chem. 46, 1590 (1954). Rosen, J. B., J . Chem. Phys. 20, 337 (1952). Tien, C., Thodos, G., A.I.Ch.E. J . 6, 373 (1959). Wehner, J. F., Wilhelm, R. H., Chem. Eng. Sci. 6, 89 (1956). RECEIVED for review December 12, 1968 ACCEPTED March 11, 1969
