Simulation of Multicomponent Multistage Vapor− Liquid Separations

07730 Distrito Federal, Me´xico. It is shown that the well-known and extensively used bubble-point methods, such as the one proposed by Wang and Henk...
0 downloads 0 Views 300KB Size
Ind. Eng. Chem. Res. 2003, 42, 175-182

175

Simulation of Multicomponent Multistage Vapor-Liquid Separations. An Improved Algorithm Using the Wang-Henke Tridiagonal Matrix Method Rosendo Monroy-Loperena† Ingenierı´a Molecular, Instituto Mexicano del Petro´ leo, Eje Central La´ zaro Ca´ rdenas 152, 07730 Distrito Federal, Me´ xico

It is shown that the well-known and extensively used bubble-point methods, such as the one proposed by Wang and Henke for simulating vapor-liquid multistage distillation problems, can be used efficiently for the simulation of other types of vapor-liquid multistage separation problems, such as absorption and reboiled-absorption processes, even when the thermodynamic properties are a strong function of the composition. Via numerical simulations, it is shown that a robust bubble-point temperature iteration scheme provides convergence stability in the whole bubble-point method. The proposed robust bubble-point temperature iteration scheme is a combination of an equation decoupling method and a simultaneous method. This hybrid bubblepoint temperature iteration scheme exploits the reliability of the equation decoupling method, which makes good convergence progress from a poor starting point. However, if the initial point in the bubble-point temperature calculation is good enough, only the equation decoupling method is used, to avoid excessive successive substitution steps when the separation factors are updated. As a consequence of the extension of the bubble-point method to simulate absorption cases, a didiagonal matrix with an off column is generated, when the interstage vapor rates are calculated. An efficient and simple algorithm derived from the Gaussian elimination method to solve the didiagonal matrix with an off column is presented. Results are reported for distillation, absorption, and reboiled-absorption processes, using the Soave-Redlich-Kwong equation of state. 1. Introduction

S Equations (summations of mole fractions)

The solution of a multistage separation problem requires mass balances, heat balances, and equilibrium conditions to be fulfilled by every stage of the process. The basic equations are often simple, but the countercurrent operation interconnects all stages. Amundson and Pontinen1 described the problem as the solution of a set of nonlinear simultaneous equations. For a system with N stages and NC components, the mathematical expressions that describe the countercurrent, multistage separation processes are derived by making material and heat balances around the jth tray of the model (see Figure 1). The stages are assumed to be in theoretical equilibrium, that is:

M Equations (material balance) Mi,j ) Vj+1yi,j+1 + Lj-1xi,j-1 - (Vj + SVj)yi,j - (Lj + SLj)xi,j + Fj zij ) 0 (1)

E Equations (equilibrium relationships) Ei,j ) yi,j - Ki,jxi,j ) 0

(2)

† E-mail: [email protected]. Tel: +52-(55)30038105. Fax: +52-(55)3003-6239.

NC

Sxj )

xi,j - 1 ) 0 ∑ i)1

(3)

NC

Syj )

yi,j - 1 ) 0 ∑ i)1

(4)

H Equations (heat balance) Hj ) Vj+1Hj+1 + Lj-1hj-1 - (Vj + SVj)Hj - (Lj + SLj)hj + FjHFj + Qj ) 0 (5) These equations are conveniently designated as the MESH equations,2 where xi,j, yi,j, and zi,j are the mole fractions of component i on stage j, respectively, in the liquid phase, vapor phase, and feed stream; Ki,j is the separation factor of component i on stage j, which is defined as yi,j/xi,j; Lj, Vj, SLj, and SVj are liquid, vapor, liquid side stream, and vapor side stream flow rates leaving stage j; and hj, Hj, and HFj are the liquid, vapor, and feed stream enthalpies, respectively. The separation factors and the enthalpy values depend on the stage temperature, Tj, liquid and/or vapor composition, and stage pressure, pj. The (large) set of N(2NC + 3) equations (that is, eqs 1-5) should be solved for the N(2NC + 3) unknowns (that is, Lj, Vj, Tj, xi,j, and yi,j). Current practice is based mainly on three strategies, namely, simultaneous correction methods3,4 (SC), 2N methods,5-8 and equation

