Ind. Eng. Chem. Process Des. Dev. 1986, 25, 184-188
184
L = liquid flow rate per unit weir width, m3/(s,m of weir width) Pe = Peclet number, dimensionless t = time, s U G = superficial air velocity, based on the total area in the column, m/s U L = mean linear velocity of the liquid across the plate, m/s Greek Letters ug2 = a2/r2,dimensionless T = mean liquid residence time, s C T ' = variance of the distribution curve about its mean, s2
Literature Cited Aleksandrov, I.A.; Sheynman, V. I.; Abroslmov, B. Z.; Svezhencev. V. S. USSR Patent 169061, 1965. Arafa, M. A,; Chekov. 0. S. Teor. Osn. Khim. Tekhnol. 1972, 6 , 343. Barker, P. E.; Self, M. F. Chem. Eng. Sci. 1962, 1 7 , 541. Bell, R . L. AIChE J. 1972, 18 (3),498. Berkovskiy, M. A.; Aleksandrov, I. A,; Skoblo, A. I. Chim. Techno/. Topliv Masei 1971, 16 (I), 32.
Blddulph. M. W. AIChEJ. 1977, 23(5),770. Danckwerts, P. V. Chem. €ng. S d . 1953, 2 , 1. Gautreux, M. F.; O'Connel, H. E. Chem. Eng. Progr. 1953, 57 (5), 232. Gerster, J. A.; Hill, A. B.; Hochgrab, N. N.; Robinson, D. G. Final Report of the University of Delaware. American Institute of Chemical Engineers, New York, 1958. Gilbert, T. J. Chem. Eng. Sci. 1959, 10 (4),243. Kafarov, V. V.; Shestopalov, V. V.; Gorenshteyn, 8. M. J. Prikl. Khim. 1969, 42 (2),368. Miganchi, T.; Imai. H.; Kubo. K.; Kogaku, K. Chem. Eng. (Jpn.) Abrid. Ed. Engl. 1967, 5 , 20. Porter, K. E.; Lockett, M. J.; Lim. C. T. Trans. Inst. Chem. Eng. 1972, 50,
91. Tassev, Zh. A.; Stefanov, Zh. S. Khim. Neft. Mashinostr. 1982, 4 , 35. Taylor, G. I.Proc. R . SOC.London. Ser. A 1953, A219, 189. Tichacek, L. J.; Barkelew, C. H.; Baron, T. AIChEJ. 1957, 3. 439. van der Laan, E. T. Chem. Eng. Sci. 1958, 7 , 187. Welch, N. E.; Durbin, L. D.; Holland, C. D. AIChE J. 1964, 10 (3), 373.
Received for review August 5, 1983 Revised manuscript received April 22, 1985 Accepted June 2, 1985
Simplified Flash Calculations for Cubic Equations of State Michael L. Micheisen Instituttet for Kemiteknik, Danmarks Tekniske
H 0jskole, Bygning 2 2 9 DK-2800 Lyngby, Denmark
An efficient flash procedure for systems described by the Cubic Equations of State with all binary interaction coefficients equal to zero is presented. I t is shown that the total set of material balances and equilibrium relations can be reduced to three nonlinear algebraic equations. Methods for efficient automatic generation of initial estimates are discussed, and a modification of the Newton-Raphson method for solving the set of equations is described. This modification is based on a constrained Gibbs energy minimization principle and effectively prevents divergence of the iterative method.
The Cubic Equations of State have found widespread engineering applications, in particular, for describing the phase equilibrium for light hydrocarbon mixtures. One reason for their popularity is that considering their simplicity, they yield a remarkably accurate description of the phase behavior. Parameters of the Cubic Equations are obtained primarily from pure component properties, but in order to increase their accuracy and applicability, it is customary to introduce binary interaction coefficients, kij,which are obtained by fitting the Equation of State to experimental data for binaries. In many cases, the resulting gain in accuracy is modest, and if binary data are missing, taking kiJ = 0 is usually a good approximation. Another incentive for neglecting the binary interaction coefficients is that this permits very significant simplifications in the calculation of phase equilibrium for multicomponent mixtures, in particular when the number of components is high. The present paper describes a new procedure for isothermal flash calculations with a Cubic Equation of State, where all binary interaction coefficients equal zero. The methodology and the computational gains with the new method correspond to those of a similar procedure for critical point calculations, described earlier (Michelsen and Heidemann, 1981). A significant drawback of the proposed method is that it is restricted to particular Equations of State and requires suboptimal parameters for these. One application will therefore be for preliminary design calculations, where the inherent errors of the thermodynamic model are acceptable. A different situation, where the method could be 0196-4305/86/1125-0184$01.50/0
useful, is for modeling mixtures containing a large number of not well-defined components,i.e., in the characterization of heavy oils (Pedersen et al., 1984a, 1984b). For such mixtures, the requirement for binary interaction coefficients equal to zero can be compensated for through judicious selection of the critical parameters for the pseudocomponents chosen in the characterization. Cubic Equations of State We consider cubic two-constant equations of state of the form a p = - -RT (1) u - b (U + 61b) (U 62b)
+
where and b2 are numerical constants and where the mixture parameters a and b are given by the mixing rules a
CC~i~jaij t
(2)
J
and b = CyJb,
(3)
J
where y is the vector of mixture mole fractions. Usually, a, = (aiLall)1/2(l - kJ, where the binary interaction coefficients k,, are small, typically zero for hydrocarbon-hydrocarbon interactions and in the range 0-0.15 otherwise. The pure component parameters a,, and b, are given by a,, =
fiAR2Tct2/Pct(1 + m,(l - (T/Tc2)1'2)}2 (4)
with 0 1985 American Chemical Society
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 1, 1986 185
mi = m ( q )
(5)
and bi = QBRTci/Pci
mol and L = (1 - P)F mol of composition y and x, respectively, the material balances and equilibrium relations are
(6)
2;
+ (1- P)x;
= pyi
Tci,PCi,and ai are the critical temperature, critical pressure, and acentric factor for component i, and the constants QA and QB are determined such that the correct critical temperature and pressure are obtained for a pure component. For two frequently used Equations of State we have Soave-Redlich-Kwong equation (Soave, 1972) 61
= 1, 62 = 0, QA = 0.427 48, ClB = 0.086 64, m(wi) = 0.480
+ 1 . 5 7 4 ~- ~0.176wi2
x.4.' 1 1 = y.4.v 1 1
= 2'1'
+ 1, 62 = -2'1' + 1, ClA = 0.45724, QB =
(15)
This set of (2N + 1)equations in the 2N + 1 unknown (x,y,P)can be reduced to three equations as follows: Multiplication of the material balances, eq 14, by a,and summation yields Cziai = PCyia, + (1 - P)Cxiai i
1
(17)
1
or
0.077 80
+ 1.542 2 6 q - 0.269 92wi2
m(q) = 0.374 64
(14)
and
Peng-Robinson equation (Peng and Robinson, 1976) 61
i = 1,2,...,N
= Pa,
+ (1 - P)al
(18)
bF = Pb,
+ (1 - P)bi
(19)
"F
and similarly
Fugacity coefficients are given by
These relations imply that given the V-phase mixture parameters a, and b,, the corresponding L-phase parameters are given by a' =
with a/ = Caijyj j
(7)
If it is assumed that all interaction coefficients kij are equal to zero, we obtain aij = (a..a..)1/2 11 11 = a.a. 1 I' with ai = ai?/z
(8)
a = C C y a j a i j = a2
(9)
a = Cyjaj
(10)
and i
l
where I
Finally, a/ = Caijyj = a p 1
(11)
"F - Pa, ~
1-P
, biz-
bbu 1-0
bF -
(20)
When Ki = 4/14? is defined, combination of eq 14 and 15 yields the well-known relation xi =
zi
1 1 + @(Ki- 1)'
A computational procedure involving only the three independent variables a,,b,, and /3 can now be formulated as follows: (i) At the given T and P, evaluate the ai and bi and calculate aF and bp Estimate values of a,, b,, and /3. (ii) Calculate ai and bl from eq 20. Evaluate from eq 12 and 13 the 4/ and 4; and the K factors, Ki = dil/@;". (iii) Calculate V- and L-phase mole fractions y and x from eq 21. (iv) Evaluate check functions el =
C(yi - xi) i
(224
When these relations are substituted into eq 7, the fugacity coefficient expression can be rewritten In
4i = 40 + qaai + qpbi
(12)
where qo, qa and qg depend on a,b, T, and P only (taking into account that the volume v is implicitly given by a,b, T , and P). The expressions for the coefficients are
Equation for Flash Calculation For F mol of overall composition z at a given temperature and pressure split into two phases containing V = PF
If these are all equal to zero, the phase distribution a t equilibrium has been determined. If not, perform an iterative correction of a,,b,, and (3 and return to step (ii). To converge eq 22, we recommend the Newton-Raphson iteration using analytical derivatives. The relation between the independent variables and the check functions appears to be complex, but analytic derivatives can be coded compactly and efficiently.
Initial Estimates Frequently good initial estimates of a,, b,, and 0 are available from a previous phase split calculation under conditions similar to the actual case. Naturally such information, if available, should be used. In other cases, it may not even be known in advance whether a given specification corresponds to two equilibrium phases. Then it is recommended to generate initial
186
Ind. Eng Chem Process Des. Dev., Vol. 25, No. 1, 1986
estimates by means of stability analysis, as suggested by Michelsen (1982). Stability analysis for a mixture of given composition z involves determining all solutions to the equation In Y,
+ In 4,(Y) - In 2, - In 4,(z) = 0
(23)
where the Y, can formally be interpreted as the mole numbers. The mixture is stable a t the specified T and P provided C,Y, 5 1 a t all solutions of eq 23. Stability of course implies that two equilibrium phases do not exist. If solutions of eq 23 with C,Y, > 1 are found, the mixture is unstable, and the Y, values can be used to generate initial estimates for a subsequent phase split calculation as described below. The set of eq 23 can be reduced to two equations involving only (a*$*),the CY and b parameters of the Y phase, as follows: (i) Make an initial estimate of CY*and b*. Evaluate the fugacity coefficients &(z) for the feed. (ii) Calculate the fugacity coefficients of the Y phase by using eq 12 and 13. Evaluate the phase compositions Y, from In Y, = In z,
+ In $,(z) - In +,(a*,b*)
(24)
Unfortunately, cases are also encountered where divergence occurs or where an excessive number of iterations is required. Step-size reduction is of little help, whereas modification of the computational procedure effectively eliminates the convergence problem. The modification utilizes a switch to solving a simpler problem whenever the Newton-Raphson method predicts an excessive correction. The problem simplification consists of fixing the value of @ and solving for the phase distribution which at fixed @ yields the minimum Gibbs energy of the system. Apparently constrained Gibbs energy minimization of this type has not been described earlier, and hence a fairly complete development will be presented. Given F mol of feed of composition z at specified T and P, we wish to determine the V-phase mole numbers (ul,uz,.....,uN) subject to the constraint V = C,u, = PF, such that the system Gibbs energy G is minimized. The minimum is found as the stationary point of the Lagrangian function Q (Fletcher, 1981)
Q
CYI((Y,- C Y * )
that is,
aQ/auI = aG/au, - A = 0, i = 1,2,.....,N
(29)
aQ/ax = - ( ~ u -, PF) = o
(25a)
1
g, = CY,(b,- b*)
(28)
1
(iii) Evaluate the check functions g, =
= G - A(Cu, - P F )
1
(25b)
Substituting dG/au, = p l u - p l i , we obtain at the solution
I
If g, and g2both equal zero, a solution of eq 23 has been determined. Otherwise iterative correction is required. For this purpose, the Newton-Raphson method with analytic derivatives is recommended. In order to verify stability (or instability), it is in practice sufficient to solve eq 23 twice, using as initial estimates (a*$*)= (al,bl)and (a*$*)= (aN,bN), respectively, where indexes 1 and N refer to the lightest and the heaviest component in the mixture. The outcome of the stability analysis can be (1)convergence to the “trivial solution” (a* = a&* = bF), (2) convergence to a nontrivial solution with C,Y, I 1, or (3) convergence to a nontrivial solution with C,Y, > 1. If the latter outcome is obtained for one or both of the initial estimates, the mixture is unstable. Converging the Phase Split Calculation. Provided that the stability analysis leads to Y vectors with C,Y,> 1, we use as initial estimates for the K factors for the phase split calculation
KL = Y J z ,
(26)
and subsequently solve the Rachford-Rice equation
Cz,(K,- 1)/(1+ P(K, - 1))= 0
(27)
1
for 6. The mole fractions of the individual phases are found from eq 21, and P and (CYJJJ calculated from composition y are used as initial estimates for the solution of eq 22. In the case where both initial estimates for the stability analysis converge to a nontrivial solution with C,Y, > 1, phase compositions using eq 26 and 27 are generated for both, and the set of compositions with the lower Gibbs energy is selected for continued iteration. With this choice of initial estimates, convergence (norm of the correction vector less than lo4) normally takes place in two to five iterations. Near ideal behavior or cases with one equilibrium phase present in small amounts lead to particularly rapid convergence.
,LIIV
(30)
- p1L = A
or
In y l + In 4,” - In x , - In
@lu
= A/(RT) = A’
(31)
Using as earlier K , = @11/41u, eq 31 becomes y , = BK,x,, with 6 = exp(A’)
(32)
The modified “phase equilibrium” problem differs from the original, eq 14-16, only in the modification of the equilibrium relation. Thus, the procedure leading to eq 22 can still be used provided we replace K , by OK, in the calculation of the mole fractions (x,,yJ. Furthermore, /3 is now a parameter of the problem, and A’ (or O ) is the new independent variable. The solution a t different values of yields A‘ in dependence of @,A’ = A’@), and the solution to the original phase equilibrium problem is the particular solution to the modified problem, a t which A‘ = 0. The theory of constrained minimization (Fletcher, 1981) also shows that at the solution of eq 29, aG/ap = FA or l/(FRT)dG/dA = A’ (33) Hence, in order to obtain the global (unconstrained) minimum of the Gibbs energy, @ must be increased when A‘ is negative and decreased when A’ is positive. Furthermore the Newton-Raphson determination of a”, b,, and A’ permits an inexpensive simultaneous determination of the derivatives of these quantities with respect to 0.These derivatives are of aid in selecting new P values and in extrapolating to the new conditions. One point of the A’(@) curve is available from the stability analysis: It is readily observed that by taking (3 = 0, Le., x , = z,, eq 31 and eq 23 are equivalent, with A’ = - In C,Y,. Therefore, calculation of A’@) could be initiated with = 0, A‘ = - In C,Y,,using a tracing technique as described in Michelsen (1980). Since we are only interested in the particular point on the curve, at which p’ = 0, the following procedure is preferable:
Ind. Eng. Chem. Process Des, Dev., Vol. 25, No. 1, 1986 187 Table I. Iteration Course for 22-Component Mixture with Initial Estimate of Lower Gibbs Energy Newton step t z actual step iter. no. P 1 0.804 -0.32 -0.15 ( t J -0.061 ( t z ) 2 0.654 -0.061 3 4
0.005 0.0
0.593 0.598
0.005 ( t z ) 0.0 ( t 2 )
Table 11. Iteration Course for 22-Component Mixture with Alternative Initial Estimate Newton step t 2 actual step a iter. no. 1 0.005 0.005 ( t 4 ) 0.005 0.013 2 0.011 ( t 4 ) 0.010 0.027 ( t 4 ) 0.032 3 0.021 4 0.056 ( t l ) 0.048 0.41 5 0.108 ( t 4 ) 0.104 0.096 (ti) 6 0.212 -0.29 7 0.150 ( t J 0.308 8 9
10
In eq 21, we replace Ki by exp(X’)Ki and evaluate (e = (el,e2,e3)T.In addition, evaluate the vectors of partial derivatives hl, h2,h3,and h4,where hl = de/&,, h2 = &/ab,, h3 = de/dp, h, = de/dA’ Linearization of eq 22 yields (34)
This set of three equations with four variables is treated as follows: The general solution to eq 3 can be written
where
rlhl + r2h2+ r4h4+ e = 0
(36)
+ s2h2 + ~4h4+ h3 = 0
(37)
and
slhl
With this’choice,eq 34 is satisfied for any value of t, and the actual value used is chosen according to the following criteria: (i) Calculate the value of t , t = t,, which minimizes
SSQ(tl)= ( A c ~ , / a , ) ~+ (Ab,/b,)’
+ AP2 + AXf2
-0.010 0.0
0.150
(t3)
-0.010 ( t z ) 0.0 ( t z )
+
Figure 1. A’(@) for example 22-component mixture.
hlAa, + h2Ab, + h3A@+ h4AA’+ e = 0
0.30
0.458 0.608 0.598
(38)
If SSQ(tl)is greater than lo9, the current iterate is far from the A’(@) curve, and refinement is required before proceeding. Set t = tl, and continue a t step (iv). (ii) When SSQ(tl)is small, the current iterate is close to the curve, and we can proceed toward the solution. Taking t = 0 corresponds to no change in p and predicts a value of A‘, Apred = A‘ + r4. From eq 33, we note that @ must be increased ( t > 0) when Xpred is negative and increased ( t < 0) when Aprd is positive. The ultimate objective is to arrive a t a solution with A’ = 0. We calculate t 2 = -r4/s4(the Newton step) and take t = t2when t 2 < A@,, and t2hA < 0, and t = t3 = -A@” sin (A’prd) otherwise, where A@- is a preselected maximal step length, e.g., 0.15. (iii) In rare cases, this selection leads to excessive changes in a, or to /3 values exceeding the bounds 0 < @ < 1. Hence, t is modified (reduced in magnitude) in the following
situations: (a) If (a, Aa, - a f ) / ( a u - at) < 0.65, choose t = t 4 where (a,- Aa&) - a f ) / ( a u- at)= 0.65. (b) If @ t < @/2,choose t = t 5 = -@/2. (c) If p + t > ( p 1)/2, choose t = t 6 = (1- @)/2. (iv) Calculate the correction vector. If convergence has not occurred, reevaluate the e and h vectors, solve for r and s, and return to step (i).
+
+
Example Problem In Figure 1is shown A’(@) (at T = 371 K, P = 393 atm) for a 22-component hydrocarbon mixture including a range of heavy components up to c65. The characterization procedure of Pedersen et al. (1984a) was used for determining critical properties for the heavy fractions. Both trial phases for the stability calculation reveal instability, the light component trial phase with C,Y,- 1 = 8.8 X and the heavy trial phase with CiYL- 1 = 5.6 X The corresponding initial phase splits yield decreases in the Gibbs free energy ( G / ( R T ) of ) -6.6 X 10“ and -9.6 X respectively. The iteration sequences using both initial estimates are shown in Tables I and 11. With the initial estimate based on the light trial phase, only four iterations are required, and the Newton step is only modified after the first iteration. The initial estimate from the heavy trial phase requires a total of 11additional iterations. Small steps are taken near @ = 0, where A’(@) is extremely steep, and actually the step length is limited by the requirement that excessive changes in a, are not permitted. In the intermediate region, where A’ decreases with increasing @,the Newton steps are either far too large of in the wrong direction, but the modifications enable the algorithm to arrive safely a t the correct solution. Although convergence is obtainable with both initial estimates, the advantage of selecting the estimate with the lower Gibbs energy is obvious in the present example. Indeed, the majority of problems do not require the auxiliary variable A’ and the precautions for step-length reduction, but these modifications of the simple NewtonRaphson method provide additional safety for difficult problems without sacrificing speed in easier problems. The computing time for the above example is about 0.024 s for the sequence of Table I (IBM 3081, FORTRAN H compiler). A sequence of calculations where the pressure is decreased in steps of 20 atm and the results of the previous step are used as an initial estimate requires only on the average 0.006 s/calculation. The time consumption is only weakly affected by the number of components
188
I&. Eng. Chem. Process Des. Dev. 1986,25, 188-197
present. A similar example with 44 components results in an increase of about 50%. Conclusion and Suggestions for Further Work A rapid and safe algorithm for isothermal flash calculations using the Cubic Equations of State, where all binary interaction coefficients are equal to zero, has been developed. The procedure is particularly attractive in situations, where the advantage of a drastic reduction in computer time expenditure is of much higher importance than the disadvantage of a slight decrease in accuracy, e.g., in reservoir simulations. For mixtures where the assumption of binary interaction coefficients equal to zero is unacceptable, e.g., in hydrocarbon mixtures containing significant amounts of carbon dioxide, similar simplifications may still be possible. It can be shown that each component having non-zero interaction coefficients with the remaining components introduces two new equations. Hence, phase equilibrium in hydrocarbon-carbon dioxide mixtures can be described by a total of five equations, irrespective of the number of hydrocarbons, provided only the hydrocarbon-carbon dioxide interaction coefficients are non-zero. Adaptation of the present procedure to such mixtures is currently being investigated.
Nomenclature mixture parameter of the Cubic Equation of State a pure component parameter, eq 4 aii mixture parameter, eq 2 ai; a/ defined in eq 7 mixture parameter of the Cubic Equation of State b pure component parameter, eq 6 bi discrepancy functions, eq 22 moles of feed discrepancy functions for stability analysis, eq 25 Gibbs energy vectors of partial derivatives of e component indexes binary interaction coefficient equilibrium factor, component i, K, = y , / x , moles of L phase number of components in mixture N P pressure, atm
Te3
P, q,,, q,,
Q qs
xr4 sI-s~
t t,-t6 T Tc V
4
V xr YI
YL 2,
critical pressure, atm constants in fugacity coefficient expression, eq 12 Lagrangian function solution to eq 36 gas constant solution to eq 37 step length step length from different criteria temperature, K critical temperature, K molar volume moles of component i in V phase total moles of V phase mole fraction of component i in L phase mole fraction of component i in V phase trial-phase mole number of component i mole fraction of component i in feed
Greek Letters modified a parameter of the Cubic Equation of State, eq 10 fraction of V phase, p = V / F P constants of the Equation of State 6,, 6, 0 defined in eq 33 x Lagrange multiplier a
V@T)
A' $1
PI w
OA,
QB
fugacity coefficient, component i chemical potential, component i acentric factor constants of the Equation of State
Literature Cited Fletcher, R. "Practical Methods of Optimization"; Wiley: New York, 1981; Vol. 2, Chapter 9. Michelsen, M. L. Fiuid Phase Epuilib. 1980, 4 , 1-10, Michelsen, M. L.; Heidemann, R. A. AIChE J. 1981, 27 (3), 521-523. Michelsen, M. L. FiuM Phase Epuiiib. 1982, 8 , 1 ~ 1 9 . Pedersen, K. S.; Thomassen, P.; Fredenslund, Aa. Ind. Eng Chem. Process Des. Dev. 1984a, 3 3 , 163-170. Pedersen, K. S.; Thomassen, P.; Fredenslund, Aa. Ind. Eng. Chem. Process D e s . D e v . l 9 8 4 b , 3 3 , 566-573. Peng, D.-Y.; Robinson, D. B. Ind. Eng. Chem. Fundam. 1976, 15, 59-64. Soave, G. Chem. Eng. Sci. 1972, 27, 1197-1203.
Received f o r review October 16, 1984 Accepted June 17, 1985
Flow Model for the Solid in a Continuous Fluidized Bed with Increase of the Cross Section in I t s Upper Zone Jos6 Corella, Rafael Bllbao, Antonio Monzen, Javler LezaOn, and Fernando Fernandez Department of Chemical Engineering, Faculty of Science, University of Zaragoza, Zaragoza, Spain
A model for the flow of the solid in a continuous and improved fluidized bed has been developed and experimentally tested. Four different flukitzed beds with an enlargement in the cross section in the upper zone of the bed have been used. The model developed and used is a combined or compartment one, and it consists of three zones: a piston flow with dispersion, a backmixing flow in series, and a dead zone. The three parameters of this model are calculated from the RTD curves experimentally obtained in the four different contactors under different operation conditions. A small increase of the gas velocity at the bed surface above the uglumtof these modified fluidized beds gives rise to a considerable increase of the dispersion number of the solid flow. An increase of the fraction of the solid in the upper zone of the contactor produces an increase of sdid circulating in piston flow. The positioning of internals in the tronconical zone of the contactor produces a diminution of dead zones and of the dispersion number.
An enlargement of the cross section in the upper zone of a fluidized bed results in a reduction of the average superficial gas velocity in that zone. So, a highly fluidized
lower zone may be followed by a less fluidized and even by a fixed upper zone. This enlargement generates a new type of solid-gas contactor which can operate in very
0196-4305/86/1125-0188$01.50/00 1985 American Chemical Society
