3792
Ind. Eng. Chem. Res. 2001, 40, 3792-3800
Efficient, Robust, and Reliable Single-Stage Vapor-Liquid Flash Calculations Rosendo Monroy-Loperena† Simulacio´ n Molecular, Instituto Mexicano del Petro´ leo, Eje Central La´ zaro Ca´ rdenas 152, 07730 Distrito Federal, Me´ xico
In this paper the problem of single-stage vapor-liquid flash calculations using equations of state is addressed. A solution scheme in the framework of complex domain is used in conjunction with a Newton-based formulation in block algebra, providing an efficient, robust, and reliable computation scheme. The new scheme minimizes the number of iterations, because it always ensures a direct-forward step toward the solution, making fewer calls to thermodynamic property evaluations during iterations and avoiding the use of pseudophases along the calculation path. Convergence to unambiguous results is assured for all cases of practical interest, because derivative discontinuities inherent to real-domain models, including the need for bounding mole fractions to lie between 0 and 1, are avoided. The efficiency and robustness of the proposed scheme are shown by calculations in the retrograde region and in the construction of phase envelopes for model mixtures. 1. Introduction Over the last 5 decades, a number of methods for solving the single-stage vapor-liquid flash (SSVLF) separation process have been developed (cf. work by King,1 Henley and Seader,2 and Heidemann3). All of the above methods exhibit definite strengths and weaknesses, depending upon the particular problem being treated. This paper describes the application of a Newton-type algorithm carried out in the complex domain. The main feature of this approach is the dependable reaching to the solution, regardless of the conditions where it is applied. Equations of state (EOSs), verbi gratia the volumetric relations between pressure, molar volume, and absolute temperature, have played a central role in thermodynamic modeling of the vapor-liquid equilibrium (VLE) of hydrocarbon fluids, especially at moderate and high pressures. With recent developments in EOS modeling, these models are becoming useful tools for the correlation and prediction of VLE of highly nonideal mixtures over broad ranges of pressure and temperature. The transformation of phase equilibrium modeling from activity coefficient models to EOSs is largely the result of the recently developed class of mixing rules that allows the use of liquid activity coefficient models in the EOS formalism. The implications of this transformation are far-reaching: for an EOS it offers a unified approach in thermodynamic property modeling. In contrast to the use of liquid excess free-energy or activity coefficient models, including the critical region, there is no need to choose an arbitrary reference state (which is especially troublesome when hypothetical states are used), and other thermodynamic properties, such as volumetric and calorimetric properties of mixtures, can also be obtained from the same EOS model (see the interesting work by Orbey and Sandler4 or the review by Ghosh5 on this topic). One of the main drawbacks of current SSVLF calculational procedures is the handling of phase specifica† E-mail:
[email protected]. Tel: +52-5333-8105. Fax: +52-5333-6239.
tions that do not actually exist based on the EOS calculations. Because the SSVLF is solved by an iteration method, it is frequently encountered that, during the course of iterations, the thermodynamic routines using an EOS are required to calculate the properties of a phase, which may not really exist, and an “appropriate” density root cannot be found according to its specifications (composition of the mixture, temperature, and pressure). The problem of converging to physically nonexistent conditions is generally known as the trivial root problem. In these situations, the problem can be addressed either in the high-level process calculations6-9 by changing the values of the iteration variables, i.e., phase composition, or in the thermodynamic routines10-16 by using pseudoroots for an unfeasible phase to evaluate the thermodynamical properties. An alternative approach to handling the single phase is to perform a stability analysis by means of the tangent plane distance (TPD).17 The TPD is simply the difference between the Gibbs free energies of a system in a single phase and in two phases, in which the amount of the second phase (that is, the trial phase) is small. Because at constant temperature and pressure the more stable of the tested phases has a lower Gibbs free energy, the single-phase state is stable if the TPD is positive. The procedure has the advantage of giving an initial estimate of the composition of the second phase if the single phase is unstable. However, its implementation is not a trivial task, especially in the retrograde region, because the TPD needs to be formulated as a minimization problem where the procedure may be sensitive to the initial guesses. A novel approach for solving phase equilibrium problems, which uses numerical methods that operate in the complex domain, was recently suggested by Lucia and co-workers.18-26 Their approach shows the following advantages: (A) the ability to find real-valued solutions from starting points in the complex domain that cannot be calculated using real initial values; (B) smoothness of the conservation equations and physical property models across phase boundaries that removes the need for ad hoc model manipulation in systems in which the
10.1021/ie000720e CCC: $20.00 © 2001 American Chemical Society Published on Web 07/27/2001
Ind. Eng. Chem. Res., Vol. 40, No. 17, 2001 3793
number of phases may change as a process parameter is varied or during iteration; (C) termination of simulation in a finite number of iterations to a real- or complexvalued steady-state solution that provides some physical insight about the solution. The results obtained by them involving phase equilibria, in order, are in regards to the following: (1) to a dynamical fixed vaporization isothermal flash process (Lucia et al.21), where they showed that both dynamical stable real-valued and dynamical stable complex-valued steady states can exist and the existence of stable complex-valued steady states is a simple consequence of the bifurcation of a pair of real-valued solutions into the complex domain; (2) to distillation calculations (Taylor et al.24), by carrying out Newton-based multistage distillation calculations entirely in the complex domain (they showed that this removes discontinuities inherent in real domain modes including conditions for bounding mole fractions to lie between 0 and 1); (3) to a fixed vaporization isothermal flash process (Gow et al.25), where VLE for refrigerant mixtures modeled by Soave’s27 modification of the Redlich-Kwong (SRK) EOS are studied, giving particular attention to root finding and root assignment of the cubic EOS to the vapor and liquid phases in the case where all roots to the EOS are complex-valued, yielding correct results, even in retrograde regions. Single-stage flash calculations may result in intrincate phase behavior, which includes multiple liquid and solid phases. Chemical reaction equilibrium may also occur. In this paper we limit the scope to nonreactive VLE because this covers the vast majority of applications in chemical processes. The main objective of the present work is to perform the SSVLF calculations in the retrograde region by carrying out a Newton-based formulation almost entirely in the complex domain, to avoid well-known problems found when the SSVLF calculations are carried out in the real domain (i.e., trivial solutions). For this purpose a previously proposed formulation of Monroy-Loperena and Duran28 is extended here to avoid a rank-deficient Jacobian, that is, a condition to apply a Newton-type iteration procedure in the complex domain to prevent periodic and chaotic behaviors.25 The approach of Monroy-Loperena and Duran28 in the real domain has shown remarkable flexibility, speed, and robustness, even for simulating multicomponent two-phase equilibria in supercritical asymmetric solutions that contain very large and very small components, such as polymers and supercritical solvents or monomers.29 The paper is organized in the following way. A mathematical model for the one-stage vapor-liquid flash is present in section 2. In section 3, the equation solution to solve simultaneously all of the material balance, equilibrium, constitutive and energy balance equations using a Newton iterative procedure in the framework of block algebra. Practical applications for retrograde condensation diagnosis, phase envelope construction, and VLE are presented in section 4. Comments about the SSVLF in the complex domain are given in section 5. Finally, the work is closed with some concluding remarks. 2. Equilibrium Single-Stage Vapor-Liquid Model In reference to Figure 1, the equations that describe a single-stage flash, for nonreactive VLE, are as follows:
Figure 1. Schematic representation of a SSVLF separation.
Component mass balance fi1 ) (1 - Φ)xi + Φyi - zi ) 0
for i ) 1, ..., N
(1)
Phase equilibrium fi2 ) Kixi - yi ) 0
for i ) 1, ..., N
(2)
Constitutive N
f3)
∑ i)1
N
xi -
yi ) 0 ∑ i)1
(3)
Energy balance f4 ) (1 - Φ)ΠL + ΦΠV - ΠF ) 0
(4)
where Φ is the ratio of the total mole flow rate of the vapor to that of the feed stream, usually named as the vapor fraction; N is the number of species in the mixture; xi, yi, and zi are the mole fractions of component i (i ) 1, ..., N), respectively, in the liquid phase, vapor phase, and feed stream; Ki is the separation factor of component i which is defined as yi/xi; and Πσ can be enthalpy, entropy, or internal energy of the stream σ where eq 4 is applied. In principle, SSVLF calculations may be specified by any set of independent variables that are consistent with the phase rule. Chemical engineering practice has evolved a few specific formulations that are most useful. In all of these formulations, the feed compositions of each component in a multicomponent mixture are specified. Following this, the phase rule stipulates that the specification of two additional independent variables completely and uniquely determines the state of the system. Six pairs of specified independent variables define the SSVLF types that are most useful in chemical engineering practice: type I, ΠF, P; type II, ΠF, T; type III, ΠF,Φ; type IV, Φ, T; type V, Φ, P; type VI, T, P. T and P respectively are the temperature and pressure in the flash drum. As stated, Π could be enthalpy, entropy, or internal energy. For instance, when the feed stream is reduced adiabatically across a valve, the problem is a type I flash with Π as enthalpy. If Π is taken as entropy, type I flash is a compressor or expander calculation. Taking Π as internal energy is a natural outcome of applying the first law of thermodynamics that often occurs in dynamic simulations. In general, Ki’s and Π’s are obtained from a thermodynamic model and xi (i ) 1, ..., N) and yi (i ) 1, ..., N), and depending on the type of problem, two variables of the set {Φ, ΠF, T, P} are treated as unknowns. Because there are N material balance equations (1), N equilib-
3794
Ind. Eng. Chem. Res., Vol. 40, No. 17, 2001
ΠF, T, and P (in that specific order); and the last row represents the corresponding derivatives of eq 4 for the energy balance with respect to xi and yi (i ) 1, ..., N), Φ, ΠF, T, and P (in that specific order). To solve the system in eq 7, the Jacobian matrix will be constructed by taking only two of the four last rows of the system, depending on the type of SSVLF problem. From the incidence matrix of Figure 2, the following equivalent block representation of the Jacobian can be defined:
J)
Figure 2. Incidence matrix of the generated Jacobian.
rium equations (2), one constitutive equation (3), and one energy balance, we end up with a system of 2N + 2 nonlinear algebraic equations.
B(v f b) ) 0
vk + b u(v bk) b v k+1 ) b
J(v bk)u b(v bk) ) -B(v f bk)
Bf )
(9)
[ ]
(10)
-g b1 -g b2
u, that is, where b u1 contains the first N rows of vector b u2 contains the changes of variables xi (i ) 1, ..., N), and b the remaining N + 2 rows of vector b u, that is, the change b1 of variables yi (i ) 1, ..., N), R, and β. Similarly, -g contains the first N rows of vector Bf and -g b2 contains the remaining N + 2 rows of vector B. f Therefore, we rewrite eq 7 as
[ ][ ] [ ] A B b u1 b g1 ) C D b u2 b g2
(11)
As a result of block lower-upper (LU) factorization of the Jacobian J in eq 8 and the required block linear algebra in eq 11, one obtains an equivalent decomposed linear system of equations:
(7)
vk where B(v f bk) is the vector of functions (5) evaluated at b and J(v bk) is the Jacobian matrix of partial derivatives of Bf with respect to b v, evaluated at b vk. One specific iterative approach for solving the system of linear equations (7) was proposed by Monroy-Loperena and Duran.28 They realized that the Jacobian matrix in eq 7 has easily identifiable matrix blocks that lead themselves to block linear algebra treatment and hence to a problem decomposition. Figure 2 shows the incidence matrix (i.e., the structure) of the partial derivatives of Bwith f respect to b v, where the first N rows represent the derivatives of equations (1) with respect to xi and yi (i ) 1, ..., N), Φ, ΠF, T, and P (in that specific order); the next N rows represent the derivatives of equations (2) for phase equilibrium with respect to xi and yi (i ) 1, ..., N), Φ, ΠF, T, and P (in that specific order); the next row represents the derivatives of eq 3 of constitutive with respect to xi and yi (i ) 1, ..., N), Φ,
[]
and
(6)
where k is the iteration index, b vk is the current solution estimate, b u(v bk) is the predicted vector changes as a vk+1 is the vector of new solution function of b vk, and b estimates. The solution step b u(v bk) at iteration k is obtained by solving the following system of linear equations:
(8)
b u1 b u2
b u)
(5)
where Bf is the vector of functions (i.e., eqs 1-4) and b v is the vector of 2N + 2 unknowns, b v ) {(xi: i ) 1, ..., N), (yi: i ) 1, ..., N), R, β}T, where R and β represent two variables from the set {Φ, ΠF, T, P} according to the type of SSVLF problem to be solved. The system of equations (5) can be solved simultaneously using a Newton iterative equation:
A B C D
where A, B, C, and D are the matrix blocks that have the dimensions of N × N, N × (N + 2), (N + 2) × N, and (N + 2) × (N + 2), respectively. A special feature of Jacobian (8) is that A is a diagonal matrix and B is an almost diagonal matrix. For the efficient solution of eq 7, exploiting the special structure identified in (8), let us define the following partition for the b u and Bf vectors:
3. Equation Solution The model defined above leads to a set of 2N + 2 algebraic equations, which can be written as
[ ]
b2 ) b g 2 - CA-1b g1 (D - CA-1B)u
(12)
b1 - Bu b2) b u1 ) A-1(g
(13)
and
In this decomposition, the variables b u2 are obtained 1 independently of b u by solving the reduced linear subsystem of N + 2 equations (12). Once b u2 has been 1 determined, b u is obtained explicitly from eq 13. The actual matrix elements of these systems of equations are shown in the appendix. Regarding the merits of the equation solution approach presented above, the following comments must be done: (i) Obtaining A-1 is simple because A is a diagonal matrix. Note that following the precedent choice of the variable vector
b v ) {(xi: i ) 1, ..., N), (yi: i ) 1, ..., N), R, β}T
Ind. Eng. Chem. Res., Vol. 40, No. 17, 2001 3795
yields A ) (1 - Φ)I, where I denotes the identity matrix, and A-1 ) (1 - Φ)-1I; thus, for the case where Φ f 1, then A goes to singular, and in the worst case, when Φ ) 1 (dew point), A is singular. This problem can been avoided by choosing the variable vector as
b v ) {(yi: i ) 1, ..., N), (xi: i ) 1, ..., N), R, β}T then the A matrix will be A ) ΦI, eliminating this problem. As a general rule, if |Φ| > |1 - Φ|, choose
b v ) {(yi: i ) 1, ..., N), (xi: i ) 1, ..., N), R, β}T Otherwise, choose
b v ) {(xi: i ) 1, ..., N), (yi: i ) 1, ..., N), R, β}T This will also avoid ill-conditioning of systems (12) and (13). (ii) Solution of eqs 12 and 13 requires less computing and is numerically more stable than that of the original system (7) because of the reduced dimensionality. 4. Practical Applications 4.1. Thermodynamics. In this study, the Soave modification27 of the Redlich-Kwong EOS is used to generate equilibrium ratios, Ki’s, mixture Π’s, and appropriate derivative properties. The Soave equation may be written as
f(Z) ) Z3 - Z2 + (A - B - B2)Z - AB ) 0 (14) where Z is the compressibility factor and A and B are given by the following mixture rules:
A)
P
1
[ ] N
9(21/3 - 1) T2 B)
2
1/3
xi ∑ i)1
-1P
3
N
Tc,ixRi
xPc,i
Tc,i
∑xi Ti)1 P
2
forward way to calculate thermodynamic properties of mixtures, especially the fugacity coefficients and their derivatives. Equilibrium ratios, 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. To avoid unnecessary operations, the equilibrium ratios are better expressed as Ki ) exp(ln φL - ln φV).33 4.2. Solving for the Compressibility Factor. One of the most important aspects of solving the SSVLF problem modeled by a single EOS is the assignment of a qualitatively correct compressibility root to each of the equilibrium phases. An incorrect root assignment can lead to erroneous thermodynamic properties and is often computationally catastrophic, causing convergence to a trivial solution or even divergence. In either case, this is viewed as algorithmic failure. The worst case is when the EOS only has one real-valued root instead of three; the two missing real-valued roots or solutions have to be seen as they bifurcate into the complex domain.18 Regarding this, we actually initialize the numerical liquid compressibility calculation by setting Z ) B + 0.1i and initializing the vapor compressibility by setting Z ) 1 + 0.1i, by using a Newton-Raphson numerical approach. This allows the iterations to enter the complex domain and, therefore, to find a complex root if one exists. If no complex root exists, Newton-Raphson’s method still finds a real-valued root. Thus, at least two distinct compressibilities will almost always be found (even if only one of them is real). 4.3. Initialization. When no initial estimates are available, the Wilson34 K separation factor expressions can be used to generate approximate values of phase properties
ln KWilson ) ln(Pc,i/P) + 5.373(1 + ωi)(1 - Tc,i/T) (18) i coupled with the Rachford and Rice equation:35
(15) N
∑ i)1
zi(KWilson - 1) i
(KWilson - 1)Φ + 1 i
(16)
c,i
Here xi is the mole fraction in either the vapor or liquid phase. According to Soave26
yi )
xRi ) 1 + (0.480 + 1.547ωi - 0.176ωi )(1 - xT/Tc,i) where Tc,i, Pc,i, and ωi are the component i critical temperature, critical pressure and the acentric factor, respectively. To obtain values of thermodynamic properties, eq 14 is first solved for Z. In the real domain when eq 14 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 14 has only one real root, we need to use pseudoroots10-16 for the unfeasible phase to evaluate the thermodynamical properties. In the complex domain, a different approach is taken and will be fully described herein. The corresponding values of Z are then used to calculate a fugacity coefficient and a departure function for each phase by the appropriate equations listed, for example, in Reid et al.,30 or even more, by use of the Helmholtz energy surface,31,32 which is an efficient and straight-
(19)
The initial guesses for the vapor and liquid compositions are obtained by
2
(17)
)0
xi )
ziKWilson i (KWilson - 1)Φ + 1 i zi (KWilson i
- 1)Φ + 1
(20)
(21)
Combined with an expression for the ideal gas heat capacity as a function of temperature, this enables us to calculate all relevant thermodynamic properties. As stated by Michelsen,36 the approximation obtained by means of the Wilson approximation is normally adequate when the vapor fraction at the solution is not very small; that is because we are assuming that
ln Ki = ln φLi = ln KWilson i The most obvious defect is the severe misrepresentation of liquid volumes, which are evaluated to be identically zero. Its main advantage is that it enables a simple and safe determination of an initial approximation to the solution.
3796
Ind. Eng. Chem. Res., Vol. 40, No. 17, 2001
continuously changing states in isenthalpic or isentropic compression and expansion. To describe the procedure to calculate the partial derivatives (∂v b/∂δ)γ, let us rewrite eq 5 by adding two new functions:
f 5 ) f2N+3 ) γ - γ*
(24)
f 6 ) f2N+4 ) δ - δ*
(25)
where γ* and δ* are the two specifications that define the type of SSVLF calculation, leading to a set of 2N + 4 algebraic equations. Now Bf is the vector of functions (eqs 1-4, 24, and 25), and b v is the vector of 2N + 4 unknowns, either
b v ) {(xi: i ) 1, ..., N), (yi: i ) 1, ..., N), R, β, γ, δ}T or
b v ) {(yi: i ) 1, ..., N), (xi: i ) 1, ..., N), R, β, γ, δ}T
Figure 3. Retrograde condensation. Table 1. Diagnosis of Retrograde Condensation (See Figure 3) point
P
T
Φ
(∂P/∂T)Φ
(∂Φ/∂T)P
A B C
PA PB ) PA PC
TA TB TC ) T B
1.00 1.00 1.00
>0 0
0 >0
Our recommended approach, to get initial estimates, is to converge the Rachford-Rice equation (19) using the Wilson K separation factors (18). The resulting estimate of temperature, pressure, or vaporization and phase composition cannot be expected to be very accurate but need to be reasonable. In the case that the results are not reasonable, we suggest that the calculation be repeated on some reasonable condition (i.e., in a lower pressure) and then these estimates be used to initialize the proposed method. On the other hand, if we need a temperature, the guess is assumed to be 0.7 times the average critical temperature:
T ) 0.7
∑ziTc,i
(22)
Equation 22 is based on the fact that, for almost pure fluids, the normal boiling temperature is approximately 0.7 times the critical temperature. If we need a pressure, the guess is assumed to be 0.3 times the average critical pressure:
P ) 0.3
∑ziPc,i
(23)
It is important to note that the initialization process is carried out in the real domain. 4.4. Retrograde Condensation Diagnosis. As stated by Asselineau et al.37 in a mixture such as that represented in Figure 3, it can be seen that both cricondetherm and cricondebar points lie on the dew curve. This means that the same pressure corresponds to two dew temperatures (isobar AB in Figure 3) and the same temperature to two dew pressures (isotherm BC in Figure 3). A distinction between the dew states, points A-C, can be made from the sign of the derivatives, as is shown in Table 1. In the same way the cricondetherm and cricondebar can be located from the value of (∂P/∂T)Φ. Knowing the derivatives is also useful for computing the evolution of a mixture through
(see section 3). R, β, γ, and δ are the elements of the set {Φ, ΠF, T, P}. For example, in a type IV problem R ) P, β ) ΠF, γ ) Φ, and δ ) T. To obtain the partial derivatives (∂v b/∂δ)γ, it is known that for any infinitesimal change of the equilibrium state, b ve, the component mass balance, phase equilibrium, eliminate constitutive and energy balance relations are preserved. We can thus write
[ ]
∂fk T [dv be] ) dfk ) 0 ∂v b
for k ) 1, ..., 2N + 2
(26)
b] is the vector of the partial derivatives of where [∂fk/∂v v. fk with respect to the vector of unknowns b Because we are calculating partial derivatives, taking γ as constant, we have
[ ]
∂f2N+3 T [dv be] ) df2N+3 ) dγ ) 0 ∂v b
(27)
and the expression of δ differential
[ ]
∂f2N+4 T [dv be] ) df2N+4 ) dδ ∂v b
(28)
b] and [∂f2N+4/∂v b] are the vectors of partial where [∂f2N+3/∂v derivatives of γ and δ with respect to the vector of unknowns b v. From eqs 26-28, after division by dδ, a linear system is obtained, where the unknowns are the partial derivatives of the equilibrium conditions
[ ][( ) ] ∂Bf ∂v b
∂v be ∂δ
γ
)C Bδ
(29)
where Bf is the vector previously defined by eqs 1-4, 24, and 25 and C B δ is
C B δ ) {0, ..., 0, 1}T
(30)
where the first 2N + 3 elements are zero and the last one is unity. The resulting system (29) also has identifiable blocks as in eq 8, but A, B, C, and D now have the dimensions of N × N, N × (N + 4), (N + 4) × N, and (N + 4) × (N
Ind. Eng. Chem. Res., Vol. 40, No. 17, 2001 3797 Table 2. Compositions (Mole Fractions) of Systems Studied component N2 C1 C2 C3 n-C4 n-C5 n-C6 n-C7
I
II
0.9685
0.5000
0.0315
0.500
Table 3. Diagnosis of Retrograde Condensation for the Mixture of Figure 4
III
point
P
T
Φ
(∂P/∂T)Φ
(∂Φ/∂T)P
0.0140 0.9430 0.0270 0.0074 0.0049 0.0010 0.0027
A B C
6500.00 6500.00 2216.14
332.46 348.33 348.33
1.00 1.00 1.00
45.1595 -64.8523 120.3444
-24.7780 0.4153 0.1621
+ 4), respectively, so we can decompose the linear system of equations to use eqs 12 and 13. It must be pointed out that the construction of the system of linear equations (29) is straightforward, because the first 2N + 2 elements of the derivative b], k ) 1, ..., 2N + 2, have been calculated vectors [∂fk/∂v in the last step of the Newton iterative search and are the same whatever the data of the previously solved SSVLF problem may be and are the kind of partial derivatives we want. It is also obvious that the vectors of derivatives only have to be calculated with respect to the data used to define the problem and the last four column vectors have to be reordered to get the partial derivatives we want according to eqs 27-29. 4.5. Phase Envelope Construction. Handling derivatives of equilibrium conditions enables a straight sweeping process to be used for the construction of the phase envelope of a mixture. The method used to construct the phase envelope is similar to that proposed by Ziervogel and Poling,38 by stepping around the phase envelope in suitably small increments of T or P. Fast convergence is obtained if one does dew- or bubble-point pressure calculations (type IV problem) when the phase envelope is flat and dew- or bubble-point temperature calculations (type V problem) when the phase envelope is steep. These conditions have been stated as follows: (i) If |d ln P/d ln T| < 2, calculate dew- or bubblepoint pressures by the type IV problem. (ii) If |d ln P/d ln T| > 20, calculate dew- or bubblepoint pressures by the type V problem. When |d ln P/d ln T| is between 2 and 20, convergence is obtained regardless of the type of problem used (IV or V). 4.6. Numerical Examples. We have carried out several numerical calculations with three representative systems (two binary and one multicomponent) to illustrate the performance of the Newton iterative procedure in the framework of block algebra in the complex domain. Table 2 shows mixtures’ properties. Purecomponent properties were taken from Reid et al.30 In all of the cases, initial guesses were calculated in accordance with the procedure described in section 4.3, and the phase envelopes were calculated using the procedure described in section 4.5. 4.6.1. Retrograde Diagnosis. The binary system corresponds to an ethane/hepthane mixture. In this case we stated the following question: Does the mixture present retrograde condensation at 6500 kPa? The normal procedure to answer this question is to apply a type V SSVLF with Φ ) 1 and P ) 6500 kPa. After we solved the problem, we needed to know if we obtained point A or B, so we used the information of the partial derivatives described in Table 1. For this case, using the initialization procedure described in section 4.3, we found first point B; to find point A, we needed to
Figure 4. Phase envelope and retrograde condensation for mixture I. Table 4. Phase Calculations for the Mixture of Figure 5 point
temp [K]
pressure [kPa]
vaporization
iterations
1 2 3 4 5 6 7 8 9 10 11 12
470 450 470 490 430 450 470 490 410 450 490 450
6950 6500 6500 6500 6000 6000 6000 6000 5000 5000 5000 3000
0.0106 0.0832 0.2039 0.5341 0.0785 0.1919 0.3321 0.6400 0.1395 0.3513 0.7941 0.6059
3 3 4 7 3 3 4 5 3 4 5 4
initialize a problem type V at a pressure of 3000 kPa and a temperature of 328.33 K, that is, 20 K less than that obtained in point B. Point C is easily obtained by a type IV SSVLF, with initial estimates at a pressure of 1000 kPa and a temperature of 348.33 K. Results are presented in Table 3 and Figure 4. As a general rule, to obtain convergence to points such as A, we recommend using estimates generated at a lower pressure with the bubble temperature obtained at the specified pressure, minus some degrees; usually 20 K are good enough to ensure convergence. 4.6.2. Vapor-Liquid Calculations. The second problem involves mixture I with different feed compositions. Note that the cricondebar point now lies on the bubble-point curve. We carried out 12 VLE point calculations using the methodology outlined in the previous sections. Results are present in Table 4 and Figure 5. The third problem is a multicomponent mixture. No problems have arisen in the solution of the 12 equilib-
3798
Ind. Eng. Chem. Res., Vol. 40, No. 17, 2001
Figure 5. Vapor-liquid equilibria calculated for mixture II.
Figure 6. Vapor-liquid equilibria calculated for mixture III.
Table 5. Phase Calculations for the Mixture of Figure 6
ibility roots to be selected for the liquid and vapor phases. This, in turn, gives unequal fugacity coefficients, resulting in the annihilation of trivial solutions. (iv) The disadvantages of storage and numerical operations are compensated by the block algebra treatment, because the number of operations needed to solve a system of N linear equations is proportional to N3. (v) The use of the complex domain to solve SSVLF calculations is justified by the fact that a given problem has a solution when in the real domain does not.
point
temp [K]
pressure [kPa]
vaporization
iterations
1 2 3 4 5 6 7 8 9 10 11 12
230 230 210 230 250 210 230 250 210 230 250 230
8100 7100 6600 6600 6600 6100 6100 6100 5100 5100 5100 3000
0.9988 0.9900 0.9705 0.9881 0.9987 0.9452 0.9873 0.9980 0.9572 0.9874 0.9972 0.9909
4 4 7 5 2 7 5 2 6 6 3 7
rium points; results for this mixture are present in Table 5 and Figure 6. As a general observation, all points obtained have real-valued solutions, in spite of the fact that the Newton iterative procedure enters in the complex domain, and the number of iterations in all of the cases are relatively small. 5. Comments on the SSVLF Calculations in the Complex Domain The real advantage of the proposed approach in the solution of SSVLF can be summarized as follows: (i) The block linear algebra treatment in Newton’s method does not generate a Jacobian matrix of rank deficiency, because at any moment periodic and chaotic behaviors are found. This is important because it is not necessary to perform a special manipulation to force a normal reduction.25,39 (ii) The proposed method can be used with confidence, because numerical difficulties that would normally occur in the real domain, as derivative discontinuities, do not occur. (iii) In the real domain, many starting values converge to a trivial solution (i.e., equal phase composition for vapor and liquid). On the other hand, in the complex domain, allowing use of the complex solutions of the EOS always allows for distinctly different compress-
6. Conclusions A new alternative calculation procedure for SSVLF equilibrium problems in the framework of the complex domain has been presented. The proposed approach preserves the desired quadratic convergence rate and requires no special modification to preserve the iteration values with a meaningful context. As a result, the EOS manipulation requires no artificial solutions. Sample calculations show that convergence to unambiguous results is assured even in the retrograde region with simple initialization procedures. Appendix. Construction of the Matrices of the Proposed Method Let us define the following system of linear equations for the system (12):
Mu b2 ) b n
(A.1)
where the elements of the matrix M are as follows. Case A. The variable vector is chosen as
b v ) {(xi: i ) 1, ..., N), (yi: i ) 1, ..., N), R, β}T then
mi,j ) xi
( )(
)
∂Ki ∂Ki Φ - δij δijKi + xi ∂xj 1-Φ ∂xj
(A.2)
Ind. Eng. Chem. Res., Vol. 40, No. 17, 2001 3799
and the elements of vector b n are
where
δij ) mi,N+1 ) xi
mi,N+2 ) xi
∂Ki
-
( ∑(
}
∑ δikKi + xi ∂x
1 - Φk)1
-
∂β
N
1
δikKi + xi
1 - Φk)1
mN+1,j ) -
mN+1,N+1 ) -
mN+1,N+2 ) (1 - Φ)
∂ΠL
) )
k
mi,N+2 ) xi
∂Ki ∂f1k
∑
∂f1k
∂R
(
∑ k)1 ∂x
)
mN+2,N+1 ) -
mN+2,N+2 ) (1 - Φ)
∂ΠL ∂β
1
∂f1k
∑ 1 - Φk)1 ∂β ∂ΠV
+Φ
N
-
∂β
1
∑
1 - Φk)1
nN+1 )
2 gN+1
-
∂β
k
nN+2 )
)
∂Kj ∂xk
-
∑
k)1
g1k
(
1
∂R
+Φ
(A.8)
mN+2,N+2 ) (1 - Φ)
∂ΠL
nj ) g2j -
∂β 1
+Φ
N
(A.16)
∂R
(A.17)
∂yk ∂β (A.18)
1
∑
∂f1k
(A.19)
Φk)1 ∂R ∂ΠV
N
1
∂ΠV ∂fk
∑ ∂y
-
∂R
k)1
(
nN+2 )
1
N
∑
(A.22)
Φk)1 ∂β ∂ΠV
N
-
∂β
1
∂ΠV ∂fk
∑ ∂y
k)1
k
∂β
)
∂Kj
1
k
(A.23) (A.24)
N
g1k ∑ Φk)1
∑
(A.20) (A.21)
∂f1k
(
-
∂R
)
N
2 gN+2
k
∂ΠL ∂ΠV ∂xj ∂yj
2 + nN+1 ) gN+1
(A.11)
N
g1k -δjk + xj ∑ Φk)1 ∂y
(A.10)
g1k
k)1
∂ΠL
(A.13)
∂xk
)
1
(A.14)
Case B. The variable vector is chosen as
b v ) {(yi: i ) 1, ..., N), (xi: i ) 1, ..., N), R, β}T then
( )(
mN+1,N+2 ) (1 - Φ)
∂ΠL
mN+2,N+1 )
(A.12)
∂fj ∂fj 1 2 2 - uN+2 g1j - Φu2j - uN+1 1-Φ ∂R ∂β
mi,j ) δijKi + xi
-δik + xi
(A.7)
After the linear system (A.1) is solved for b u2, b u1 is obtained by
u1j )
) )
∂Ki ∂f1k
N
Φk)1
k
(A.25)
∂ΠV
(A.26)
∂yk
N
g1k ∑ 1 - Φk)1 N
2 gN+2
-
∂β
1
∑ k)1 ∂x
g1k δjkKj + xj 1
1
mN+1,N+1 )
∂ΠL ∂fk
(
N
∂Ki
mN+2,j ) (1 - Φ)
(A.9)
and the elements of vector b n are
nj ) g2j -
( ∑(
∑ -δik + xi ∂y
Φk)1
∂R
k
∂ΠV ∂ΠL mN+2,j ) Φ ∂yj ∂xj N
-
∂Ki ∂f1k
N
mN+1,j ) 1/Φ
1
∂ΠL ∂fk
N
-
1
∂R
(A.4)
(A.6)
1 - Φk)1 ∂R
+Φ
∂Ki
∂xk ∂β (A.5)
N
1
(A.3)
∂R
1 1-Φ
∂ΠV
∂R
mi,N+1 ) xi
∂Ki ∂f1k
N
1
∂R ∂Ki
{
if i ) j if i * j
1 0
)
∂Ki ∂Ki 1-Φ -δij + xi ∂xj 1 ∂yj
(A.15)
u1 is After the linear system (A.1) is solved for b u2, b obtained by
u1j )
(
1
)
1
∂fj ∂fj 1 1 2 2 gj - (1 - Φ)u2j - uN+1 - uN+2 Φ ∂R ∂β
(A.27)
Literature Cited (1) King, J. C. Separation Processes, 2nd ed.; McGraw-Hill Book Co.: New York, 1980. (2) Henley, E. J.; Seader, J. D. Equilibrium-Stage Separation Operations in Chemical Engineering; Wiley: New York, 1981. (3) Heidemann, R. A. Computation of High-Pressure Phase Equilibria. Fluid Phase Equilib. 1983, 14, 55. (4) Orbey, H.; Sandler, S. I. Modeling Vapor-Liquid Equilibria. Cubic Equations of State and Their Mixing Rules; Cambridge University Press: Cambridge, U.K., 1998. (5) Ghosh, P. Prediction of Vapor-Liquid Equilibria Using Peng-Robinson and Soave-Redlich-Kwong Equations of State. Chem. Eng. Technol. 1999, 22, 379. (6) Anderson, T. F. Molecular Thermodynamics of VaporLiquid and Liquid-Liquid Equilibria. Ph.D. Dissertation, University of California, Berkeley, 1978. (7) Poling, B. E.; Grens, E. A., II; Prausnitz, J. M. Thermodynamic Properties from a Cubic Equation of State: Avoid Trivial
3800
Ind. Eng. Chem. Res., Vol. 40, No. 17, 2001
Roots and Spurious Derivates. Ind. Eng. Chem. Process Des. Dev. 1981, 20, 127. (8) Veeranna, D.; Husain, A.; Subrahmanyam, S.; Sarkar, M. K. An Algorithm for Flash Calculations using an Equation of State. Comput. Chem. Eng. 1987, 11, 489. (9) Pasad, G. V.; Venkatarathnam, G. A Method for Avoiding Trivial Roots in Isothermal Flash Calculations Using Cubic Equations of State. Ind. Eng. Chem. Process Des. Dev. 1999, 38, 3530. (10) Mills, M. B.; Wills, M. J.; Bhirud, V. L. The Calculations of Density by the BWRS Equation of State in Process Simulation Context. AIChE J. 1980, 26, 902. (11) Gundersen, T. Numerical Aspects of the Implementation of Cubic Equations of State in Flash Calculations Routine. Comput. Chem. Eng. 1982, 6, 245. (12) Jovanovlc, S.; Paunovlc, R. Generating Appropriate Density Values from Cubic State Equation to Avoid False Unit K Values. Application to Distillation Problems. Ind. Eng. Chem. Process Des. Dev. 1984, 23, 801. (13) 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. (14) Gosset, R.; Heyen, G.; Kalitventzeff, B. An Efficient Algorithm to Solve Cubic Equations of State. Fluid Phase Equilib. 1986, 25, 51. (15) Stateva, R. P.; Tsvetkov, S. G.; Boyadjiev, C. B. An Approach Organizing the EOS Model in Process Simulators. AIChE J. 1990, 36, 132. (16) Zhao, E.; Saha, S. Applications of Complex Domain in Vapor-Liquid Equilibrium Calculations Using a Cubic Equation of State. Ind. Eng. Chem. Res. 1998, 37, 1625. (17) Michelsen, M. L. The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1. (18) Lucia, A.; Guo, X.; Richey, P. J.; Derebail, R. Simple Process Equations, Fixed-Point Methods, and Chaos. AIChE J. 1990, 36, 641. (19) Lucia, A.; Xu, J. Global Minima in Root Finding. In Recent Advances in Global Optimization; Floudas, C. A., Pardalos, P. M., Eds.; Princeton University Press: Princeton, NJ, 1992. (20) Lucia, A.; Taylor, R. Complex Iterative Solutions to Process Model Equations? Comput. Chem. Eng. 1992, 16S, 387. (21) Lucia, A.; Guo, X.; Wang, X. Process Simulation in the Complex Domain. AIChE J. 1993, 39, 461. (22) Sridhar, L. N.; Lucia, A. Process Analysis in the Complex Domain. AIChE J. 1995, 41, 585. (23) Lucia, A.; Wang, X. Complex Domain Process Dynamics. Ind. Eng. Chem. Res. 1995, 34, 202. (24) Taylor, R.; Achuthan, K.; Lucia, A. Complex Domain Distillation Calculations. Comput. Chem. Eng. 1996, 20, 93.
(25) Gow, A. S.; Guo, X.; Liu, D.; Lucia, A. Simulation of Refrigerant Phase Equilibria. Ind. Eng. Chem. Res. 1997, 36, 2841. (26) Lucia, A. Complex Domain Chemical Process Simulation in Theory and in Practice. Ind. Eng. Chem. Res. 2000, 39, 1713. (27) Soave, G. Equilibrium Constants from a Modified RedlichKwong Equation of State. Chem. Eng. Sci. 1972, 27, 1197. (28) Monroy-Loperena, R.; Duran, M. A. The Flash Problem: Solution by a Structured Approach. Avances en Ingenierı´a Quı´mica; Academia Mexicana de Investigacio´n y Docencia en Ingenierı´a Quı´mica: San Luis Potosı´, SLP, Me´xico, Sept 1987. (29) Cheng-keng, C.; Duran, M. A.; Radosz, M. Phase Equilibria in Polymer Solutions. Block-Algebra, Simultaneous Flash Algorithm Coupled with SAFT Equation of State, Applied to SingleStage Supercritical Antisolvent Fractionation of Polyethylene. Ind. Eng. Chem. Res. 1993, 32, 3123. (30) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases & Liquids, 4th ed.; McGraw-Hill Book Co.: New York, 1986. (31) Mollerup, J. M. Thermodynamic Properties from a Cubic Equation of State, Report SEP 8601, Department of Chemical Engineering, DTH: Lyngby, Denmark, 1986. (32) Mollerup, J. M.; Michelsen, M. L. Calculation of Thermodynamic Equilibrium Properties. Fluid Phase Equilib. 1992, 74, 1. (33) Shah, M. K.; Bishnoi, M. K. Multistage multicomponent calculations using thermodynamic properties evaluated by the SRK/PR equation of state. Can. J. Chem. Eng. 1978, 56, 478. (34) 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, 1968; Paper 15c. (35) Rachford, H. H.; Rice, J. D. Procedure for use of electronic digital computers in calculation flash vaporization hydrocarbon equilibrium. J. Pet. Technol. 1952, 4, 9. (36) Michelsen, L. M. State function based flash specifications. Fluid Phase Equilib. 1999, 158-160, 617. (37) Asselineau, L.; Bogdanic, G.; Vidal, J. A Versatile Algorithm for Calculating Vapor-Liquid Equilibria. Fluid Phase Equilib. 1979, 3, 273. (38) Ziervogel, R. G.; Poling, B. E. A simple method for constructing phase envelopes for multicomponent mixtures. Fluid Phase Equilib. 1983, 11, 127. (39) Lucia, A. Complex Domain Chemical Process Simulation in Theory and in Practice. Ind. Eng. Chem. Res. 2000, 39, 1713.
Received for review August 3, 2000 Accepted May 29, 2001 IE000720E