10.1021/ie0108898 CCC: $25.00 © 2003 American Chemical Society Published on Web 12/10/2002

176

Ind. Eng. Chem. Res., Vol. 42, No. 1, 2003

The observations we report in this work are the following: (1) A BP method can handle vapor-liquid separations where the separation factors are strongly functions of the composition, the case usually found when the thermodynamic properties were calculated with an equation of state. This fact is possible through the use of a robust bubble-point temperature iteration scheme. (2) Furthermore, a BP method can handle wideboiling mixtures and may be a good candidate to solve absorption and reboiled-absorption problems without losing its stability. For the absorption case, an easy modification of the procedure to calculate the vapor flow rates is presented. The solution procedure, suggested by Friday and Smith5 and developed in detail by Wang and Henke,2 is used in this work. However, others BP methods can be used (see, for instance, work by King10). The improved BP method is applied in three practical cases, that is, distillation, absorption, and reboiled-absorption, using the Soave11 modification of the Redlich-Kwong equation of state. 2. General Strategy of Mathematical Solution of the Bubble-Point Method

Figure 1. Schematic representation of a general countercurrent cascade of N stages.

decoupling methods1,2,5,9 (ED). The criteria to select which method to use usually is based in the thermodynamic method to be applied rather than the numerical methods involved. That is, for vapor-liquid separations involving strongly nonideal mixtures, the SC approach appears to combine stability and computational speed in the best way. For vapor-liquid separations with ideal or only mildly nonideal solutions, the 2N approach is effective for handling feed mixtures with any sort of volatility characteristics. For vapor-liquid separations with ideal solutions, where the separation factor values are weakly functions of the composition, the ED arrangement can further accelerate the computation. Using the fact that the computing time for even the most general of the methods is small for one or a few solutions, computing time may become excessive when a very large number of solutions is to be made, as in optimization of the design of large chemical processes or dynamical studies. In these cases ED methods can be appropriate if stability is provided in these types of methods. On the other hand, ED methods are in agreement with the chemical engineers’ way of thinking because the model is divided into submodels relating the distinct physical processes and linked together according to a flowsheet topology. However, these procedures are reliable, easy to assemble, and usually robust.

In the BP methods, the MESH equations are grouped by types. With this approach, stage temperatures and flow rates are guessed. The M equations, eq 1, are combined with the E equations, eq 2, to form the first subset of equations. These equations are linearized by holding the values of the flow rates and separation factors invariant and are then solved componentwise for stage composition, using the Thomas algorithm (see, for instance, work by Chapra and Canale12). Using the calculated composition, the S equations, eq 4, and the H equations, eq 5, are then solved separately for the new values of tray temperatures and flow rates. The entire procedure is repeated until all equations are satisfied (cf. the works of Wang and Henke2 and Seader and Henley13). This computational procedure may be summarized in four steps, as given in Figure 2. Referring to Figure 2, note that steps 1, 2, and 4 are carried out in an explicit way to get the calculated variables (that is, Vj, Lj, and xi,j), and actions to delimit the variables between physical meaning values are straightforward. However, the functions used to obtain the tray temperature and the vapor compositions (step 3) are of implicit nature, and a numerical strategy is needed to find them. Therefore, a robust bubble-point temperature iteration scheme is needed to provide stability to the whole BP method. In the following subsections, the construction of a robust iteration scheme will be presented to obtain the temperature in each stage (solution of the E equations, eq 2, and the S equations for the vapor phase, eq 4), and some issues about the procedure used to calculate the vapor rates (solution of the H equations, eq 5, and the S equations for the liquid phase, eq 3) will be presented and discussed. 2.1. Bubble-Point Temperature Calculation. As mentioned, the temperature in each stage is calculated in the BP method by means of a bubble-point temperature calculation. The necessary bubble-point temperature equations are obtained from the E equations, eq

