INVARIANT IIMBEDDING AND COUNTERCURRENT MULTISTAGE OPERATIONS DAVID M. KOENIG 898 Oxford St., Worthington, Ohio 43086
The concept of invariant imbedding is applied to staged countercurrent separation processes and finite difference equlations solvable for outlet quantities as a function of the total number of stages and some properly chosrtn input parameter are derived. A straightforward computational procedure is presented and illustrated in two examples. The first, gas absorption with constant flow rates and linear equilibrium relation, yielded to both the analytical and numerical approach, allowing verification of the derived equations and the numerical results. The numerical solution of the second, a popular solid-liquid extraction problem, allowed comparison with the conventional graphical approach.
H E procedure of converting certain types of two-point Tboundary value problenis into initial value problems by invariant imbedding (Bellman et al., 1960) has been applied to nuclear reactor analysis (Itring, 1962) nonlinear filtering (Lee, 1968), and continuous countercurrent transfer operations (Koenig, 1967). In each application the s p t e m of interest is conventionally described by ordinary differential equations. However, staged systems are conventionally described by finite difference equations to which the concepti of invariant imbedding have not pet been applied. This article shows how the finite difference equations awociated with staged systems can be reformulated by invariant imbedding and how the resulting equations can be solved. Process design problem often fall into two categories. Sometimes the concentrations and flow rates of the incoming streams plus the desired outlet concentrations are specified and the number of equilibrium stages necessary to effect the separation is determined. We call this the “common” design problem. This problem is analogous to a n initial value problem and is dealt with in detail by most unit operations texts. On the other hand, sometimes the number of stages is given along with the incoming stream conditions and the outlet values are to be found. We call this the “existing process,’ problem. This problem is analogous to a two-point boundary value problem and usually requires trial and error solution procedures. We show that both problems can be solved with equal ease by invariant imbedding. Constant Flow Rate Syrtem
Occasionally, in operations such as gas absorption, the flow rate of the gas stream, V , and that of the liquid stream, L , can be taken as constants (Foust et al., 1960). Then a mass balance over the kth stage on the transferrable component gives v y k f Lxk = v y k + l + Lxk-1 (1) Since the streams leaving the kth stage are in equilibrium, we have yk = E ( % ) (2 ) In the common design problem we would try to find the number of stages, N , necessary to produce a specified outlet value, XN or y1, given inlet values L, 20, V , y ~ + 1 . In the existing process problem the roles of X.Y and N would be reversed. For the special case nhere the equilibrium relation is linear,
we can write Equation 2 as yk = Kxk
Equations 1 and 3 can be combined to form a second-order difference equation :
For the common design problem values would be prescribed for xo and either y1 or XN. For the existing process problem values for xo and y ~ + 1would be prescribed. The solution to the latter problem where xo = 0 and y ~ + 1= b is
x/;= ( b / K ) (1 - A k ) / ( l- AN+’)
(5 )
N-here A is the absorption factor, L/KV. Imbedding Formulation
In the invariant imbedding approach we are not interested in internal values of concentration but only in the outlet value. Therefore let us define T N ( ~ . v + ~ )as the value of ZN for a specified value of y ~ + 1 . To derive an equation solvable for TN we set k = N in Equation 1:
VYN4-LTN( Y N + ~ )= VYN+I4- LX.v-1
(6)
Now, if XN is r x ( y ~ + l )then , X N - ~ must be rN-l(yN)-i.e, if TN ( y ~ + is~ the ) outlet value for a specified y ~ + 1in an N-stage system, then r ~ - l ( y ~is) the outlet value for a specified y~ in a N - 1 stage system which can be considered as imbedded in the original N-stage system. Thus, combining Equations 6 and 2 with the new interpretation of X N and letting b = Y N + ~ we get
VE[TN(b)]+ LTN(b) = Vb+ LrN-l(E[TN(b)ll
(7)
subject to the condition T O @ ) = 20. This initial condition simply says that the outlet value of a zero-stage system i s the same as the inlet value. This then is the invariant imbedding problem solvable for r ~ ( bas) a function of N and b. Equation 7 is a finite difference equation of order l subject to one condition, whereas Equation 4 is second-order and is subject to conditions a t two points. NOW,for the special case of a linear equilibrium relation, Equation 7 becomes
VKTN( b )
+ LTN( b )
VOL.
8
=
NO.
Vb f LTN-~[KTN (b)]
(8)
3
537
AUGUST
1969
If we assume separability-i.e., if w‘ ti. r N ( b ) = t N b Equation 8 becomeq, after canceling b and dividing by L, [( V K I L )
+13s
=
(V/L)
+ Ktn.-lf,y
(l/K) (1 - AX)/(I - A”+’)
which agree5 with Equation 5 the validity of Equation 7.
if
k = N , thus demonstrating
In cases where aqiarability not be applied, we must solve Equation 7 as it stands. The main difficulty is due to the second term on the right-hand side, which requires knowledge of rx-1 for a value in equilibrium with r.,-(b). This problem is overconie by the fc?owing procedure : 1. Let AT = 1; then Equation 7, applying the initial condition, becomes
+ Lr, ( b ) = V b + Lxo
(11) 2. Solre for T I @ ) in Equation 11 for several choices of b, sag’ bl, bs, . . .,b , , in the range of interest. The trial and error numerical rnet,hod used for solving Equation 11, once b has been specified, is discussed below. 3. Pass a least squares polynomial, P M ( ~ of ) , order M through t’he values 9-1 ( b l ), r1(b2“, . . ., TI ( b m ) . Numerical experiments showed that sufficient accuracy was obbained with M = 1 and rn = 3 or 4. 4. Let AT = 2; then Equation 7 becomes
+ ~ r . ? ( b ) vb + ~ T ~ p 3 [ r ~ ( b ) l(12) j
V E [ T ~ ( ~ ) ]
=
5. Solre for rz ( b ) for several choices of b, using P2ql evaluated at E[rz(b)]iri place of T ~ ( E [ T z ( in ~ )Equation ]] 12. 6. Repeat step 3 with A‘ = 2. 7. Repeat step 4 wit’h AT = 3.
8. Continue this cycle for larger N .
To solve Equation 11 or 12 we used a simple trial and error routine. An initial guess for ~ , \ r ( b ) , say T’, and a step length s are chosen. Next, [ d ( r ’ ) I 2 is calculated where d ( r ) = I.’E(r)
+ Lr - ‘1’2, - LP.&f[E(T)]
-4 trial step is made and [ d ( r ” ) y is calculated where T” = s. The values of [d(r”jlZ and [d(r’)I2 are compared
T’
+
and a decision is made as to the modification of S. If [ d ( r ” ) ] ’ is greater than or equal to [ d ( r ’ ) I 2 , the step is considered uiisuccessful and s is decreased in magnitude and reversed in sign. If this is not true, the step is considered successful and s can be increased if desired. This procedure continues until [cl (?”)I2becomes “sufficiently” small (less than 0.0001 in our numerical experiments). Best results were obtained when s was divided by - 3 after a n unsuccessful trial and doubled otherwise. The results of a typical computation with A = 1.23 are presented in Table I. Using a time-shared GE 235 computer programmed in BASIC, approximately 3 seconds was consumed in the compilation and execution. The agreement between the calculated values and the true values (using Equation 5 ) demonstrates the validity of the coniputat,ional procedure. Solid-Liquid Extraction
this problem we have three components: the solute (No. l ) , the solvent (No. 2 ) , and the inerts (No. 3). By 111
538
I&EC
FUNDAMENTALS
X N
or
~N(YTV+I)
S o . of Stages in System. A’
Inlet Value, 1/Ntl
Calcd.
True
1
2 4
0.9974 1.989
0.9945 1.989
2
2 4
1.319 2.640
1.321 3.642
3
2 4
1,485 2.955
1.481 2.962
(10)
Computational Procedure
VECT1( b ) ]
Gas Absorption Results Outlet Value,
(9 ) vhich is a Riccati difference equation (Mickley et al., 1957), subject to the condition t o = x&. If we solve Equation 9 for the case where xo = 0, we obtain t,?. =
Table 1.
passing the V stream, originally rich in component 2, counter to the L stream, originally rich in a mixture of components 1 and 3, we wish to reduce the concentration of component 1 in t h L stream. I n the common design problem we would be giveu the L stream inlet flow rate, Lo, and concentrations, x10,X ~ Othe , V stream inlet flow rate, V N + and ~ , the concentrations, y 1 ~ + 1 ,yZh’+l, and a specified outlet concentration, X ~ N . We would then be asked for the number c: stages, N , necessary to effect the separation. The roles of N and r l would ~ be reversed in the existing process problem. The system is described by equations resulting from a n over-all mass balance over the kth stage
+
Lk-1
Vk+l=
Lk
+ Vk
tv o-component mass balances
Ytk+iVk+i= &kLk
Zllt-iLx-1-k
+
(13)
Yrkvk
(1% b )
where i = I, 2, a n equilibrium relation yLk = E ( % )
plus two equation5 relating compositions in the under- and overflowing streams YZk
=
0(Ylk)
x2k =
u(%)
These six equations apply to each of the N stages. I n addition, for each of the 2 N 2 streams we specify hat the mass fractions must add up to unity, so in total the .’T-stage ( 2 N 2 ) or 8 N 2 equations. system is described by 6 N Since these equations are written in terms of 8N 9 unknowns (including the quantity N ) , a correctly posed problem must specify values of seven properly chosen quantities, as is the case in the two problems mentioned above. The explicit forms of the underflow, equilibrium, and overflow relations must now be developed. If p is defined as the mass of solution (containing components 1 and 2 ) adhering t o 1 pound of inerts in the underfloming L stream, we can write
+ + +
p=
(21k
+
+
+
x2k)/x3k
which, n-hen combined with xlk
+ + XZk
23k
= 1
gives the underflow relation: Z2k =
U(Xlk)= P/(P+1) - ZY,
(15)
If we specify that the overflowing V stream contains no inerts-i.e., that yak = O--vve can write the overflow relation as y2k = o(Ylk) = 1 - Y l k (16) Finally, if the equilibrium condition for the streams leaving the kt,h stage requires that the concentration of component 1 in the overflowing “clear” solution equal that of the entrapped
use Equation 19e but rather
solution in the underflowing L stream, we can write Ylk/ (Ylk
-k YZk) = x l k /
(zlk
+
and, using Equation 151,obtain Ylk=
E:(xlk)= % ( p +
I)/p
(17)
Embedding Formulation
As before, we are interested only in outlet values, so we introduce four new variables as follows:
:= xw
rw-1 (ym) E
TN ( y z ~ + i )
LN RN-i ( Y Z N )
RN (ym+i)
XW-i
LN--1
and apply them to Equations 13 and 14a, b with k = N :
+ VN (18%) (c)RN-i (c)+ yiN+iViv+i rhi @ ) R N( b ) + YWVN (1%) xz.v-iRiv-i(C) + hViv+i = QNRN( b ) + CVN (18c) RN-i ( c ) -tV N += ~ RN ( b )
TN-I
=
where b = yZN+l and c = YZN. These three equations are augmented by the following set taken from Equations 15, 16, and 17 evaluated a t k == N and N - 1:
.VIN+I
= 1-
b
c = 1 - Y1N xZN =
U[rN ( b ) ]
(19b) (19c) (19d)
(c)] (19e) Thus if it were necessary, the six secondary variables V N ,c, y ~ + 1 y, ~ X Z, N + ~ ,and Q N - ~ could be algebraically eliminated f r o p the above eight equations, giving two first-order nonlinear finite difference equations in r~ ( b ) and RN ( b ) subject to two initial conditions : XZN-I
ro ( b ) =
= UCm-1
210
= 1
-0
z2k)
RO( b ) = LO
The variable V N +is~not eliminated because it can be treated as a constant once its value is specified in a correctly posed problem. Computational Procedure
To simplify matters we eliminate VN between Equations 18a and 18c and solve for RN ( b ) , obtaining
RN(~)C = ( ~ ~ - i - c ) R i ~ - i ( c ) (+b - ~ ) T ' i v + i l / ( ~ C~)-
- $10 - 230
Also TO( c ) and Ro( c ) would be replaced by $10 and LO. 2. Pass a least squares polynomial P M ( ~of) order M through the rl (bl), . . ., rl (bm) values and another polynomial QM ( b ) through the R1 ( b l ) , . . ., RI (bm) values. 3. Let N = 2 and solve for r z ( b ) and Rz(b) by systematic trial and error for the several b's, using PM and &M evaluated at c to replace TI ( c ) and R1( c ) . 4. Repeat step 2 for the r z ( b i ) and Rl(bi) calculated in step 3. 5. Continue the process and obtain r3 (bi) and Rs (bi), etc. This cycle of calculation would stop when rN(b) matched the desired outlet value in the common design problem or when N matched the number of stages in the existing process problem. Numerical Example
The following problem is one of the most often solved in the literature on solid-liquid extraction (Badger and Banchero, 1955, p. 351). The given values are LO= 2850 pounds, xl0 = 0.281, $20 = 0.01754, VN+I = 1330 pounds, ~JZN+I= 0.985. It is desired to find N such that X ~ N= 0.04. The solute is oil, the solvent is benzene, and the inert component is meal. For this system p is 8 weak function of xlk, so a least squares polynomial of order 1 was constructed to fit the data: p = 0.489 0.464~1k
+
The values of b chosen were 0.24625, 0.4925, 0.73875, and 0.985, the last value being the one given in the problem. After calculating the four TN ( b ) and RN ( b ) for each stage, two least squares polynomials of order 1 were passed through the values; the polynomials being used to estimate rN--l(c) and RN-l(c) in the calculations for r ~ ( b and ) R N ( ~ ) The . results of the calculations are shown in Table 11. The values of 0.046 and 3020 obtained for rd(0.985) and R4(0.985) agree closely with the values of 3018 and 0.04 presented by Badger and Banchero (1955, p. 353), where the problem is solved graphically. Using a time-shared GE 235 computer (12microsecond fixed point add time) the calculations (including compilation) consumed approximately 5 seconds. The value with which d2 was compared was 0.01. Subsequent runs with smaller values showed no significant change in the results.
Table II. Solid-liquid Extraction Results
(20 ) Now the problem is similar to the constant flow case-namely, to find by systematic trial and error a value of rN ( b ) that will make d2 sufficiently small, where d comes from Equation 18b: d = rN @ ) R N( b )
+ YILNVN- rN-1 (C)RN-l ( c ) -
Outlet Values No. of Stages, N
Inlet Value.
1
0.24625 0.4925 0.73875 0.985
0.322 0.256 0.193 0.134
R N ( Y 2 N +I) 3276 3214 3156 3101
2
0.24625 0.4925 0.73875
0.310 0.233 0.159
3265 3193 3124
3
0.2465 0.4925 0.73875 0.985
0.303 0.221 0.140 0.062
3258 3182 3107 3035
4
0.24625 0.4925 0.73875 0.985
0.299 0.213 0.128 0.046
3254 3174 3096 3020
YZN+1
ylN+lVN+1
Thus, given a trial value of rN(b), Equations 19a, b, c, d, e , are used to calculate t'he five secondary variables Y ~ N ylN+1, c, X Z N , and ~ 2 ~ 7 - 1 . Then Equations 20 and 18a yield RN(b) and V X ,after which ti2 is calculated and compared to the previous best value and a decision is made whether to alter the step size. The over-all procedure can be outlined as follows: 1. Let N = 1 and solve for r1 ( b ) and R1 ( b ) by systematic trial and error for several choices of b, say bl, bz, . . ., b,, in the range of interest. Of course in calculating $20 we would not
VOL.
0
?" ( Y Z N + I )
NO. 3 A U G U S T 1 9 6 9
539
Summary
The concept of invariant imbedding when applied to countercurrent equilibrium stage problems yields equations which can be solved for outlet concentrations and flow rates as functions of the number of stages and some properly chosen input quantity. The advantages of this approach are at least threefold : The solution procedure is easily programmable, obviating the need for graphical procedures which are often subject to error; one solution yields outlet values for systems with any number of stages (and any value of the chosen input quantity) , thus providing the designer an opportunity to make a n optimal choice of design parameters; and the existing process problem conventionally solvable by uncertain timeconsuming trial and error methods (graphical or otherwise) now becomes 110 more difficult than the common design problem. Nomenclature
A = absorption factor ( L / K V ) b, bi = inlet concentration of V stream c
= Y2N
d
= measure of error = component index
i k K Lk
(YN+.I
or y 2 ~ + 1 )
= stage index = equilibrium constant = flow rate of L stream leaving
kth stage
No. of values of bi chosen
m
=
M N
= order of polynomial = total number of stages = polynomial of order M = polynomial of order M
Phi QM r
= outlet concentration of
R
= outlet flow rate of L stream = rN(b)/b
LN
L stream
Vk = flow rate of V sheam leaving kth stage = mass fraction of the ith component in L stream leaving kth stage = mass fraction of the i t h component in stream leaving Ytk kth stage
v
literature Cited
Badger, W. L., Banchero, J. T., “Introduction to Chemical Engineering,,’ McGraw-Hill, New York, 1955. Bellman, R., Kalaba, R., Wing, G., Proc. Natl. Acad. Sci. 46, 1646-9 (1960).
Coulson, J hi. Richardson. J. F., “Chemical Engineering,’’ Vol. 11, Pergamon Press, New York, 1956. Foust, A., et al., “Principles of Unit Operatiom,” p. 77, Wiley, New York, 1960. Koenig, D. hl., Chem. Eng. 74, 181 (September 1967). Lee. E. S.. IXD.ENQ.CHEM.FUNDAMENTALS 7. 164 (19681. Rlikkley, H. S., Sherwood, T. K., Reed, C. E., “A plied Mathematics in Chemical Engineering,” p. 324, %IcGraw-Hill, New Tork, 1957. Wing, G. XI., “Introduction to Transport Theory,” Wiley, New Tork, 1962. RECEIVED for review July 10, 1968 ACCEPTED March 21, 1969
SPATIAL DISTRIBUTION OF PARTICLES IN A SUSPENSION LYNN M. LEEK,’ PRADEEP DESHPANDE, AND JAMES R. COUPER Department of Chemical Engineering, University of Arkansas, Fayetteville, Ark. 72701 The distribution of solid particles in a two-phase system was studied, using a spectrophotometer and transparent phases of different absorptivities. The solid phase consisted of uranium-impregnated borosilicate glass cylinders. A gel composed of 85 volume % glycerol and 15 volume % benzyl alcohol was used as the liquid phase. The gelling agent was 1.5 weight % hydroxyethylcellulose. The beta distribution gave an excellent fit to the data, as evidenced by a mean P r ( X 2 2 C) by chance of 0.746. The parameters of the beta distribution are a function of the size of the solid particles.
H E X solid particles are added to a homogeneous liquid, Wthe physical properties of the liquid are changed. Some properties of the resultant suspension, such as density and heat capacity, are easily calculated as the weighted average of component properties. Other properties, such as thermal and electrical conductivity, are highly dependent upon the spatial distribution of the particles in the liquid and canriot be calculated accurately from solid and liquid conductivities and bulk particle volume fraction. Inaccuracies increase with increased difference between liquid and solid conductivities and with increased particle volume fraction. Inveqtigators (Bruggeman, 1935; Jakob, 1950; Jefferson, 1 Present address, Humble Oil and Refinery Co., Baytown, Tex. 77520
540
l&EC
FUNDAMENTALS
1955; Lees, 1898; Loeb, 1954; McAllister and Orr, 1964; Maxwell, 1904; Meredith and Tobias, 1960; Nahas and Couper, 1966; Rayleigh, 1892; Russel, 1935; Woodside, 1958) who have characterized particle distribution with well-defined arrays of cubes, spheres, or cylinders have witnessed this discrepancy in every case. Others (Baxley, 1966; Baxley and Couper, 1966; Fricke, 1924; Hamilton, 1960; Hamilton and Crosser, 1962; Reymond, 1964; Roblee, 1958; Tsao, 1961) who have more realistically characterized particle spatial distribution have found that their predictive expressions contain variables that are difficult to evaluate. However, investigators who have considered particle spatial distribution have not taken data specifically on particle spatial distribution or particle description. Instead, they have proposed various particle characterizations and have incorporated these into
