The computational problems in reservoir fluid systems are mainly in the critical region and in liquid-liquid (LL), vapor-liquid-liquid (VLL), and higher-phase equilibria. The conventional methods to perform phase-equilibrium calculations with the equality of chemical potentials cannot guarantee a correct solution.
In this study, we propose a simple method to calculate the equilibrium state by direct minimization of the Gibbs free energy of the system at constant temperature and pressure. We use the simulated annealing (SA) algorithm to perform the global minimization. Estimates of key parameters of the SA algorithm are also made for phase-behavior calculations. Several examples, including (1) VL equilibria in the critical region, (2) VLL equilibria for reservoir fluid systems, (3) VLL equilibria for an H2S-containing mixture, and (4) VL-multisolid equilibria for reservoir fluids, show the reliability of the method.
Consider the multicomponent-multiphase flash at constant temperature and pressure sketched in Fig. 1. The equilibrium state (the right side of Fig. 1) consists of np phases; each Phase j consists of n1j,n2j,n3j,. . . nncj, moles. From the second law of thermodynamics, the equilibrium state is a state in which the Gibbs free energy of the system is a minimum. The minimum of Gibbs free energy is a sufficient and necessary condition for the equilibrium state.
At constant temperature and pressure (note that all calculations will be performed at this condition), the Gibbs free energy of the system in Fig. 1 can be written as
where Gj is the Gibbs free energy of Phase j, and G is the total Gibbs free energy of the system. When G is minimized with respect to nij (i=1, 2, . . ., nc; j=1, 2, . . ., np) subject to the following constraints:
material balance of Component i,
the non-negative mole number of Component i in Phase j,
The optimized values, ni(i=1, 2, . . ., nc; j=1, 2, . . ., np) are the mole numbers of the equilibrium state. The global minimization with the constraints is difficult to implement; as a consequence, direct minimization of the Gibbs free energy has not been widely applied.
The equality of chemical potentials of each species in all phases is often used to perform the phase-equilibrium calculations:
The number of equations in Eq. 4 is nc×(np-1), plus nc material-balance equations given by Eq. 2; a total of nc×np equations are provided. The mole numbers nij (i=1, 2, . . ., nc; j=1, 2, . . ., np) of the equilibrium state are determined by solving these nc×np nonlinear equations. The widely used solution methods are the successive substitution method through phase-equilibrium constants Ki (i=1, 2, . . ., nc) and direct application of the Newton method. Both approaches require an initial guess and work quite well for VL equilibria except in the near-critical region. In the critical region, the successive substitution becomes intolerably slow and the Newton method may fail when the initial guess is not close to the true solution. In LL and VLL equilibria, both methods may compute false solutions. The falseness is because Eq. 4 is only a necessary condition for an equilibrium state.1 The tangent-plane-distance (TPD) approach has been introduced to recognize the false solution.1,2 The concept of stability analysis is used to derive the TPD.
Suppose w is a given overall composition. The mathematical expression of the TPD function is
where D(u) is the distance function between the Gibbs free energy surface and its tangent plane at composition w. When D(u) is minimized with respect to ui(i=1, 2, . . ., nc) subject to
Equations 5b and 5c
the optimized value, D*, provides the stability analysis of the mixture at composition w. If D* 0, the system is absolutely stable; if D*<0, the system is unstable. The optimized composition u* is a good approximation of the incipient phase composition.
The application of TPD criterion improves the reliability of conventional-phase equilibrium by providing a guideline to judge that the mixture is absolutely stable. When unstable, a good initial composition u* strengthens the convergence of the Newton or the successive substitution methods. Unfortunately, the solution to Eq. 5 is also an optimization problem with constraints. Michelson2 has solved the problem by locating the stationary points of the TPD function. This approach needs to solve (nc-1) nonlinear equations. A good initial guess is required to avoid the trivial solution. Because not all stationary points can be found with this method, phase stability cannot always be guaranteed.3 Later, we will give an example of a CO2-crude system for which the approach of locating the stationary points misses the true solution in spite of its novelty and strengths.
Several methods have been proposed to improve the calculation of the TPD function. These include homotopy-continuation,4 branch and bound3 and differential geometry, and the theory of differential equations.5