Ind. Eng. Chem. Res., Vol. 42, No. 1, 2003 177

Figure 2. Information flow block for a bubble-point method.

2, and the S equations for the vapor phase, eq 4. Thus, for each stage, we have NC

ψj(Tj) ) 1 -

∑ i)1

Ki,jxi,j

(6)

where ψj(Tj) represents a bubble-point temperature function in stage j to obtain the stage temperature, Tj. The stage temperature is updated iteratively assuming that the separation factors, Ki,j, are constant. Then the vapor mole fractions are evaluated from

yi,j ) Ki,jxi,j

(7)

The separation factors are updated between iterations in a process that could be written as k+1 V L k ) ln Ki,jk - ln(fˆi,j /fˆi,j) ln Ki,j

(8)

V L and ˆfi,j are the where k is the iteration index and ˆfi,j fugacities of component i in the vapor and liquid phases, respectively. Convergence is obtained when successive separation factors, Ki,j, are sufficiently close to each other. In general, the process of solving simultaneous equations through successive substitutions can be represented by the expression

t

k+1

k

) f(t )

(9)

which is presumed to have solution t* that satisfies

t* ) f(t*)

(10)

The behavior of the process depends on the eigenvalues of the matrix

E ) [∂fi/∂tj]t)t*

(11)

Convergence to t* is only possible if all of the eigenvalues of E are less than unity in absolute value (see, for instance, work by Ortega and Rheinboldt14). As shown by Michelsen,15 the successive iterative process becomes very slow near a critical point, because two eigenvalues of the relevant E matrix tend toward +1 at such points. One relatively simple procedure to overcome oscillatory nonconvergence of successive substitution in general is to apply a “damping factor” to the computed change in the separation factors between iterations (see,

for instance, work by Heidemann and Michelsen16). The modified process could take the form k+1 k V L k ln Ki,j ) ln Ki,j - R ln(fˆi,j /fˆi,j)

(12)

where R is the damping factor. The iteration scheme composed of eqs 6, 7, and 12 will be identified as a successive substitution bubble-point temperature (SSBPT) iteration scheme. As mentioned, the efficiency of successive substitution decreases when eigenvalues are close to unity, where overshooting makes frequent use of successive substitution steps necessary. Methods that directly take into account their composition derivatives of the separation factors are preferable under these conditions. That is, an iterative sequence of a global Newton-type approach, for solving NC + 1 equations formed by the E equations, eq 3, and the S equation for the vapor phase, eq 4, with NC + 1 unknowns, formed by the stage temperature, Tj, and the vapor mole fractions, yi,j (i ) 1, ..., NC). These equations, expressed as a function set, are

gi ) yi,j - Ki,jxi,j

for i ) 1, ..., NC NC

gNC+1 ) 1 -

yi,j ∑ i)1

(13)

Then, the system of NC + 1 equations (eq 13) is solved using a Newton-type iterative procedure. This iteration scheme will be identified as a simultaneous correction bubble-point temperature (SCBPT) iteration scheme. The SSBPT and SCBPT for solving the system of nonlinear algebraic equations for the bubble-point temperature have certain limitations. They have also certain desirable features. The SSBPT, in addition to its simplicity, is also a stable method.16 The SSBPT, unlike SCBPT, does not require good initial estimates of the separation factors. The classical Newton method has quadratic convergence properties. However, because of the overshoot, the SCBPT may fail to converge when the initial estimate is not a good one. In Newton’s method, various approaches have been suggested to adjust the step length and then to alleviate the overshoot problem. These adjustments do not guarantee convergence of Newton’s method in the early iterations. The combination of the SSBPT and SCBPT is a good choice and has the desirable features of both. In this approach, the SSBPT comprises the first few iterations, and later, when a switching criterion is met, SCBPT is used.

178

Ind. Eng. Chem. Res., Vol. 42, No. 1, 2003

