6028
Ind. Eng. Chem. Res. 2001, 40, 6028-6038
GENERAL RESEARCH Modeling and Analysis of the Isothermal Flash Problem and Its Calculation with the Simulated Annealing Algorithm N. Henderson,* J. R. de Oliveira, Jr., H. P. Amaral Souto, and R. Pitanga Marques Instituto Polite´ cnico-UERJ, Caixa Postal 97282, Nova Friburgo, Rio de Janeiro, CEP 28601-970 Brazil
This paper describes a new formulation of the vapor-liquid equilibrium problem. The proposed model reveals a great mathematical similarity between the stability test problem and the isothermal flash problem, particularly in regions of interest in reservoir engineering near phase boundaries and in the vicinity of critical points. This similarity is mathematically studied and then explored to find a means of minimizing the Gibbs free energy in a more efficient way. The Peng-Robinson cubic equation of state was used to represent the phase behavior of multicomponent mixtures. The simulated annealing (SA) algorithm was used to carry out isothermal flash calculations and phase stability tests. The SA algorithm is essentially an iterative random search procedure whose moves are adapted by some suitable mechanism. It is robust and easy to implement. The algorithm was tested for several hydrocarbon mixtures. The few experimental data available in the literature were compared with theoretical predictions. The proposed methodology produced high-quality solutions. 1. Introduction Isothermal flash calculations are of particular importance in petroleum and chemical engineering. Miscible gas floods of oil reservoirs cannot be adequately modeled without a detailed flash phase behavior description. Furthermore, distillation, adsorption, and extraction are among the multistage processes where phase equilibrium calculations at specified values of P and T are required. For many years, the formulation of flash phase equilibrium was almost always restricted to the pioneering work of Rachford and Rice,1 who used the equilibrium ratios (K values) as the iteration variables. However, the slow convergence rate of the successive substitution method for nonideal mixtures prompted authors to devise algorithms to accelerate convergence. Mehra, Heidemann, and Aziz2 devised three algorithms that improved the convergence rate. Gupta, Bishnoi, and Kalogerakis3 proposed an acceleration scheme based on the Crowe and Nishio4 general dominant eigenvalue method that used the liquid and vapor compositions as the iteration variables. Newton-type methods have been employed as an alternative procedure to solve the isothermal flash problem. They are used to solve the nonlinear equations for equalizing the fugacities of both the vapor and liquid phases. The Newton-Raphson method was employed by Fussel and Yanosik5 in their MVNR algorithm for flash calculations. Boston and Britt6 used Broyden’s method to update the Jacobian matrix. Powell’s hybrid method was used by Neghiem, Aziz, and Li.7 The radius of convergence of the Newton-type methods is relatively small when compared to that of the successive substitution * Author to whom correspondence should be addressed.
method; hence, a good initial guess is required before convergence can be achieved. However, Gautman and Seider8 were perhaps the first to call the community’s attention to another approach. The minimization of the Gibbs free energy eventually proved to be a more suitable method for tackling the isothermal flash problem. Baker, Pierce, and Luks9 developed a self-consistent method for determining whether a state of equilibrium is stable. They calculated the Gibbs free energy surface and the tangent plane with the aid of an equation of state. If the tangent plane is above the Gibbs free energy surface then the system is unstable. In other words, the system coexists in two or more phases. Conversely, if the tangent plane is below the Gibbs free energy surface then the system is stable and exhibits only a single phase. At about the same time, Michelsen10 provided the first efficient implementation of the tangent plane criterion. Trangenstein11 discussed the difficulties of the successive substitution method near the phase boundary and convincingly demonstrated the superiority of the minimization of the Gibbs free energy approach. Nagarajan, Cullick, and Griewank12 reformulated the stability test by substituting the Gibbs free energy with the Helmholtz free energy as the thermodynamic potential to be minimized and obtaining the solution to the phase equilibrium problem by using highorder initial estimates. Sun and Seider13 employed a homotopy-continuation algorithm, whereas McDonald and Floudas 14-16 used the global optimization method and the branch-and-bound algorithm for stability analysis. In addition, McDonald and Floudas emphasized the importance of a thermodynamically stable solution with respect to perturbations in any or all of the phases in the search for the true phase equilibrium solution. To guarantee convergence to a global solution, they claimed that a global optimization approach must be used.
10.1021/ie001151d CCC: $20.00 © 2001 American Chemical Society Published on Web 11/08/2001
Ind. Eng. Chem. Res., Vol. 40, No. 25, 2001 6029
Independently, de Oliveira Jr.17 and Pan and Firoozabadi18 used the simulated annealing (SA) algorithm for direct minimization of the Gibbs free energy. They were probably the first to use the SA algorithm for the phase equilibrium calculations. Pan and Firoozabadi18 presented an excellent outline of the potentialities of the SA algorithm to obtain high-quality solutions to several phase equilibrium problems, including the solid-liquid equilibrium of heavy hydrocarbon fractions commonly encountered in petroleum reservoirs. Zhu and Xu19 applied the SA algorithm for stability analysis of the liquid-liquid equilibrium of highly nonideal systems. More recently, Zhu, Wen, and Xu20 proposed an enhanced simulated annealing algorithm to carry out the stability test of multicomponent systems at high pressures. The isothermal flash problem is, indeed, a complex minimization problem and cannot be properly solved by classical optimization techniques. In the present work, we solve the isothermal flash problem with the SA algorithm, which, in several difficult optmization problems, has proved to be a powerful tool. The simulated annealing algorithm is a stochastic method that does not make use of the Hessian matrix and is free from the ill-conditioning problems associated with systems of equations. Our methodology first tests the phase for stability and then, when necessary, minimizes the Gibbs free energy. Here, we use a change of variables to describe, from a new viewpoint, the equilibium calculations as well as the stability test problem. Both problems are formulated as linearly constrained minimization problems, and in contrast to Trangenstein’s work, the same variables are used. These are the intensive parameters. Given a vapor-liquid mixture consisting of r components, we select as the independent variables (r - 1) mole fractions of the chemical components in a given phase (vapor or liquid), previously chosen, and the phase mole fraction in the mixture. Thus, the number of independent variables is equal to the number of components present in mixture, which is, according to Gibbs’s phase rule, the least number of variables that can be used in twophase equilibrium calculations. The other parameters are considered to be dependent and are calculated from the independent variables by the imposed linear equality constraints. This choice of variables leads to two totally scaled problems with just r variables, where all parameters range between 0 and 1. Moreover, it simplifies the use of the SA algorithm, leading to an economy of computational effort. In addition to the aspects mentioned above, the proposed thermodynamic model reveals a great mathematical similarity between the stability test problem and the isothermal flash problem. This similarity is studied mathematically and then explored to find a means of minimizing the Gibbs free energy in a more efficient way. 2. Thermodynamic Modeling The formulation of the isothermal flash problem appears as the minimization of an objective nonlinear function with linear inequality constraints for the conservation of mass and the nonnegativity of phase compositions suggested by Henderson.21 For a system at a temperature T and pressure P with global molar fraction for each component, zi, i ) 1, ..., r, the vapor-liquid equilibrium problem is divided into two subproblems: the stability test problem, which will
determine the number of phases present in the equilibrium state, and the flash problem, which will determine the quantity of each component in each phase and the amount of each phase present in the mixture. Once the thermodynamic formulation is posed, the selection of the minimization algorithm for finding the Gibbs free energy global minimum (the true equilibrium solution) is the next crucial step. de Oliveira Jr.17 discussed the extreme usefulness of the simulated annealing algorithm in the phase equilibrium context, but the approach was applied here to the isothermal flash calculations only. 2.1. The Isothermal Flash Problem. We consider here only the isothermal flash, which refers to calculation of the quantities and compositions of the vapor and liquid phases making up a two-phase system in equilibrium, typically a vessel with T and P values fixed for a given feed (z1, ..., zr). We formulate the isothermal flash model using Callen’s thermodynamic parlance.22 The equilibrium state is defined as the state of minimum Gibbs free energy among the manifold of states of constant temperature and pressure. For a vapor-liquid system, the total Gibbs free energy is given by r
G)
∑ i)1
r
(v) n(v) i µi
+
(l) n(l) ∑ i µi i)1
(1)
where n(f) i , f ) v, l, represents the numbers of moles of (l) component i in the vapor n(v) i or liquid ni phase. The chemical potential (f) (f) (f) µ(f) i ) µi (n1 , ..., nr ), i ) 1, 2, ..., r and f ) v, l (2)
is a homogeneous function of degree 0, and for this reason, if we divide each extensive parameter by n(f) ) r ∑i)1 n(f) i , the total number of moles of phase f ) v, l, the same value for function µ(f) i is obtained. Thus, we can write eq 2 in terms of intensive parameters as (f) (f) (f) µ(f) i (n1 , ..., nr ) ) µi
(
)
n(f) n(f) 1 r (f) (f) , ..., ) µ(f) i (x1 , ..., xr ) (3) (f) (f) n n
(f) (f) is the molar fraction of the ith where x(f) i ) ni /n component in the liquid or vapor phase. As demonstrated in Appendix A, if n ) n(v) + n(l) is the total number of moles in the mixture, then the total molar Gibbs free energy, g ) G/n, can be written in the form
r-1 (l) g ) g(x(l) 1 , ..., xr-1, l))
(l) (v) (l) lx(l) ∑ i [(µi - µi ) - (µr i)1 r
(l) (v) µ(v) r )] + l(µr - µr ) +
ziµ(v) ∑ i i)1
(4)
(v) where µ(l) i and µi have the functional forms described by eqs A.8 and A.9, respectively, and l ) n(l)/n is the molar fraction of the liquid phase. The chemical potential of the ith component can be expressed by
(f) (f) (f) µ(f) i ) RT[lnφi + ln(Pxi )] + θi (T)
(5)
is the where R is the universal gas constant, φ(f) i fugacity coefficient of the ith component in phase f and
6030
Ind. Eng. Chem. Res., Vol. 40, No. 25, 2001
θi(T) is an integration parameter that depends only on the temperature and on the nature of the ith component; it is simply a basis at T and some fixed pressure, often 1 atm or 1 bar (see Smith and Van Ness23). If all chemical components are present in both phases, then the variables x(f) i take values strictly between 0 and 1
r-1 (l) δ(x(l) 1 , ..., xr-1, l))
r
(v) l (µ(l) r - µr ) +
x(l) ∑ i < 1 i)1
To model the liquid and vapor phases, the PengRobinson24 cubic equation of state and the classical mixing rules were utilized. As T is constant, we can ignore the parameter θ(f) i for all f ) v, l and i ) 1, 2, ..., r in the computation of the function g. This is equivalent (v) to taking θ(l) i ) θi for all i ) 1, 2, ..., r. According to the Gibbs free energy minimum principle, at the equilibrium state the values of the intensive parameters of a system are those that minimize the total molar Gibbs free energy of the system under the same specified values of temperature and pressure. Therefore, the equilibrium state of a two-phase system undergoing a process at constant temperature and pressure can be stated as follows: Given T, P, and zi, i ) 1, 2, ..., r, find the values of (l) x(l) 1 , ..., xr-1, and l that minimize r-1
...,
l) )
(v) (l) (l) µ(v) i ) µi (x1 , ..., xr-1, l)
(8)
(l) µ(l) i ) µi (z1, ..., zr) ) constant
(9)
and
r-1
(l) xr-1 ,
(l) (v) (l) (v) lx(l) ∑ i [(µi - µi ) - (µr - µr )] + i)1 r
(v) l (µ(l) r - µr ) +
ziµ(v) ∑ i i)1
(l) The function δ(x(l) 1 ,..., xr-1, l) is our modified tangent plane distance function. Therefore, if
(l) δ(x(l) 1 , ..., xr-1, l) g 0
(10)
for all x(l) i and l subject to
0 < x(l) i < 1, i ) 1, 2, ..., r - 1 r-1
x(l) ∑ i < 1 i)1
(11)
0
