Ind. Eng. Chem. Res. 1989, 28, 793-803
793
Zhang, X. Particle Classification for Dilute Suspensions Using Inclined Settlers. M.S. Thesis, University of Colorado, Boulder, 1988.
de Sedimentation des Suspensions dans les Recipients Inclines.
Keijo J . Med. 1937, 8, 256-296. Ponder, E. On Sedimentation and Rouleaux Formation. Q. J . Exp. Physiol. 1925, 15, 235-252.
Received for review August 9, 1988 Revised manuscript received January 30, 1989 Accepted February 17, 1989
Schaflinger, U. Influence of Nonuniform Particle Size on Settling Beneath Downward-Facing Walls. Int. J. Multiphase Flow 1985, 1 1 , 783-796.
Analysis and Algorithms for Multistage Separation Processes Lakshmi N. Sridhar and Angelo Lucia* Department of Chemical Engineering, Clarkson University, Potsdam, New York 13676
Multistage equilibrium separation processes involving binary homogeneous mixtures are analyzed. For several sets of specifications, it is shown that these process models admit a unique steady-state solution without the need for a constant molar overflow assumption or other energy balance or phase equilibria simplifications. The steady-state model equations are analyzed directly, and the role of homogeneity in determining solution uniqueness is clearly identified. The analysis is constructive and gives insight into the behavior and development of algorithms for solving multistage separation process model equations in both the binary and multicomponent cases. To illustrate this, convergence results are established for direct substitution and certain Newton-like accelerated direct substitution methods. Also, a new generalized Broyden acceleration technique is developed and used in the context of a sum-rates algorithm to solve several example problems. Numerical results show that the proposed acceleration method is comparable t o Newton acceleration and clearly outperforms traditional Broyden acceleration. Other properties predicted by the analysis are also illustrated. 1. Introduction In a recent paper, Doherty and Perkins (1982) provide a very comprehensive review of research concerned with the mathematical analysis of homogeneous distillation processes. By homogeneous distillation, we mean one in which both the liquid and vapor on any stage of the column are stable with respect to phase splitting (see Van Dongen et al. (1983) and Lucia (1986)). The paper by Doherty and Perkins focuses on the issues of uniqueness and asymptotic stability of steady-state solutions to a variety of distillations. Although many dynamic models are discussed, it is the constant molar overflow (CMO) process that forms the basis for most of the useful qualitative results in this area, including those of Lapidus and Amundson (1950), Acrivos and Amundson (1955), Rosenbrock (1960, 1962), and Doherty and Perkins (1982). In particular, Doherty and Perkins consider CMO multistage columns involving binary homogenous mixtures and use both a simplified dynamic model (without energy balances) and classical linear stability analysis to prove that any distillation of this type has a unique and stable steady-state solution. Unfortunately, in many practical applications, the assumption of constant molar overflow can be a poor one. Furthermore, analysis of a dynamic model provides little knowledge of the convergence behavior of algorithms used for the steady-state simulation of multistage separators. Therefore, in this work, the steady-state model is analyzed directly and no simplifying assumptions are made with regard to energy balance effects. Constructive proofs are presented, and algorithms based on these proofs are implemented and compared with existing methods. This manuscript is organized in the following way. Useful perturbation and eigenvalue-eigenvector results for feed and temperature perturbations to a single-stage, isobaric flash process are presented in section 2. In section 3, a multistage analogue of the isothermal, isobaric (T,P)
* Author
to whom all correspondence should be addressed.
0888-588518912628-0793$01.50/0
process involving binary homogeneous mixtures is analyzed. Here the column temperature and pressure profiles are specified, and the stages of the separator are decoupled. By use of the single-stage results, the separation process model is put back together in the form of a classical fixed-point iteration, which is shown to be globally convergent to a unique solution in the two-phase region using the contraction mapping principle. Other, more practical specifications are analyzed in section 4. These specifications, which include fixing the reboiler and condenser duties or specifying the reflux ratio and total bottoms flow rate, are analyzed by constructing nonlinear maps using the implicit function theorem and subsequently showing that these maps are one-to-one and onto. In section 5, insights provided by the analysis are used to construct Newton-like accelerated direct substitution techniques for solving separation process model equations. Convergence and rates of convergence of these algorithms are also established. Many numerical examples are presented in section 6. The main conclusions of this research are that 1. binary separation processes involving homogeneous mixtures admit unique steady-state solutions under a variety of specifications; 2. for the specifications that were studied, the assumptions of constant molar overflow and simplified phase equilibrium are unnecessary, and solution uniqueness is solely determined by the homogeneous nature of the mixture; and 3. direct analysis of steady-state separation process models gives useful insights concerning the construction of reliable and efficient computer algorithms for solving these model equations. 2. Useful Single-Stage Results
In this section, useful perturbation results for singlestage, isobaric flash processes involving homogeneous mixtures are collected because they are needed to analyze the multistage separation process models that appear throughout the remainder of the paper. 0 1989 American Chemical Society
794 Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989
--
.
give the conditions My = y and (I - M)y = 0 , where y is the equilibrium vapor-phase composition. Thus, both matrices M and I - M, which are 2 X 2 in the binary case, have eigenvalues of 0 and 1 and share the same two eigenvectors; those eigenvectors are the equilibrium liquidand vapor-phase compositions x and y , respectively. These relationships, which can be expressed compactly in the form
"I
t
Mx=O
Figure 1. Flash process.
Consider the single-stage flash process shown in Figure 1. Assuming equilibrium stage behavior, one set of model equations for this process is given by /J- -
gv = 0
f-I - v = 0
(phase equilibrium)
(1)
(conservation of mass)
(2)
FHF- L H L - VHv + Q = 0
(energy balance)
(3)
where pL and pv are vectors of component chemical potentials in the liquid and vapor phases, respectively, and f , 1 , and v are vectors of feed, liquid, and vapor phase component molar flow rates. Furthermore, F , L , and V are the total molar flow rates of the feed, liquid outlet, and vapor outlet streams, respectively, and HF, HL,and H V are the corresponding stream molar enthalpies. Finally, Q is the heat duty to the flash drum. T and p are the equilibrium temperature and pressure, respectively, and n, is the number of components in the mixture. Note that the quantities pL, pv, f , I , and v are vectors of length n,. At constant temperature and pressure, the energy balance can be removed from the analysis. Differentiation of the phase equilibrium equations gives the condition
V2G LAl - V2G 'Av = 0
(4)
where V2GL and V2GV are the Hessian matrices of the Gibbs free-energy function for the liquid and vapor phases, respectively. For feed perturbations, say A f , conservation of mass implies the simple relationship
AI = h f - Av
(5)
Substitution of eq 5 in eq 4 and some algebra give Av = (V2GL
+ V2GV)-lV2GLAf
(6)
Furthermore, the matrix V2GL + V2GV is guaranteed to be invertible at any nonazeotropic equilibrium point if the mixture is homogeneous (see Lucia (1986)). Finally, the definition M = (V2GL + V2GV)-'V2GL and eq 5 imply that
AV = MAf
A I = (I - M ) h f
(7)
which gives the changes in vapor and liquid component molar flow rates as a function of the feed perturbation at the given temperature and pressure. Equation 7 is one of the results that we will use many times in the subsequent sections of this paper. Straightforward application of the Gibbs-Duhem equation for the liquid and vapor phases shows that the matrices M and I - M both have eigenvalues of 0 and 1 in the binary case. To see this, note that the Gibbs-Duhem equation for the liquid phase implies that M x = 0, which in turn implies that (I - M ) x = x , where x is the equilibrium liquid-phase composition. Similarly the GibbsDuhem equation for the vapor phase and the fact that 1 - M = 1 - (V2G L
+ V2G V)-1V2G L (V2G
=
+ V2G V)-1V2G
(8)
(I - M)x = x
My=y (I - M)Y = 0
(9)
will also be used frequently in subsequent sections of this paper. In order to include energy balance effects in the analysis of the multistage process, the changes in equilibrium phase molar flow rates as a function of arbitrary changes in equilibrium temperature for a single-stage are also needed. These results, which are given in Lemma 4 in Lucia (1986, p 1766), are simply restated here for convenient reference. For an arbitrary scalar perturbation in equilibrium temperature, say AT, the changes in equilibrium vapor and liquid molar flows at constant pressure are given by
+ V2GV)-lAH AI = - ( A T / T ) ( V 2 G L + V2GV)-'AR AV
= (AT/T)(v~G~
(10)
where A H = Hv - HL and HLand Hv are vectors of component partial molar enthalpies for the liquid and vapor phases, respectively. Equation 10 will be used repeatedly in analyzing multistage processes involving condenser and reboiler duty specifications. Finally, we note that the results in this section carry over to the multicomponent case with little or no modification. In particular, eq 7,9, and 10 apply to the multicomponent case with no modification. Furthermore, it is possible to establish that all of the eigenvalues of M and I - M must lie in the closed interval [0,1] for multicomponent homogeneous mixtures. This fact can be proved by using a simple similarity transformation. 3. Multistage Binary Separation Processes In this section, we show that any binary homogeneous separation process, in which the temperature and pressure profiles are specified, has a unique steady-state solution. As highlighted in section 1, to establish the result, we 1. decouple the stages of the separator, 2. use the single-stage results from the previous section to reassemble the multistage model in the form of a fixed-point iteration, and 3. show that this fixed-point iteration is globally convergent to a unique solution in the two-phase region using the contraction mapping principle (Ortega and Rheinboldt, 1970, p 120). Consider then the multistage separation process shown in Figure 2, in which the temperature and pressure at each stage are specified. The stages are numbered from top to bottom, and although only one feed is shown in the figure, the results that are established remain valid for any number of feeds and/or separator equipped with a total or partial condenser. The stages of the separator can be decoupled by "tearing" the vapor input stream to each stage, 1 through ns - 1. Application of the straightforward perturbation analysis in section 2 (i.e., eq 7 ) to each stage in the separator gives the equations j = 1, ..., ns (11) Avj = Mj[Avj+l + AI,-J Alj = (I - Mj)[Avj+l + AIjJ
j = 1, ..., ns (12)
Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989 795 Chart I 0 M3
0
1
... 0
0 ... 0
T
x, Y
Figure 3. Tzy diagram for homogeneous binary mixture.
v T8 *Pa
18
Q8
-y z
Figure 2. Multistage separation process.
with appropriate modification for the end stages (Le., AZO = 0 and Avns+l = 0 ) . Elimination of AZ,, j = 1, ..., (ns - 1)by using eq 12 in eq 11shows that the following fixed-point iteration, in the variables vi, j = 2, ..., ns, can be constructed
Xk+’ = G(Xk) (13) where X = (vZT,vST,..., vWqTand where G’, the Jacobian matrix of the fixed-point iteration, is defined in Chart I. Each of the block matrices Mj and I - M j is n, X n, or 2 X 2 in the binary case. Note that eq 13 is a linear equation in the case of binary mixture and a nonlinear equation in the multicomponent case. To show that the fixed-point iteration defined by eq 13 converges to a unique solution, we apply the contraction mapping principle (Ortega and Rheinboldt, 1970, p 120). One way to do this is to show that the spectral radius of the Jacobian matrix, p(G’), satisfies
p(G’) = maxl,isn lAi(G’)l 1 (15) To establish this bound, first note that, under the assumption that a solution exists, the compositions on any stage of the separator must satisfy yj = aj-’xj-’+ (1 - aj-l)y,-l j = 2, ..., ns (16) xj =
Oj+lXj+l
+ (1 - fij+l)Yj+l
I’= 1,
(see Figure 3). For binary, heterogeneous mixtures, ai-’ and Pj+’ need not be unique at stage temperatures below the heterogeneous azeotropic temperature (see Figures 11 and 17 in Van Dongen et al. (1983)),and this observation clearly indicates that homogeneity is a necessary condition for guaranteeing solution uniqueness. Define the transformation P by
.*a,
(ns - 1) (17)
where ajw1, fij+l E (0,l). For binary, homogeneous mixtures, the values of aj-l and fij+l in eq 16 and 17 are unique
0 0
0 y3
0
... 0 ... 0
wp
0
... 0
w3 ... 0
0
0
P=
where wj = yj-’ - ( yj-lTyj/yjTyj )yj,j = 2, ..., ns and y j denotes the vector of vapor-phase compositions on stage j in the separator. Note that P is invertible, and this gives the similarity transformation P-IG’P =
[
A
@
0
0
]
where the matrix A is given by
J where aj, Oj E (0,l) for each j as before and where each block of the matrix P-lG’P is (ns - 1) X (ns - 1). Note that A is a nonnegative matrix. Equation 19 shows that the eigenvalues of G‘ are the same as the eigenvalues of the two (ns 1) X (ns - 1) diagonal blocks of P-’G’P. Thus, one-half of the eigenvalues of G’ are zero and p(G’) = p ( A ) . Application of
-
796 Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989 energy balance equations
Lemma 2.5 in Varga (1962, p 31),
gives the desired upper bound on p(A) because eq 20 shows that each column of the nonnegative matrix A has a sum that is strictly less than one for cyI, bj E (0,l). Consequently, d G ' ) = p(A) < 1 (22) We are now in a position to use the contraction mapping theorem, which is stated here for reference. Contraction Mapping Theorem (Ortega and Rheinboldt, 1970, p 120). A mapping G:D C R" Rn is contractive on a set Do C D if there is a k < 1 such that IlG(x) - G(y)II Ikllx- yII for all x, y E Do. Suppose G is contractive on a closed set Do C D and GDo C Do. Then G has a unique fixed point. To apply this theorem to the fixed-point iteration defined by eq 13, let Do be the set intersection of the twophase regions for each stage in the column and set k = p(G'). If a two-phase solution to the column model exists, then Do must be a closed set. Furthermore, the condition p(G') < 1 and Lemma 2.2.8 (Ortega and Rheinboldt, 1970, p 44) establish that G is contractive on Do. Consequently, the iteration defined by eq 13 converges to a unique fixed point in Do, and any separator involving a binary homogeneous mixture and fixed temperature and pressure profiles has a unique solution, provided those specified temperature and pressure profiles admit any solution.
-
4. Other Specifications Specifying the column temperature and pressure profiles has both analytical and computational value. As shown in the previous section, it provides a useful way of analyzing the number of solutions to binary separation processes using classical fixed-point theory. Furthermore, the analysis is constructive, and as we show in subsequent sections, this is important in understanding the convergence behavior of certain equation-tearing algorithms for multistage separation calculations (e.g., the Sujata sumrates algorithm (Sujata, 1961)). However, specifying the temperature and pressure profiles is neither the usual way nor the only way to fix the operation of a multistage separator. It is more common to specify the adiabatic operation of all intermediate stages in the separator and to fix two additional end conditions such as the reboiler and condenser duties or the reflux ratio and total bottoms flow rate. Other end specifications are possible. In this section, we study the steady-state behavior of binary multistage separators involving the specifications described above (i.e., adiabatic intermediate stages with fixed reboiler and condenser duties or fixed reflux ratio and total bottoms flow rate). To do this, nonlinear maps from a set of specifications for which the solution is unique to the set of specifications of interest are analyzed by using the implicit function theorem. Uniqueness is preserved if the nonlinear map is one-to-one or, equivalently, if the Jacobian matrix of the map is invertible throughout the domain of interest. We begin with the specifications of reboiler and condenser duties. 4a. Reboiler and Condenser Duties. Define the set of solutions in which the temperature and pressure profiles in the separator are fixed as the set of ( T Q ) solutions and the set in which the heat duty to each stage and pressure profile are fixed as the set of (Q,P)solutions. Note that the stage energy balance equations define a nonlinear map from the set of ( T Q )solutions to the set of ( 8 9 )solutions (see Figure 4). It is this map whose Jacobian matrix must
/
)
(0,P)solutions
Figure 4. Nonlinear energy balance mapping.
be shown to be invertible throughout the domain of interest. Consider the binary separation process shown in Figure 2 . The energy balance for the j t h stage in the process is given by
Ej(lj-1, vj+l,l j , vj, TI, 7 ' 2 , TnJ = F j H r + Lj-lHkl + Vj+1Hy+1- LjHiL - V,H.' I 1 *a*,
+ Qj = 0 (23)
ciljj,
where Lj = V, = Civij,and all other quantities are defined as before. Application of the implicit function theorem to eq 23 gives AQj = -Lj-l(dH~-l/dTj-l)ATj_lns
Hk1[ E (bLj-l/aTk)ATk]- Vj+l(dHy+l/dTj+l)AT,+l k=l
Hy+1 [ (d v;+l/d Tk)A Tk] + Lj(dHjL/a Tj)AT; -t k
H;L[c(dLj/dTk)ATk] + Vj(aHj"/dTj)ATj + k
HjV[C(dvj/aTk)ATk](24) k
Note that the term dLj/dTkand dVj/dTkappear in eq 24 because of the material balance interlinking between adjacent stages in the separator and that a change in any stage temperature affects the total liquid and vapor flow rates throughout the separator. From eq 24, it is readily apparent that in order to determining the quantity AQj, it is necessary to determine dLj/dTkand dVj/dTk. Using the single-stage results (i.e., eq 7 and 10) in section 2, the Neumann lemma (Ortega and Rheinboldt, 1970, p 45), the inverse nonnegativity of the matrix V2GL+ V2GV (Lucia, 1986, p 1768), and certain continuity properties, it is possible to establish the following conditions aLj -aT, =(
eo j l k >O j > k
aVj -dTk = (
>O j 5 k k
(25)
provided there are no retrograde effects or supercritical components throughout the separator. A derivation of the conditions given by eq 25 and all of the associated assumptions are presented in detail in Appendix A; however it is important to note that the property of inverse nonnegativity can only be guaranteed for homogeneous mixtures and, again, illustrates the pivotal role of homogeneity in determining solution uniqueness in the binary case. We are now in a position to show that the nonlinear map from the set of (2'9)solutions to the set of (QQ) solutions preserves uniqueness. Let Q = E ( T ) denote the set of stage energy balance equations for the separator. Application of the implicit function theorem gives A Q = JAT (26) where A Q T = (AQ1,AQ2, ..., AQnS),ATT = (AT1,AT2, ...,
Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989 797
ATn ), and the elements of the Jacobian matrix, Jjk = dQ.JbTk,are defined through eq 24 for j = 1, 2, ..., ns. rfo establish that the Jacobian matrix in eq 26 is invertible, note that aQj/dTk = -sj-l,k[Lj-l(aHfil/aTk)] - Hfil(aLj-i/aTk) s j + l , k [ vj+l(aHT+l/ a Tk)1 - Hy+l(dVj+l/aTk) + bjk[Lj(dHjL/dTk)+ vj(aHjV/dTk)]+ HjL(aLj/aTk)+ Hjv(a vj/dTk) (27) where
6jk
by application of eq 24 and 29. This, together with eq 30, 31, and 32, implies the condition
aQk/aTk
> I(&!k-i/aTk)+ (aQk+i/aTk)l
which shows that the matrix S is diagonally column dominant and therefore invertible. Moreover, because the quantities aQj/dTkfor j # k - 1, k , k + 1 are relatively small with respect to the elements in the kth column of S, it follows that p(S-'T) < 1 and that I + S-'T is invertible. Therefore,
is the Kronecker delta defined by J-1
,a,
1 j = k = 0 j # k
Choose a reference state such that HjL < 0 Hjv > 0 j = 1, 2,
(28)
..., ns
(29)
Use of eq 24 and 29 in eq 27 with j = k and the fact that dH/aT > 0 shows that
aQk/aTk = -Hk-l(aLk-l/aTk) - Hx+i(avk+i/aTk)+ Lk(dHkL/dTk)+ vk(aHkV/aTk)+ HkL(dLk/dTk) HkV(dvk/aTk)> 0 (30)
+
Note that each term in eq 30 is positive. Furthermore, it is possible to show that
dQk+l/dTk -Lk(aHkL/dTk)- HkL(aLk/aTk)Hx+z(avk+z/aTk)+ Hk+,(aLk+l/aTk) Hx+l(aVk+l/aTk)< 0 (31) and dQk-l/dTk = -Hk-z(aLk-,/aTk) vk(aHkV/dTk)- HkV(dVk/dTk)+ Hk-l(aLk-,/aTk) + Hx-,(aVk-l/dTk) < 0 (32)
A detailed derivation of the conditions given by eq 31 and 32 can be found in Appendix B. Finally, using arguments similar to those in Appendix B, it can be shown that for j # k - 1, 12, and k + 1,
(37)
= (I + S-'T)-'S-'
(38)
and this implies that the map from the set of (TQ) solutions to the set of ( Q , P )solutions perserves uniqueness. 4b. Reflux Ratio and Total Bottoms Flow Rate Specification. For separation processes like distillation, it is common to define the operation of the column by fixing the pressure on all stages, by fixing the heat duty to all intermediate stages (usually Qj = 0 for j = 2, ..., (ns - 1)but this is not necessary), and by specifying the reflux ratio and total bottoms (or distillate) flow rate denoted by R and B, respectively. In fact, this set of specifications is often referred to as the standard set of specifications. In this section, it is shown that this standard set of specifications gives a unique solution. Uniqueness is established by defining a nonlinear map from the set of (8s) solutions to the set of solutions for which the reflux ratio, bottoms flow rate, and heat duty to all intermediate stages are fiied and by producing a contradiction if the Jacobian matrix of that map exhibits singularity. Define the latter set of solutions as the set of ( R , B ) solutions. Let AQj = 0, j = 2, ..., (ns - 1). Also let AR = 0 and AB = 0. Since AB = 0, the total mass balance around the separator gives
AVl = 0 (39) because the feed to the separator is fixed. Straightforward differentiation of the defining equation for the reflux ratio
R = Ll/Vl
(40)
A R / R = AL1/L1 - AVl/Vl
(41)
implies the condition which together with eq 39 and the fact that AR = 0 gives AL1=0
are relatively small quantities in comparison to the tridiagonal terms of the Jacobian matrix, even though they can be either positive or negative. Therefore, the Jacobian matrix in eq 26 can be written in the form
J = S + T = S(I + S-lT) (34) where S is an invertible tridiagonal M matrix and p(S-lT) < 1. The matrix S is invertible because it is diagonally column dominant. To see this, consider the energy balance around the envelop enclosing stages k - 1, k, and k + 1, which is given by Lk-2Hk-2 - vk-1Hx-1 - Lk+lHk+l + vk+2Hx+2 + Qk-i
Qk + Qk+i = 0 (35) Straightforward differentiation of eq 35 with respect to Tk shows that
a(Qk-1
+ Qk + Qk+l)/dTk = -Hk-z(aLk-,/aTk) Hl-i(avk-i/aTk) Hk+i(aLk+i/aTk)Hx+2(aVk+z/aTk)> 0 (36)
(42)
Moreover, the mass and energy balance around stage 1 implies that both
AV2 = 0
(43)
AQ1 = 0
(44)
and If AQns = 0, then we are done because this means that the Jacobian matrix of the nonlinear map from the set of (QQ) solutions to the set of (R,B) solutions is nonsingular. Therefore, suppose AQ,, # 0. Equations 39 and 42-44 together with the fact that AB = 0 and the energy balance equation for the envelop enclosing stages 2 through ns imply the condition AQ, # 0 for a t least one j = 2, ..., (ns- 1) (45) However, this is a contradiction since the specifications demand that all intermediate stage heat duties remain fixed. Consequently, AQ, = 0, the map from the set of ( Q P )solutions to the set of (R$) solutions is nonsingular, and the standard set of specifications admits a unique solution.
798 Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989
5. Algorithms In this section, we show that the preceding analysis can be used to help understand the behavior of certain existing algorithms for steady-state multistage separation process calculations. It can also be used to construct new algorithms. In particular, we show that the multistage TP analysis in section 3 gives insight into the convergence of the inner loop of sum-rates methods for absorbers, strippers, and distillations of wide boiling mixtures (see Sujata (1961)). Furthermore, it can be used to construct new reliable and efficient Newton-like acceleration methods for sum-rates algorithms. The analysis in section 4, on the other hand, can be used to understand the convergence behavior of the outer (temperature) loop in the sum-rates method (see Sujata (19611, Friday and Smith (1964), and King (1971)). 5a. Sum-Rates Algorithms. The sum-rates method for absorber, stripper, and wide boiling distillation calculations was originally proposed by Sujata in 1961 and is certainly well-known by now. The reader is referred to the paper by Sujata (1961) as well as the book by King (1971, p 529) for a description of the algorithm. For the purposes of our discussion, it suffices to state that these algorithms use fixed temperature and pressure profiles in an inner loop and solve the phase equilibrium and mass balance equations for the column (in any number of ways and to a variety of desired accuracies). The outer loop consists of using the energy balance equations and stage heat duty specifications to adjust the temperature profile in the inner loop. Several choices of inner loop tear stream variables are possible, but certainly one set of variables is the vapor component molar flow rates for each stage j = 2 , ..., ns. Thus, the inner loop of the sum-rates method can be written in the form of the fixed-point iteration given by eq 13, where the associated Jacobian matrix, G ’ ( X ) ,is given by eq 14. The outer loop can be handled by using eq 26 and 34. It is convenient to interpret direct substitution and accelerated direct substitution methods in the context of solving the system of equations F ( X ) = X - G ( X )= 0 (46) The values of Xk+l are defined by the iteration Xk+’= Xk - Jk -‘F(Xk) (47) where the Jacobian matrix J k is given by Jk = 1 - G ’ ( X k )
(48)
Approximation of the matrix Jk by the identity matrix (by assuming G’ = 0 ) at each iteration corresponds to direct substitution, which is usually characterized as extremely reliable but slow. Evaluation of the Jacobian matrix Jk at each iteration defines Newton acceleration and gives a local quadratic rate of convergence; however, in this context, Jacobian evaluations can be expensive because derivative information is most commonly available in implicit form. Other acceleration techniques such as Wegstein (1958) acceleration, dominant eigenvalue methods (Orbach and Crowe, 1971; Crowe and Nishio, 1975), and traditional Broyden acceleration (Boston and Britt, 1978) can also be used to speed convergence. 5b. Direct Substitution. The analysis in section 3 shows that direct substitution of the inner loop is globally convergent to a unique solution for any well-defined, binary homogeneous multistage separation process in which the temperature and pressure profiles are specified. The analysis is not restricted to stripper and absorber calculations but does show that, when applied to binary mixtures, the sum-rates method is globally convergent.
Furthermore, the Linear Convergence Theorem shows that the rate of convergence is R linear with a convergence factor equal to p(G’(X*)),where X * denotes the fixed point. Proof of this theorem can be found in Ortega and Rheinboldt (1970, p 301) and requires no modifications for this application. Unfortunately, the convergence results for direct substitution do not carry over to the multicomponent homogeneous case. Simple examples can be constructed for which p ( G ( X * ) )> 1and thus show that direct substitution cannot be guaranteed to converge in the multicomponent case without additional restrictions. 5c. Newton Acceleration. The analysis in section 3 can be easily extended to show that Newton-accelerated direct substitution is convergent to a unique solution by applying Lemma 10.1.7 in Ortega and Rheinboldt (1970, p 304) to the Newton fixed-point iteration defined by eq 47. In this case, we define GN(X)as the right-hand side of eq 47, and straightforward differentiation of eq 47 shows that G’N(X) = I - ( I - G’)-’[I- G ’ ] = 0 (49) because G’ is a constant matrix in the binary case. In practice, this means that Newton’s method converges in one acceleration step in the absence of rounding error. 5d. Quasi-Newton Acceleration Methods. Because of the fact that p(G’(X*))> 1 for direct substitution can occur, we have also considered a variety of Newton-like acceleration methods for sum-rates algorithms. In this section, it is shown that the analyses in sections 2 and 3 suggest quasi-Newton acceleration methods that have outstanding convergence characteristics in the binary and ternary cases. In order to approximate J - l , we build approximations to the blocks M,, j = 1, 2, ..., ns by using a generalized Broyden updating strategy and by constructing G’ according to eq 14. Each block M, is subject to the GibbsDuhem and secant constraints given by eq 9 and 11, respectively. Thus, the following least change secant updating problem (see Dennis and Schnabel (1979)) can be defined min //AJ- A,IIF:A,Z, = R,
(50)
where A, and A, are successive approximations to M, = (V2GIL + V2Glv)-’V2GJL and where Z, and R, are n, x 3 matrices defined by Z, = [ v I k +I,k+’, l, sik] R, = [vik+l, 0, y I k ] (51) where s,k = (Av,k-l + AIJ-lk)and ylk= AvJk. The solution to this problem is given by
A, = A,
+ (R,
-
A,Z,)(ZJTZ,)-’ZJT
(52)
Those familiar with quasi-Newton updating will recognize the update defined by eq 52 as a generalized Broyden formula. It is a least change rank p correction to A, and reduces to the usual Broyden update for p = 1. In our applications, p = 2 or 3 in the binary or multicomponent case, respectively. For binary mixtures, Newton and generalized Broyden accelerations are identical. This is true because each block of G ’ ( X )is constant and uniquely determined by its associated pair of Gibbs-Duhem equations, and thus the Jacobian matrix produced by the generalized Broyden update is the true Jacobian matrix. For ternary mixtures, it is possible to show that the generalized Broyden acceleration technique is quadratically convergent. In this case, each block of G ’ ( X ) ,while not constant, is uniquely determined by its associated pair of Gibbs-Duhem equations and secant condition. Consequentlv, the matrix
Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989 799 Table I. Example Problem: Wide Boiling Binary Mixture Problem Definition vapor ideal condenser partial liquid ideal feed sat liq, stage 2 no. of 4 (1) n-butane, 40 kmol h-' (2) n-decane, 30 kmol h-' stages Dressure 4.013 X 10' Pa starting values and solution vo, kmol h-' v*, kmol h-' stage T,K 1 334 25 10 35.68 0.0494 2 370 25 10 36.84 0.0545 25 10 11.48 0.6215 3 400 3.53 0.5695 4 430 25 10
produced by the generalized Broyden update converges to the true Jacobian matrix in the limit, and the asymptotic rate of convergence is quadratic. 6. Preliminary Numerical Results In this section, some preliminary numerical results are presented that compare direct substitution and a variety of Newton-like acceleration techniques on several example problems. Also, some of the properties predicted by the analysis are illustrated. 6a. Problem-Solving Results. For the purpose of illustration, the convergence methods are compared in the context of an inner loop of a sum-rates method, in which the temperature and pressure profiles in the separator are specified. Thus, the corresponding iteration is defined by eq 47. The acceleration techniques studied include Newton and traditional Broyden accelerations, as well as the new generalized Broyden acceleration technique developed in this work. These convergence methods are defined as follows: 1. Direct Substitution (DS). Jk = I for each iteration k in eq 48. 2. Broyden Acceleration (BA). The matrix Jk-l is approximated directly using the inverse Broyden update and traditional secant information. 3. Newton Acceleration (NA). Each block of the matrix G'(X) in eq 14 is evaluated using analytical derivative information. The Jacobian matrix is computed by using eq 48. 4. Tridiagonal Newton Acceleration (TNA). Only the tridiagonal blocks of G'(X) in eq 14 are evaluated. 5. Generalized Broyden Acceleration (GBA). Each block of G'(X) is approximated by using eq 52. 6. Tridiagonal Generalized Broyden Acceleration (TGBA). Only the tridiagonal blocks of G'(X) are approximated by using eq 52. Note that, for methods NA, TNA, GBA, and TGBA, the matrix G'(X) or its tridiagonal approximation must be constructed and the corresponding Jacobian matrix must be inverted, whereas for traditional Broyden acceleration the inverse Jacobian matrix is approximated directly. There is also a provision in our computer program that allows any specified number of direct substitution steps to be used prior to acceleration or between successive acceleration steps. The foregoing convergence methods were tested on a set of separation problems involving both wide boiling and narrow boiling binary and multicomponent mixtures. Problem definitions and solutions can be found in Tables I-V. All physical properties were calculated using the data and models given in Prausnitz et al. (1980). For Broyden acceleration, two initial direct substitution steps were used, and subsequent acceleration was used every fifth iteration. This strategy of delaying has been recommended by others
Table 11. Example Problem: Narrow Boiling Binary Mixture Problem Definition vapor nonideal condenser partial liquid nonideal feed sat liq, stage 4 8 (1)methanol, 40 kmol h-' no. of (2) water, 60 kmol h-' stages pressure 1.013 X lo5 Pa starting values and solution y o , kmol h-' v *, kmol h-' stage T,K 1 2 3 4 5 6 7 8
339.05 339.71 340.98 343.61 344.68 347.98 356.77 367.75
3.28 211.93 211.50 190.17 180.62 164.96 113.27 40.90
0.91 10.14 17.50 32.54 40.93 57.79 99.95 160.98
37.86 217.92 214.49 194.16 184.62 164.96 113.27 40.90
0.91 10.24 19.50 35.53 40.93 57.79 99.95 160.98
Table 111. Example Problem: Wide Boiling Ternary Mixture Problem Definition vapor ideal feed sat liq, stage 2 liquid ideal (1) propane, 40 kmol h-' no. of stages 8 (2) n-butane, 30 kmol h-' (3) n-decane, 30 kmol h-' pressure 4.013 X lo5 Pa condenser Dartial starting values and solution y o , kmol h-' v*, kmol h-' stage .'7 K 1 2 3 4 5 6 7 8
310 350 370 390 420 440 470 495
25 25 25 25 25 25 25 25
20 20 20 20 20 20 20 20
10 10 10 10 10 10 10 10
39.98 40.38 4.96 2.20 0.93 0.33 0.15 0.04
29.24 30.22 8.78 8.30 6.72 3.93 2.73 1.15
0.02 0.46 0.22 0.39 0.87 0.97 1.91 2.52
Table IV. Example Problem: Narrow Boiling Ternary Mixture Problem Definition vapor non idea 1 feed sat liq, stage 1 liquid nonideal (1)methanol, 25 kmol h-' (2) benzene, 40 kmol h-' no. of stages 6 (3) n-hexane, 35 kmol h-' Dressure 1.00 x io5 Pa starting values and solution stage T,K vo,kmol h-' v*, kmol h-' 2 3 4 5 6
330.0 332.0 334.0 336.5 339.0
15.0 7.5 2.0 0.5 0.1
10.0 7.5 2.5 1.0 0.5
12.5 7.5 2.5 1.0 0.5
2.83 1.22 0.67 0.37 0.15
2.20 1.17 0.82 0.65 0.41
2.46 1.28 0.88 0.69 0.43
(see Crowe and Nishio (1975),Westerberg et al(1979), and Michelsen (1982)) because the corrections in the unknown variables have a natural tendency to align with the eigenvector associated with the dominant eigenvalue of the Jacobian matrix. This causes near collinearity of the steps and usually results in secant information that is of little use. For Newton and generalized Broyden accelerations, on the other hand, acceleration was used on every iteration. Finally, the convergence tolerance for the loop was IlVk" - vkl12 < and all calculations were performed in double-precision arithmetic on an IBM 4341 computer. A summary of the results for five separation problems is given in Table VI. The results in Table VI show that the proposed generalized Broyden acceleration technique is competitive with Newton acceleration. It usually requires the same or only a few iterations more than Newton acceleration and clearly outperforms traditional Broyden acceleration.
800 Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989 Table V. ExamDle Problem: Wide Boiling Four-ComDonent Mixture Problem Definition vapor ideal pressure 4.013 X lo5 Pa liquid ideal condensor partial 10 feed sat liq, stage 2 no. of stages
feed
(1) propane, 25 kmol hW1 (2) n-butane, 25 kmol h-' (3) isobutane, 25 kmol h-l (4) n-decane. 25 kmol h-l
starting values and solution 330 350 370 385 400 420 440 460 480 495
1
2 3 4 5 6 7 8 9 10
v*, kmol h-I
vo, kmol h-'
T,K
stage
20 20 25 25 25 25 25 25 25 25
30 30 20 20 20 20 20 20 20 20
20 20 10 10 10 10 10 10 10 10
24 24 24 24 24 24 24 24 24 24
24.99 25.09 2.57 0.97 0.42 0.18 0.08 0.03 0.02 0.01
24.56 24.74 6.07 4.86 4.08 3.28 2.30 1.54 0.92 0.37
24.78 24.92 5.03 3.46 2.55 1.85 1.20 0.76 0.44 0.18
0.13 0.46 0.22 0.28 0.39 0.61 0.82 1.07 1.37 1.18
Table VI. Comparison of Various Convergence Methods iterations mixture
no. of stages
c4, G o MeOH, H,O cs, c4, ClO
4
8 8 6 10
MeOH, C6H6, n-C6Hl, cs,n-C,, i-C4, c,o
a
BA -
DS 9 50" 35 13
9 29 26 13 28
44
NA
TNA
GBA
TGBA
2
4 8 6 5 9
2 3 7 5 13
4 8 7 5 14
3 4 4 4
Failed to reach convergence tolerance in 50 iterations.
Table VII. Coefficients of Convex Combinations stage CY P 1 0.0444 2 0.0570 0.4401 3 0.1176 0.8170 4 0.8931
perturbations AT, = 0.01 K for k = 1-4, the Jacobian matrix in eq 26 is given by
In order to reduce the storage requirements associated with both Newton and generalized Broyden accelerations, we studied the effect of approximating G ' ( X ) ,and therefore J(X), by its block tridiagonal part. The results in Table VI show that this is, in general, a useful trade-off, especially for problems involving wide boiling mixtures. Usually only a few additional iterations are required to reach reasonable accuracy. 6b. Illustration of Properties Predicted by Analysis. Here we show that the bounds on the spectral radius given by eq 21 are easily satisfied by separation process problems involving binary mixtures. The first example problem in Table VI is used for illustrative purposes. For this problem, the coefficients of the convex combinations given by eq 16 and 17 are shown in Table VII. The corresponding matrix A (see eq 20) is given by
From this, it is easy to verify that the tridiagonal part (S) of the Jacobian matrix is a diagonally column dominant M matrix and is invertible. Furthermore, for J = S + T, the eigenvalues of S-'T are (-0.0589,0.0183, -1.849 X 9.368 X lo*). Consequently, p(S-'T) = 0.0589, and the Jacobian matrix given by eq 54 is invertible.
A =
1
0.0249 0.0036 0.0017
0.9430 0.0104 0.0050
0 0.8824 0.0126
]
-
69.75 J = lo2 -29.02 -0.824 -0.739
-
(53)
The eigenvalues of A are (0.1570, 0.5144 + 0.0761i, 0.5144 - 0.0761i) and thus p ( h ) = 0.5905. Note that two of the
eigenvalues of A are complex. It is easily verified that 0.0302 = minlSJSn (Chi,)< p(A) < 1=1
n
maxlSJSn(CA,) = 0.9584 1=1
Finally, it is illustrated that the Jacobian matrix in eq 26 can be written in the form given by eq 34 and that p(S-'T) < 1. For problem 1 (see Table I) and temperature
-
-75.44 306.36 -230.91 0.0045
0.005 -32.25 197.91 -165.65
-4.56 -4.56 -26.73 158.37 -
(54)
Acknowledgment This work was supported by the Office of Basic Energy Science, U.S. Department of Energy, under Grant DEFG02-86ER13552.
Nomenclature A = quasi-Newton matrix approximation E = energy balance function f = vector of feed component molar flow rates F = total feed flow rate G L, G v, G = total Gibbs free energy of liquid phase, of vapor phase, of fixed-point function H , H = molar enthalpy, partial molar enthalpy J = Jacobian matrix I = vector of liquid component molar flow rates L = total liquid flow rate M = operator P = pressure Q = heat duty R = real vector space, right-hand side of matrix constraint equation T = temperature v = vector of vapor component molar flow rates V = total vapor flow rate x = vector of liquid mole fractions X = unknown variables
Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989 801 y = vector of vapor mole fractions or right hand side of secant
condition Z = left-hand side of matrix constraint equation
tk
= Apflk
-tk = - A v k X k
+ VkGk - Yk)
(A6)
+ Lk(& - xk)
(A7)
Greek Letters a,fl = coefficients of convex combination A = nonnegative matrix 6 = Kronecker delta A = eigenvalue p = vector of chemical potentials
where AT, > o implies AVk > 0. Equations A6 and A7 complicate the analysis because of the presence of the quantities and 4 . However, because, for example, t k is a positive vector, the approximation tk A p f l k (A8)
Subscripts i = component index j = stage index
must be qualitatively correct; otherwise, a contradiction that t k is not a positive vector would occur. The same is true for ? .
Superscripts F = feed k = iteration counter L = liquid V = vapor
-tk
Appendix A: Effect of Temperature on the Total Vapor and Liquid Flow Rates in a Multistage Separator In this appendix, we show that the conditions given by eq 25 are satisfied in the absence of any retrograde effects throughout the separator. Application of eq 7 and 10 to each stage in the separator and some algebra give the equation A v = (I - G')-lr
(A9)
-AvkXk
Consequently, the second term in eq A6 and A7 can be dropped, simplifying the analysis without affecting its qualitative outcome. However, while eq A8 and A9 can be used for qualitative purposes, they will in general give poor quantitative approximations. Let ATk > 0 and ATj = 0 for j # k. From eq 9,17, A2, A8, and A9, it follows that
I
0
0 Yh -(1 - @k+l)Yk+l -(1 - @k+l)@k+lYk+Z
r=
(AI)
where A v = ( A v , ~ AvQT, , ..., A v , ~ ) ~where , G' is defined by eq 14, where r is the vector
Because p(G') < 1, the Neumann lemma (p 45 in Ortega and Rheinboldt (1970)) gives r=
(I - G')-' = I
(All) Application of the transformation P-l, where P is defined by eq 18, together with the use of eq A10 and A l l gives the equation A V = [I + A + Az + ...I R (AW
J and where
t k = (V2GkL+ V2Gk")-'(&" -
nkL)-
+ G' + (G'), + (G')3+ ...
(A3)
Tk
and the subscript k denotes a stage index. Thus, eq A1 accounts for the change in vapor component molar flows throughout the separator for perturbation in any stage temperature through eq A2 and A3. For any binary homogeneous mixture, the matrix V2GkL V2GkVis inverse nonnegative (see p 1768 in Lucia (1986)). Furthermore, in the absence of retrograde effects and supercritical components in the separator, the vector - BkLis a positive vector. Consequently for AT, > 0, the vector tk, which represents the change in the vapor component molar flow rates for stage k,is a positive vector. Also note that the corresponding change in liquid molar flow rates for stage k is - t k , which is a negative vector (see p 1766 in Lucia (1986)). It is more convenient to express the vectors tk and - t k in the form
where AV = (AV,, AV,, ...,AV,,JT,A is defined by eq 20, and the vector R is given by
+
nk"
tk = -tk
p&k
- vflk
= i k 2 k - LkXk
Approximation of each element of the matrix-vector product in eq A12 by the first nonzero term, which is also the dominant term, shows that
(-44) (A5)
where p k , i k , j j k , and are the total vapor flow rate, liquid flow rate, vapor composition, and liquid composition associated with stage k that result from the temperature perturbation AT,. For sufficiently small AT,, it follows that
(1 - a k - Z ) ( l - ak-1) AV =
Apk
!
1 + ak-1(1 - @k) - (1 - ak)(l - @ k + l ) - ak-l@k)(l- @ k + l ) - ak(1 - @ k + d 2 (1 - ak+l)bk+l(l - @k+2)
(A141
.
802 Ind. Eng. Chem. Res., Vol. 28, No. 6, 1989
Observe that for j = 2, 3, ..., ( k - 1)
Similar application of eq 11 and 12 to stage ns - 2 and use of eq B5 and B6 gives
k-1
avj/aTk
= A v k ( n , ( l - ai))> 0 I =I
because all avk/aTk
(A151
E (0,l). Furthermore, = A V k ( a k - i ( l - Pk) + a k + P k + l ( l - a h ) ) > 0
= [I - (1- Mns-~)M,B-iBns-i,ns-21-1(I =
M”s-2)A1n,-3
ai
(-416) since
G s - 2
(Yk-1, a k , P k
E (0,l). Finally, for j > k
- Mns-2)Alns-3
Ans-2,ns-3(I
039)
and Avns-2
=
+ Mns-lBns-l,ns-~Ans-~,ns-31 X (1 - Mns-2)AIns-3 = MnS-2Bns-2,ns-3Alns-3
MnS-2[1
Continuing in this manner, we obtain the expression Avk+2
since all a, /3 E (0,l). Consequently, eq A15-Al7 show that
= =
Mk+2Bk+2,k+lAlk+l Mk+PBk+P,k+lAk+l,k(I
- Mk+l)Alk
(B11)
However, Taylor series expansion of Alk in eq B11 shows that A l k = ( t k - L k ) X k + J!/k(kk - x k ) = h L k X k + L k ( k k - x k ) (BW where the quantities t k and k k are the total liquid flow rate and composition associated with stage k at Tk ATk. Use of eq 17 and subsequent substitution of eq B12 in eq B11 gives
+
The partial derivatives d L j / d T k can be analyzed in a similar manner and show that
=
Avk+?Yk+2
[GkPk+l
Lk(bk+l
+
- Pk+l)lMk+2Bk+2,k+lAk+l,kXk+l
where the quantity in brackets is a scalar quantity. The results in hppendix A show that for small perturbations ATk > 0, L k < L k . Moreover, it can be shown that P k + l > P k + l and that tkPk+l LkPk+l (B14
which establishes the condition given by eq 25.
Appendix B: Structure of the Jacobian Matrix for the ( T , P ) to ( Q ,P ) Mapping In this appendix, we show that
This follows from the fact that L k a 1 / T k and @k+l 0: T k for small ATk Consequently, the scalar quantity in eq B13 satisfies
+ Lk(Pk+l - Pk+l)
UkPk+l
;z
0
(B15)
and
aQk+i/dTk
=
- Hl+2(aVk+2/dTk)l
[Hl+l(aVk+l/aTk)
+ Hk+i(aLk+i/
[-HkL(aLk/aTk)
+
Application of eq 11 to stage k Appendix A show that Avk+lYk+l
=
AVk+i
- Lk(aHkL/aTk)
VTk)]
The last term in eq B3 is negative. We show that the terms in brackets are also negative. Consider then the first term in brackets and begin by comparing the quantities aVk+l/aTkand dVk+z/dTk. To do this, we use eq 11 and 12 recursively and traverse the column from stage ns to stage k + 1. For stage ns,we have Avns
=
(B4)
MnsAlns-1
For stage ns - 1 application of eq 11and 12, use of eq B4 and some algebra give the relationships Alns-1
=
Ans-l,ns-dI
- Mns-JAlnS-2
035)
and Avns-1
=
Mns-1Bns-l,ns-2ALns-2
(B6)
=
Mk+l[Alk
=
Mk+l[ukXk
= k‘ [
(B3)
-
(B16)
0
aVk+z/aTk
Taking partial derivatives of the energy balance equation for stage k + 1 with respect to Tk gives
+ 1 and the results in
4- A v k + , l
- (‘_kPk+l
+ LkAXkI +
Lk(Pk+l
- Pk+l))lYk+l
Since the term in parentheses in eq B17 goes to zero faster than u k as ATk 0, it follows that
/ a Tk I >> l a v k + 2 / a Tk I
l a v k +1
0318)
Moreover, if Hi+l and Hi+2are of the same magnitude, then eq B18 implies
IH’a;1(8v k + 1/
Tk)
I > I H 1 + 2 (av k + 2 / d T k ) I
(B19)
This together with eq 25 and 29 shows that the first term in brackets in eq B3 is negative. Similar arguments can be used to show that the second term in eq B3 is also negative and consequently that the inequality given in by eq B1 is valid. Finally, the inequality given by eq B2 can be established in exactly the same manner.
Literature Cited
where Ans-1,ns-2
= [I - (1 - Mns-i)Mnsl-l
037)
Acrivos, A.; Amundson, N. R. Application of Matrix Mathematics to Chemical Engineering Problems. Znd. Eng. Chem. 1955, 47, 1533-1541.
and Bns-l,nS-z
= 1 + MnsAns-1,ns-2(I
- Mns-J
(B8)
Boston, J. F.; Britt, H. I A Radically Different Formulation and Solution to the Single-Stage Flash Problem. Comput. Chem. Eng. 1978, 2, 109-121.
Ind. Eng. Chem. Res. 1989,28, 803-808 Broyden, C. G. A Class of Methods for Solving Nonlinear Simultaneous Equations. Math. Comput. 1965,19, 577-594. Crowe, C. M.; Nishio, M. Convergence Promotion in the Simulation of Chemical Process-The General Dominant Eigenvalue Method. AZChE J . 1975,21, 528-533. Dennis, J. E.; Schnabel, R. B. Least-Change Secant Updates for Quasi-Newton Methods. SIAM Rev. 1979,21, 443-459. Doherty, M. F.; Perkins, J. D. On the Dynamics of Distillation Processes. I V Uniqueness and Stability of the Steady State in Homogeneous Continuous Distillations. Chem. Eng. Sci. 1982,37, 381-392. Friday, J. R.; Smith, B. D. An Analysis of the Equilibrium Stage Separations Problem - Formulation and Convergence. AZChE J . 1964, IO, 698-706. King, C. J. Separation Processes; McGraw-Hill: New York, 1971. Lapidus, L.; Amundson, N. R. Stagewise Absorption and Extraction Equipment-Transient and Unsteady State Operation. Znd. Eng. Chem. 1950,42, 1071-1078. Lucia, A. Uniqueness of Solutions to Single-Stage Isobaric Flash Processes Involving Homogeneous Mixtures. AZChE J. 1986,32, 1761-1770. Michelsen, M. L. The Isothermal Flash Problem. Part 11. PhaseSplit Calculation. Fluid Phase Equilib. 1982, 9, 21-40. Orbach, 0.;Crowe, C. M. Convergence Promotion in the Simulation of Chemical Processes with Recycle-The Dominant Eigenvalue Method. Can. J . Chem. Eng. 1971,49, 509-513. Ortega, J. M.; Rheinboldt, W. C. Iterative Solution of Nonlinear
803
Equations in Seueral Variables; Academic Press: New York, 1970. Prausnitz, J. M.; Andersen, T. F.; Grens, E. A.; Eckert, C. A.; Hsieh, R.; O’Connell, J. P. Computer Calculations for Multicomponent Vapor-Liquid and Liquid-Liquid Equilibria; Prentice-Hall: Englewood Cliffs, NJ, 1980. Rosenbrock, H. H. A Theorem of Dynamic Conversion for Distillation. Trans. Znst. Chem. Eng. 1960, 38, 279-287. Rosenbrock, H. H. A Lyapunov Function with Applications to Some Nonlinear Physical Problems. Automatica 1962,1, 31-53. Sujata, A. D. Absorber-Stripper Calculations Made Easier. Hydrocarbon Process. Pet. Ref. 1961, 40, 137-140. Varga, R. S. Matrix Zteratiue Analysis; Prentice-Hall: Englewood Cliffs, NJ, 1962. Van Dongen, D. B.; Doherty, M. F.; Haight, J. R. Material Stability of Multicomponent Mixtures and the Multiplicity of Solutions to Phase Equilibrium Equation. 1. Nonreacting Mixtures. Znd. Eng. Chem. Fundam. 1983,22, 472-485. Wegstein, J. H. Accelerating Convergence of Iterative Process. ACM Comm. 1958, I , 9-13. Westerberg, A. W.; Hutchison, H. P.; Motard, R. L.; Winter, P. Process Flowsheeting; Cambridge University Press: Cambridge, England, 1979. Received for review July 29, 1988 Revised manuscript received February 11, 1989 Accepted February 17, 1989
Magnetic Filter for Solids: Theory and Experiment Shoichi Kimura and Octave Levenspiel* Chemical Engineering Department, Oregon State University, Corvallis, Oregon 97331
This paper introduces a new type of filter for removing ferromagnetic solids from a flow stream, one that freezes the particles to the tube wall by magnetic forces alone. Basic theory is developed t o tell how to find the parameters of any given system, and experiments are run to verify the theory. This is a report on a new type of filter that is designed to remove ferromagnetic solids from a slurry flowing in a tube or pipe. It acts by creating a strong magnetic field in the flow channel so as to freeze the solids to one side of the flow channel. Solids build up there, and eventually the flow channel becomes so restricted that no more solids can be captured. It is then time to turn off the magnetic field, which releases the captured solids to a collection chamber. The whole operation is then repeated. Asema (1988) presently offers such filters for removing suspended magnetic fines from industrial wastewaters. Another use for this device is for recovering ferromagnetic catalyst particles or catalyst particles having ferromagnetic cores from fluid streams leaving slurry reactors. At the beginning and close to the end of an operating cycle, a small fraction of the flowing solids may escape the magnetic filter. Hence, if all the solids are to be recovered, the magnetic filter will have to be followed with an ordinary filter. Nevertheless, with this magnetic filter, the duty of the mechanical filter will be greatly reduced because the magnetic filter can hold a large amount of solid with only a relatively small increase in pressure drop. In general, low pressure drop is a special advantage of the magnetic filter over its mechanic61 counterpart. This follows from the fact that it has no mechanical obstruction in its flow channel.
Theory Consider a pair of activated electromagnetic pole pieces snugly fitted to a paramagnetic tube through which is flowing a slurry containing ferromagnetic particles, as shown in Figure 1. With a high enough field intensity, 0888-5885/89/2628-0S03$01.50/0
particles can be arrested at the pole pieces, and since the magnetic flux density between the pole pieces is highest when the air gap is smallest, particles collect there first. Once a layer of particles is formed, the magnetic flux density increases because of the increased permeability of this solid layer. The solid layer grows progressively while reducing the cross-sectional area of the flow channel. This process continues until the magnetic forces just balance the net frictional and gravitational forces on the particles. This progression is shown in Figure 2. This section develops the basic theory which relates the capture current and the quantity of solids captured with the pole piece geometry, particle size, system geometry, and operating conditions. Balance of Forces. The pertinent forces influencing the action of a magnetic filter are the magnetic force, gravity, and hydrodynamic drag. Let us focus attention on a spherical particle that is touching a layer of already frozen particles. The net force tending to sweep the particle downstream is the result of gravity and drag, given by
where Cd is the drag coefficient and u is the mean velocity of the slurry. In eq 1,j = +1 for upflow of slurry, and j = -1 for downflow of slurry. Opposed to this is the vertical component of the magnetic force, F,, acting on the particle, or
FZ = CfF, 62 1989 American Chemical Society
(2)