In this work, we suggest to begin with the SSBPT. In the SSBPT to update the stage temperature, we use Muller’s method (see, for instance, work by Chapra and Canale12), as suggested by Wang and Henke,2 because of its reliability and lack of need of derivatives, making the search of the stage temperature considerably quicker (bubble-point temperature). If, in the course of the calculations, when the separation factors are updating by successive substitutions, this process takes more than 10 iterations, we switch to the SCBPT, taking the last values of the stage temperature and vapor mole fractions of the SSBPT as the initial values for the SCBPT. This hybrid scheme exploits the reliability of the SSBPT to find the stage temperature (bubble-point temperature) starting from a poor point that could be a common situation in the first iterations of a BP method. Note that, if the starting point is good enough, the SCBPT will not be used. 2.2. Vapor Rate Calculation. In the BP methods, to calculate the vapor rates, the liquid rates are expressed as functions of the vapor rates by an overall material balance of all stages, from the top through the jth stage (see Figure 1) j

Lj ) Vj+1 +

∑ (Fm - SVm - SLm) - V1

(14)

where

R1 ) H2 - h1 µ1 ) h 1 - H 1 γ1 ) SV1(H1 - h1) + F1(h1 - HF1 ) - Q1

{

βj ) Hj+1 - hj

Rj ) hj-1 - Hj µj ) hj - hj-1 γj ) Fj(hj - HFj ) + SVj(Hj - hj) - Qj + j-1

(hj - hj-1)

SLm) + Fj(hj -

HFj )

∑ (Fm - SVm m)1

+ SVj(Hj - hj) + V1(hj-1 - hj) Qj}/(Hj+1 - hj) (15)

In distillation problems, when the distillate rate and the reflux rate are specified (a common case used in simulation), V1 (vapor distillate) and V2 are fixed, so eq 15 is only used to get V3 to VN. In reboiled-absorption problems it is usual to specify V1, and therefore eq 15 is used to get V2 to VN. However, in absorption problems V1 is unknown and it needs to be calculated, so we cannot use eq 15 directly, but we can cast a linear system of equations as follows:

[

β1 R2 β2 R3 β3

l Rj βj

l RN-2

]

µ1 µ2 µ3 l µj l βN-2 µN-2 RN-1 βN-1 µN-1 RN µN V2 γ1 V3 γ2 V4 γ3 l l Vj-1 ) γj l l VN-1 γN-2 VN γN-1 V1 γN

[ ][ ]

(16)

}

The solution of system (16) for {V} can be easily obtained by use of a simple algorithm derived from the Gaussian elimination method. In this algorithm, two auxiliary quantities, γ*j and µ*j , are calculated by first evaluating γ1 and µ1 and advancing forward with j increasing, that is

γ*1 ) γ1

Equation 14 is used twice in the H equation, eq 5, to eliminate Lj-1 and Lj. After arrangement j-1

∑ (Fm - SVm - SLm) m)1

j ) 2, ..., N

m)1

Vj+1 ) {Vj(Hj - hj-1) + (hj - hj-1)

j ) 2, ..., N - 1

{

µ*1 ) µ1

* γ*j ) γj - Rjγj-1 /βj-1 * * µj ) µj - Rjµj-1/βb-1

}

j ) 2, 3, ..., N

(17)

Then, values of Vj are calculated by first evaluating V1. Thus

V1 ) γ*N /µ*N Vj ) (γ*j-1 - µ*j-1V1)/βj-1

j ) 2, 3, ..., N (18)

3. Practical Applications Three study cases were taken from Seader and Henley13 to show the reliability of a robust bubble-point scheme coupled with a BP method, such as the one proposed by Wang and Henke:2 (1) a typical distillation problem, where it is suggested that it be solved by the original Wang-Henke method2 because of the narrow boiling properties of the mixture involved; (2) an absorption problem, where Seader and Henley13 suggest to solve it by the Burningham-Otto sum rates method,9 because a BP method will fail because of the fact that the calculation of the stage temperature by a normal bubble-point temperature iteration scheme is sensitive to liquid-phase composition and the stage energy balance, eq 5, and is much more sensitive to stage temperatures than to interstage flow rates; and (3) a reboiled-absorption problem, where Seader and Henley13 suggested that it be solved using the SC method proposed by Naphtali and Sandholm,3 because this type of separator is like an absorber or stripper in one section and a fractionator in another section. The characteristics of these study cases are presented in Table 1. In this study, the Soave modification11 of the RedlichKwong equation of state is used to generate separation factors and mixture vapor and liquid enthalpies. The

Ind. Eng. Chem. Res., Vol. 42, No. 1, 2003 179 Table 1. Characteristics of the Study Cases

top pressure [kPa] reflux drum pressure [kPa] reboiler pressure [kPa] stage pressure drop [kPa] stages condenser type vapor top product [kmol/h] liquid top product [kmol/h] bottom product [kmol/h] liquid side streams [kmol/h] stage number vapor side streams [kmol/h] stage number heat extraction [kJ/h] stage number reflux ratio feeds methane (C1) [kmol/h] ethane (C2) [kmol/h] propane (C3) [kmol/h] n-butane (C4) [kmol/h] n-pentane (C5) [kmol/h] n-hexane (C6) [kmol/h] aborbent oil [kmol/h] temperature [K] pressure [kPa] stage number

example 1

example 2

example 3

1654.23 1640.05 1675.50 1.38 16 partial 6.818 2.272

2757.18

2757.18

0.0 6

0 13

350.000 1.364 3 16.818 13 211011.2 3 7.5 2

1 1.136 6.364 8.636 2.273 0.227

0.227 2.727 8.182 13.636 2.045

349.78 2067.94 6

383.11 1895.63 9

Soave equation may be written as

f(Z) ) Z3 - Z2 + (A - B - B2)Z - AB ) 0 (19) where Z is the compressibility factor and A and B are given by the following mixture rules:

A)

1

p

9(21/3 - 1) T2

B)

[ ] NC

xi ∑ i)1

TcixRi

(20)

ci

21/3 - 1 p NC Tci xi 3 Ti)1 pci



2 72.727 168.181 109.091 11.364 2.273

1 0.000 0.000 0.000 0.068 1.073

2 72.727 168.181 109.091 11.364 2.273

74.623 305.33 2757.18 1

0.000 313.67 2757.18 6

226.132 305.33 2757.18 1

0.000 313.67 2757.18 7

The pure-component physical properties required in the calculation were taken from work by Reid et al.18 To initialize the separation factors, the Wilson19 correlation was used, that is

()

ln Ki ) ln

2

xp

1 0.000 0.000 0.000 0.023 0.354

(21)

Here xi is the mole fraction in either the vapor or liquid phase. According to Soave11

xRi ) 1 + (0.480 + 1.547ωi - 0.176ωi )[1 2

(T/Tci)1/2] (22)

where Tci, pci, and ωi are the component i critical temperature, pressure, and acentric factor, respectively. To obtain values of thermodynamic properties, eq 19 is first solved for Z. When eq 19 has three real roots at fixed composition, temperature, and pressure, the largest represents the vapor phase and the smallest represents the liquid phase, but when eq 19 has only one real root, we need to use a pseudoroot for the unfeasible phase to evaluate the thermodynamical properties. In this work we use the heuristic and iterative procedure developed by Mathias et al.17 to calculate the density roots of the Soave equation of state. The corresponding values of Z are then used to calculate a fugacity coefficient and departure function for each phase by the appropriate equations listed, for example, in work by Reid et al.18 Separation factors, Ki’s, are then given by Ki ) φLi /φVi , where φi is the fugacity coefficient of component i in either the liquid (L) or vapor (V) phase.

p ci p

(

+ 5.373(1 + ωi) 1 -

)

Tci T

(23)

The initial set of vapor rates was established based on the assumption of constant molal interstage flows. Temperatures were initialized by assuming a linear variation between a top temperature estimate and a bottom temperature estimate. The following convergence criterion based on successive sets of the stage temperatures was considered: N

τ)

(Tkj - Tk-1 )2 e 0.01N ∑ j j)1

(24)

With the above considerations, the following numerical simulation was carried out. 1. Distillation. The simulation of the column with initial temperature values of 300 K for stage 1 and 450 K for stage 16 took seven iterations. The complete set of results is presented in Figure 3. Figure 4 shows the temperature and interstage vapor profiles at the initialization step, after the first iteration, and in the solution. Note that in this case after the first iteration the values of the profiles are so close to the ones of the solution that only seven iterations are required to solve the problem. 2. Absorption. Initial assumptions for the top-stage and bottom-stage temperatures were 305.33 and 313.67 K, respectively (corresponding to the liquid feed stream and vapor feed stream, respectively). For the simulation of this case, the absorbent oil was taken as n-nonane. This column took 11 iterations to reach the solution. The complete set of results is presented in Figure 5. According to Figure 6, it can be seen that after the first iteration the temperature and vapor profiles are far

180

Ind. Eng. Chem. Res., Vol. 42, No. 1, 2003

Figure 5. Specifications and results for the absorption case.

Figure 3. Specifications and results for the distillation case.

Figure 6. Temperature and interstage vapor profiles for the absorption case.

Figure 4. Temperature and interstage vapor profiles for the distillation case.

away from the solution, but this situation is corrected after the second iteration where the temperature and vapor profiles tend toward the solution. Seader and Henley13 present the solution of this column based on the Burningham-Otto9 solution scheme, assuming that the thermodynamic properties are independent of phase composition, using temperature polynomials. Seader and Henley13 reported seven iterations to reach the

solution accepting a convergence criterion of τ ) 0.0217; this criterion was reached with the proposed scheme at iteration 9. 3. Reboiled Absorption. Initial assumptions for the top-stage and bottom-stage temperatures were 338.67 and 449.78 K, respectively. This case took 15 iterations to reach the solution. The complete set of results is presented in Figure 7. According to Figure 8, after the first iteration, the temperatures of stages 1-7 are close to the ones of the solution; however, temperatures in stages 8-13 are far from the ones of the solution. This type of behavior evidences that the top section of this column is like a distillation case and the bottom section behaves as an absorption case and it will be very

Ind. Eng. Chem. Res., Vol. 42, No. 1, 2003 181

Figure 7. Specifications and results for the reboiled-absorption case.

Figure 8. Temperature and interstage vapor profiles for the reboiled-absorption case.

difficult to reach the solution. Seader and Henley13 present the solution of this column based on the Naphtali-Sandholm3 solution scheme, assuming that the thermodynamic properties are independent of the phase composition, using temperature polynomials. Seader and Henley13 reported seven iterations to reach the solution and also present the product compositions using the Soave modification11 of the Redlich-Kwong equation of state, but they did not comment if the convergence

Figure 9. Convergence patterns for the distillation, absorption, and reboiled-absorption cases.

properties are affected by introducing this type of calculation. Figure 9 shows the convergence patterns for the three cases studied, showing that a robust bubble-point temperature iteration scheme, as proposed in this work, integrated in a BP method, such as that proposed by Wang and Henke,2 is a stable enough method to solve all kinds of vapor-liquid multistage separations, because its ability of convergence by approaching to the ultimate solution is monotonic and without either oscillations or divergences. On the other hand, the ED methods, like the BP method, can be seen as two nested iteration loops. One disadvantage of these schemes is that an increase of the number of iterations in the inner loop may exist by means of a decrease of the number of iterations in the outer loop. Figure 10 shows the number of calls to the physical property package required by the bubble-point temperature iteration scheme in each iteration of the BP method. In general, as the solution is reached, the number of iterations in the bubble-point temperature calculation is decreased. The reported number of calls to the physical property package include the update of the separation factors. This feature results in a method in which in each iteration the robust scheme uses as the initial point the previous values of the BP method and in which the iterative step could move toward the solution. 4. Conclusions It has been shown that a judicious application of modern constructive phase equilibria techniques constitutes a useful tool to understand, systematize, and improve the existing ad hoc multicomponent multistage vapor-liquid separation simulation methods. By doing this, it has been found that a robust bubble-point

182

Ind. Eng. Chem. Res., Vol. 42, No. 1, 2003

Literature Cited

Figure 10. Relative calls to the physical property package based on the first iteration, for the distillation, absorption, and reboiledabsorption cases.

temperature iteration scheme coupled with a BP method is stable for wide boiling point mixtures, well beyond the point where different zones are presented in a column. Furthermore, the incorporation of a robust bubble-point temperature iteration scheme to a BP method has the advantage of yielding solutions to all different types of tower processes, including both absorption and distillation operations. Finally, a robust bubble-point temperature iteration scheme coupled with a BP method gives numerically stable solutions with a decrease in the number of iterations in the bubble-point temperature calculation as the solution is reached. Furthermore, the improved BP method is stable in the case where the thermodynamic properties are a strong function of the composition of the mixture, as is the case when an equation of state is used.

(1) Amundson, N. R.; Pontinen, A. J. Multicomponent Distillation Calculations on a Large Scale Computer. Ind. Eng. Chem. 1958, 50, 730. (2) Wang, J. C.; Henke, G. E. Tridiagonal Matrix for Distillation. Hydrocarbon Process. 1966, 45 (8), 155. (3) Naphthali, L. M.; Sandholm, D. P. Multicomponent Separation Calculations by Linearization. AIChE J. 1971, 17, 148. (4) Goldstein, R. P.; Stanfield, R. B. Flexible Method for the Solution of Distillation Design Problems Using the NewtonRaphson Technique. Ind. Eng. Chem. Process Des. Dev. 1970, 9, 78. (5) Friday, J. R.; Smith, B. D. An Analysis of the Equilibrium Stage Separation ProblemsFormulation and Convergence. AIChE J. 1964, 10, 698. (6) Tomich, J. F. A New Simulation Method for Equilibrium Stage Processes. AIChE J. 1970, 16, 229. (7) Orbach, O.; Crowe, C. M.; Johnson, A. I. Multicomponent Separation Calculations by the Modified Form of Newton’s Method. Chem. Eng. J. 1972, 3, 176. (8) Boston, J. F.; Sullivan, S. L. A New Class of Solution Methods for Multicomponent, Multistage Separation Processes. Can. J. Chem. Eng. 1974, 52, 52. (9) Burningham, D. W.; Otto, F. D. Which Computer Design for Absorbers? Hydrocarbon Process. 1967, 46 (10), 163. (10) King, J. C. Separation Processes, 2nd ed.; McGraw-Hill Book Co.: New York, 1980. (11) Soave, G. Equilibrium Constants from a Modified RedlichKwong Equation of State. Chem. Eng. Sci. 1972, 27, 1197. (12) Chapra, S. C.; Canale, R. P. Numerical Methods for Engineers: with Programming and Software Applications; McGrawHill: New York, 1998. (13) Seader, J. D.; Henley, E. J. Separation Process Principles; Wiley: New York, 1998. (14) Ortega, J. M.; Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables; Academic Press: New York, 1970. (15) Michelsen, M. L. The Isothermal Flash Problem. Part II: Phase-split Calculation. Fluid Phase Equilib. 1982, 9, 21. (16) Heidemann, R. A.; Michelsen, M. L. Instability of Successive Substitution. Ind. Eng. Chem. Res., 1995, 34, 958. (17) Mathias, P. M.; Boston, J. F.; Watanasiri, S. Effective Utilization of Equations of State for Thermodynamic Properties in Process Simulation. AIChE J. 1984, 30, 182. (18) (18) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases & Liquids, 4th ed.; McGraw-Hill Book Co.: New York, 1986. (19) Wilson, G. A Modified Redlich-Kwong Equation of State. Application to General Physical Data Calculations. Presented at the American Institute of Chemical Engineers National Meeting, Cleveland, OH, 1968; Paper 15c.

Received for review November 1, 2001 Revised manuscript received August 23, 2002 Accepted October 8, 2002 IE0108898