SciELO - Scientific Electronic Library Online

vol.25 issue3High-pressure cloud point data for the system glycerol + olive oil + n-butane + AOTPerformance prediction and validation of equilibrium modeling for gasification of cashew nut shell char author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand




Related links


Brazilian Journal of Chemical Engineering

Print version ISSN 0104-6632On-line version ISSN 1678-4383

Braz. J. Chem. Eng. vol.25 no.3 São Paulo July/Sept. 2008 



Phase stability analysis of liquid-liquid equilibrium with stochastic methods



G. NagataniI; J. FerrariI; L. Cardozo FilhoII; C. C. R. S. RossiII; R. GuirardelloIII; J. Vladimir OliveiraI; M. L. CorazzaI, *

IDepartment of Food Engineering, URI, Campus de Erechim, Phone: +(55) (54) 3520-9000, Fax: +(55) (54) 3520-9090, Av. Sete de Setembro 1621, CEP: 99700-000, Erechim - RS, Brazil. E-mail address:
IIDepartment of Chemical Engineering, Maringá State University, UEM, Av. Colombo 5790, CEP: 87020-900, Maringá - PR, Brazil
IIIDPQ-FEQ, College of Chemical Engineering, State University of Campinas, UNICAMP CEP: 13081-970, Campinas - SP, Brazil




Minimization of Gibbs free energy using activity coefficient models and nonlinear equation solution techniques is commonly applied to phase stability problems. However, when conventional techniques, such as the Newton-Raphson method, are employed, serious convergence problems may arise. Due to the existence of multiple solutions, several problems can be found in modeling liquid-liquid equilibrium of multicomponent systems, which are highly dependent on the initial guess. In this work phase stability analysis of liquid-liquid equilibrium is investigated using the NRTL model. For this purpose, two distinct stochastic numerical algorithms are employed to minimize the tangent plane distance of Gibbs free energy: a subdivision algorithm that can find all roots of nonlinear equations for liquid-liquid stability analysis and the Simulated Annealing method. Results obtained in this work for the two stochastic algorithms are compared with those of the Interval Newton method from the literature. Several different binary and multicomponent systems from the literature were successfully investigated.

Keywords: Phase stability; Liquid-liquid equilibrium; Stochastic algorithms; Thermodynamic models.




The phase split phenomenon in liquid mixtures is related to thermodynamic stability of the system and comprises important information for the projection, modeling and simulation of industrial processes that involve liquid phase separation. Phase stability analysis allows determination of the exact number of coexisting phases in stable equilibrium and also provides estimation of phase compositions, affording a suitable initialization for phase equilibrium calculations.

Phase stability analysis of multicomponent mixtures is usually investigated by the Gibbs tangent plane distance criterion (Baker et al., 1982; Michelsen, 1982). For a given temperature and pressure, the necessary and sufficient condition for a phase with composition z to be stable is that the Gibss free energy surface of the mixture is not intercepted by the tangent hyper plane associated with this surface in x. To fulfill this condition, the Gibbs tangent plane distance function, TPD(x), must be nonnegative for any acceptable x. Therefore, it is possible to examine whether a phase is stable by minimization of the TPD(x) function, subject to mass balance constraints. If the tangent plane distance function at the global minimum point has a value greater than or equal to zero, then the analyzed phase is thermodynamically stable because the TPD(x) function is also nonnegative for all x vectors contained in the permissible region. If the TPD(x) function is negative at its global minimum point, then the tangent plane lies above the Gibbs free energy surface, and hence the examined phase is unstable and will be subdivided into new phases.

For the computational evaluation of this test, two distinct approaches can be employed: (i) resolution of an algebraic nonlinear set of equations that represents the stationary points of the Gibbs tangent plane distance function and (ii) direct minimization of the TPD(x) function, taking into account the constraints imposed by the mass conservation principle. In both cases, the application of conventional mathematical models presents difficulties, such as dependence on an arbitrary initialization or the possibility that the solution converges to a trial root or to a local minimum. Since multiple roots may be found, especially when liquid phases are under inspection, the first strategy requires a mathematical method able to find all roots for the set of equations. The second approach, in contrast, requires a reliable global optimization method.

The majority of numerical methods applied to phase equilibrium calculations has local convergence characteristics resulting in only local solutions and are very sensitive to initial guesses. Some work is available in the current literature regarding application of deterministic techniques for global optimization, such as homotopy continuation, branch-and-bound and Interval Newton analysis (Sun and Seider, 1995; McDonald and Floudas, 1995; Souza et al., 2006; Lima et al., 2006). However, just a few studies make use of global stochastic techniques, like Simulated Annealing and Genetic Algorithm, for phase equilibrium and phase stability calculations (Zhu and Xu, 1999; Zhu et al., 2000; Rangaiah, 2001).

In this context, the aim of this work is to provide a systematic comparison between the two mentioned strategies for phase stability analysis: Simulated Annealing (Press et al., 1992), wich is a global optimization procedure, and a subdivision method (Smiley and Chun, 2001; Corazza et al., 2007) – henceforth denominated SubDivNL – employed for solution of the stationary points on TPD of liquid mixtures modeled by the NRTL activity coefficient model (Renon and Prausnitz, 1968). Comparison of the approach used in this work with Interval Newton analysis is also provided.



The Gibbs free energy for liquid mixtures at low or moderate pressures can be modeled in terms of activity coefficients. Thus, the Gibbs tangent plane distance for a mixture with n number of components at specified temperature T and pressure P is given by (Baker et al., 1982):

where γ and are the activity coefficient of component i of tried and tested phases, respectively, and x = (x1, ..., xn) and z = (z1, ..., zn) represent the composition of tried and tested phases, respectively.

Since a phase is stable if TPD(x) > 0, the function represented in Eq. (1) must be minimized with respect to all possible compositions, i.e.:

with 0 < xi < 1 and i = 1, ..., n.

As discussed previously, this optimization problem can be solved by direct minimization of the objective function (Eq. 1) or computation of the stationary points of the TPD(x) function, wich are determined by the first-order derivatives with respect to the (n – 1) independent mole fractions. The stationary points of TPD(x) can be obtained by solving the following set of algebraic nonlinear equations, subject to the mass balance restriction:

When the activity coefficients are estimated using the NRTL model, this set of equations may have multiple roots in addition to the trivial solution x = z.

Simulated Annealing

Simulated Annealing (SA) is a stochastic global optimization technique suitable for detecting a global minimum hidden among many local minima (Press et al., 1992). This method explores the analogy with the way that, when cooled slowly, a metal forms a crystalline structure with minimum energy. Metropolis et al. (1953) proposed that the thermal equilibrium of an annealing process could be simulated by the Boltzmann probability distribution.

The Simulated Annealing algorithm is based on the Monte Carlo method and makes use of an objective function φ(x) to replace energy E. The vector x is defined in the n-dimensional space of the optimization variables and represents a particular configuration of the system. A parameter TA, analogous to the annealing temperature, is employed for evaluation of the Boltzmann factor. In accordance with Faber et al. (2005), a new vector, xk+1, is randomly generated from the earlier state xk. If φ(xk+1) < φ(xk), then the new vector is unconditionally accepted, i.e., xk+1 = xk = xopt; otherwise, the acceptance probability is determined by the Metropolis criterion.

The progressive decrease in annealing temperature reduces the probability that a more energetic state be selected. As suggested by Faber et al. (2005), the annealing schedule can be set by for j = 1, ..., ζ, in which 0 < α < 1. However, just as excessively fast cooling may produce a structure with a higher energetic content, an inappropriate change in the annealing temperature can result in convergence to a local minimum. Furthermore, it is also essential to start from a conveniently high temperature in order to allow a more detailed exploration of the solution space of the problem.

Over the past years the SA algorithm has been applied to several types of phase equilibrium calculation problems involving phase stability tests and flash equilibrium calculations by minimizing the Gibbs free energy (Zhu and Xu, 1999; Zhu et al., 2000; Rangaiah, 2001; Teh and Rangaiah, 2003; Corazza et al., 2004; Souza et al., 2004; Bonilla-Petriciolet, 2006) and parameter estimation of thermodynamics models (Henderson et al. 2001; Moura et al., 2005; Carvalho Jr., 2006; Franceschi et al., 2006). Preliminary studies indicated that the best values for phase stability calculations are = 1000.0, α = 0.9 and z = 10000. In this work, the SA algorithm used was that presented in Press et al. (1992), where the probability acceptance criterion used is Pr(φ(xK+1)) > xr, Pr is the Boltzmann probability distribution and xr is a random number [0,1].

Subdivision Algorithm

According to Smiley and Chun (2001) the subdivision algorithm consists in a given set of algebraic nonlinear equations F(x) and a primitive interval (hyper rectangle) R Rd, where d is the problem dimension, and the aim of a subdivision algorithm is to determine all values x* = {x Rd|F(x) = 0} through systematic subdivisions and trials in R. The first partition of R to congruent subrectangles is by division of all coordinates (variables) into two equal parts so as to obtain 2d subrectangles Rij from one rectangle in the Rd space (the index j identifies the subrectangles produced at a subdivision level i).

For each subdivision level, the subintervals obtained are tested for the existence of roots in order that only those containing one or more solutions are maintained. The selection criterion is based on calculation of the Euclidean norm of the set of equations:


then the subrectangle Rij is retained or else, it is discarded (y is the vector of the randomly sampled points), and ak and bk are the lower and upper bonds of unknowns.

After a finite number of subdivisions, a conventional method (e.g. the Newton-Raphson method) is employed to determine the roots, taking as initial values the middle points of the last subintervals. Since the solutions are confined to the retained subintervals, the convergence of the adopted method becomes more efficient and reliable (Smiley and Chun, 2001).

A recent study, a subdivision algorithm for phase equilibrium calculations of systems at high-pressures modeled with an equation of state was proposed. The results obtained were quite satisfactory (Corazza et al., 2007). It may be interesting to note that, to our knowledge, an application of a subdivision algorithm to liquid-liquid equilibrium calculations is not available in the literature. In this work, from preliminary tests, the appropriate maximum number of subdivisions for phase stability tests was determined to be six.



To validate the algorithms used in this work (Simulated Annealing and SubDivNL subdivision method), some results presented in the literature for liquid-liquid stability calculations with the NRTL model were selected. The investigated systems include binary, ternary and multicomponent mixtures. All the calculations were performed in a PC Pentium 4 2.66GHz with 512MB RAM.

Case 1: citric acid (1) + butan-2-ol (2) system

In Table 1 the results of phase stability tests for the citric acid (1) + butan-2-ol (2) binary system at 25ºC and different overall compositions are presented. The binary parameters of the NRTL model were taken from Gecegormez and Demirel (2005), which were presented by Lintomen et al. (2001). One can see in Table 1 that, for this system, solutions obtained in this work are in perfect agreement with the roots found by Gecegormez and Demirel (2005), which used an interval analysis algorithm. It can also be observed that the SubDivNL algorithm had lower computational effort in terms of CPU time that the SA algorithm (our results) and was shown to be much less time-consuming than the Interval Newton (IN/GB) algorithm reported on Gecegormez and Demirel (2005).

Case 2: n-pentanol (1) + 2,2-dimethylbutane (2) system

Results of phase stability analysis obtained in this work for the n-pentanol (1) + 2,2-dimethylbutane (2) system at 25 ºC and several global compositions are presented in Table 2. Results from the SA and SubDivNL algorithms are compared with solutions reported by Gecegormez and Demirel (2005). The NRTL parameters, which refer to a VLE set of parameters previously reported in the literature, were taken from Demirel and Gecegormez (1991).

It can be observed in Table 2 that for the two overall compositions, z = (0.05, 0.95) and z = (0.20, 0.80), results obtained in this work are different from those presented by Gecegormez and Demirel (2005). The TPD curves for mixtures containing respectively 0.05 and 0.20 in mole fraction of n-pentanol are represented in Figures 1 and 2. As clearly shown in these figures, the nontrivial solutions reported by Gecegormez and Demirel (2005) are not true solutions for the phase stability test, so the system has only one phase in both cases.





Moreover, the solutions presented in the literature are not stationary points on the TPD surface. In fact, the Euclidean norms of the set of equations calculated for such roots are far from zero: for z = (0.05, 0.95) and x = (0.1294, 0.8706), the value found is 0.0229, and for z = (0.20, 0.80) and x = (6.5086 × 10-2, 0.9349), the result is 0.0292.

For all overall compositions tested for this system, our calculations only indicated the existence of a one-phase, homogeneous system, which just corroborates phase equilibrium data reported by Sayegh and Ratcliff (1976).

Nevertheless, one should call attention to the fact that the use of a set of parameters coming from VLE to LLE calculations is not recommended (Poling et al., 2000). In this work, however, a parameter set from VLE (Demirel and Gecegormez, 1991) for LLE computations is applied just to enable numerical comparison with the results of Gecegormez and Demirel (2005).

Case 3: water (1) + butyl glycol (2) system

The NRTL parameters for this system, which are presented in the literature (Demirel and Gecegormez, 1991) and refer to a set of parameters from VLE, were taken from Gecegormez and Demirel (2005). Solutions of the phase stability test in relation to LLE for the water (1) + butyl glycol (2) system at 5ºC and several global compositions are presented in Table 3. As depicted in Figures 3 and 4 the nontrivial roots reported by Gecegormez and Demirel (2005) for the water overall molar fractions of 0.05 and 0.10, respectively, are not true roots. According to the results from the SA and SubDivNL algorithms presented in this work, this system has only one stable phase at these compositions.





Case 4: water (1) + citric acid (2)+ 2-butanol (3) system

The LLE binary parameters of the NRTL model for this system were used as presented by Gecegormez and Demirel (2005), originally reported by Lintomen et al. (2001). Table 4 contains the results of phase stability analysis for this ternary system at 25ºC and different overall compositions. Solutions found in this work are compared with the roots presented by Gecegormez and Demirel (2005).

As shown in Table 4, some divergences between our results and those presented in the literature can be verified. For the global composition of z = (0.10, 0.05, 0.85), the composition with the lowest TPD value found in this work was x = (8.4640 × 10-2, 0.4743, 0.4411) with TPD = -0.1609, while this root is not reported in the literature. Also, for the global composition z = (0.05, 0.10, 0.85) the SubDivNL and SA algorithms found an unstable composition, not reported in the literature. It can also be verified from Table 4 that for the last third overall composition, a second root, found in the literature, was not found by the SA and SubDivNL algorithms.

Another discrepancy was verified for the two overall compositions: a) z = (0.25, 0.05, 0.70) and b) z = (0.30, 0.05, 0.65) that correspond to mass fractions of a) z = (0.0682, 0.1456, 0.7862) and b) z = (0.0855, 0.1520, 0.7625), respectively. While unstable roots from the phase stability test are presented in the literature, indicating the formation of two liquid phases, application of the SA and SubDivNL algorithms did not result in the occurrence of separate liquid phases (see Table 4). In order to elucidate these results, the ternary experimental diagram for this system is presented in Figure 5. The experimental data were obtained from the literature (Lintomen et al., 2001), and as one can see, the two feed points (a) and (b) lie in a homogeneous region, which is in complete agreement with the results found in this work.



Case 5: acetonitrile (1) + benzene (2) + n-heptane (3) system

Table 5 contains the results of phase stability analysis for this ternary system at 45ºC and different overall compositions. Solutions found in this work are compared with the roots presented by Gecegormez and Demirel (2005). As shown in this table, application of the subdivision algorithm indicated the existence of a third root for all compositions analyzed. This additional root however is not presented in the work of Gecegormez and Demirel (2005). For this problem, the Simulated Annealing algorithm had a better performance in terms of CPU time than the subdivision method.

Case 6: n-propanol (1) + n-butanol (2) + benzene (3) + water (4) system

Solutions of the phase stability test for this system are presented in Table 6. The binary interaction parameters for the NRTL model were taken from Tessier et al. (2000). Aiming at reducing the computational effort, the maximum number of subdivisions was set to four. It can be verified that both algorithms employed in this work were able to find the same solutions as those presented in the literature (Tessier et al., 2000). The CPU time values listed in this table confirm the advantageous performance of the Simulated Annealing algorithm for systems with larger number of components.

Case 7: n-propanol (1) + n-butanol (2) + benzene (3) + ethanol (4) + water (5) system

Table 7 depicts the results of the thermodynamic stability test for this multicomponent system. For this system the NRTL parameters were also taken from Tessier et al. (2000), whose results are compared in this work. The maximum number of subdivisions was again set to four. Although the Simulated Annealing algorithm had a better performance in terms of CPU time than the SubDivNL, solutions obtained with both mathematical methods employed in this work are in agreement with those reported by Tessier et al. (2000), who used the Interval Newton method.



In this work thermodynamic phase stability with two stochastic global algorithms was examined for binary, ternary and multicomponent liquid mixtures using the NRTL activity coefficient model. Results of application of the subdivision algorithm to liquid-liquid equilibrium calculations show that it was able to find all roots of phase stability tests for the systems investigated with relatively short CPU times for binary and ternary systems. Simulated Annealing was shown to be reliable and robust for systems with larger numbers of components. In a general sense, the two algorithms are relatively simple, provide reliable solutions for phase stability analysis, are easy to implement and require relatively low computational effort.




problem dimension



nonlinear equation



system of nonlinear equations



multidimensional rectangle



lower coordinate of the "child" rectangles at "k" dimension



upper coordinate of the "child" rectangles at "k" dimension



lower coordinate of the "parent" rectangles at "k" dimension



upper coordinate of the "parent" rectangles at "k" dimension



nonrandom two liquid



subdivision level



vector of random point in Rij and vapor phase mole fraction



middle point of the rectangle "j" at the subdivision level "i"'



vector of liquid phase mole fraction



vector of global mole composition






initial annealing temperature



sample number in Rij



maximum subdivision number



maximum coverage number



tangent plane distance


Greek Letters


coordinate of the middle point of "parent" rectangle at "k" dimension



activity coefficient of component "i"



subdivision parameter of the rectangle "j" at the subdivision level "i"



matrix of zero and one elements




pure component



number of components in the mixture and nth components of the system




The authors thank CENPES/PETROBRAS, FAPERGS and CNPq for the scholarships and financial support of this work.



Baker, L. E., Pierce, A. C. and Luks, K. D., Soc. Pet. Eng. J., 22, 731 (1982).         [ Links ]

Bonilla-Petriciolet, A., Vazquez-Roman, R., Iglesias-Silva, G. A. and Hall, K. R., Performance of stochastic global optimization in the calculation of phase stability analyses for nonreactive and reactive mixtures, Ind. Eng. Chem. Res., 45, 4764 (2006).         [ Links ]

Carvalho Jr., R. N., Corazza, M. L., Cardozo-Filho, L. and Meireles, M. A. A., Phase equilibrium for (camphor+CO2), (camphor+propane), and (camphor +CO2+propane), J. Chem. Eng. Data, 51, 997 (2006).         [ Links ]

Corazza, M. L., Corazza, F. C., Cardozo Filho, L. and Dariva, C., A subdivision algorithm for phase equilibrium calculations at high pressures, Braz. J. Chem. Eng, In press (2007).         [ Links ]

Corazza, M. L., Filho, L. C., Oliveira, J. V. and Dariva, C., A robust strategy for SLV equilibrium calculations at high pressures, Fluid Phase Equilibria, 221, 113 (2004).         [ Links ]

Demirel, Y. and Gecegormez, H., Simultaneous representation of excess enthalpy and vapor-liquid equilibrium data by NRTL and UNIQUAC models, Fluid Phase Equilibria, 65, 111 (1991).         [ Links ]

Faber, R., Jockenhövel, T. and Tsatsaronis, G., Dynamic optimization with Simulated Annealing, Computers & Chem. Eng., 29, 273 (2005).         [ Links ]

Franceschi, E., Kunita, M. H., Rubira, A. F., Muniz, E. C., Corazza, M. L., Oliveira, J. V. and Dariva, C., Phase behavior of binary and ternary systems involving carbon dioxide, propane, and glycidyl methacrylate at high pressure, J. Chem. Eng. Data, 51, 686 (2006).         [ Links ]

Gecegormez, H. and Demirel, Y., Phase stability analysis using Interval Newton method with NRTL model, Fluid Phase Equilibria, 237, 48 (2005).         [ Links ]

Henderson, N., Oliveira Jr., J. R., Amaral Souto, H. P. and Marques, R.P., Modeling and analysis of the isothermal flash problem and its calculation with the Simulated Annealing algorithm, Ind. Eng. Chem. Res., 40, 6028 (2001).         [ Links ]

Lima, E. R. A., Tavares, F. W., Biscaia, E. C., Jr., Utilização de Algoritmos Heurísticos de Otimização para Cálculo de Equilíbrio de Fases. In: XVI Congresso Brasileiro de Engenharia Química – COBEQ, Santos (2006).         [ Links ]

Lintomen, L., Pinto, R. T. P., Batista, E., Meirelles, A. J. A. and Maciel, M. R. W., Liquid-liquid equilibrium of the water+citric acid+short chain alcohol+tricaprylin system at 298.15 K, J. Chem. Eng. Data, 46, 546 (2001).         [ Links ]

McDonald, C. and Floudas, C. A., Global optimization for the phase stability problem, AIChE J., 41, 1798 (1995).         [ Links ]

Metropolis, N., Rosenbluth, A. W. and Rosenbluth, M. N., Equation of state calculations by fast computing machines, J. Chem. Phys., 21, 1087 (1953).         [ Links ]

Michelsen, M., The isothermal flash problem. Part I. Stability, Fluid Phase Equilibria, 9, 1 (1982).         [ Links ]

Moura, L. S., Corazza, M. L., Cardozo-Filho, L. and Meireles, M. A. A., Phase equilibrium measurements for the system fennel (Foeniculum vulgare) extract + CO2, J. Chem. Eng. Data, 50, 1657 (2005).         [ Links ]

Poling, B.E., Prausnitz, J.M. and O'Connell, J.P., The Properties of Gases and Liquids, 5th ed., McGraw Hill, New York (2000).         [ Links ]

Press, W. H., Flannery, B. P., Teukolsky, S. A. and Veterlling, W. T., Numerical Recipes in FORTRAN - The Art of Scientific Computing, 2nd ed., Cambridge University Press, New York (1992).         [ Links ]

Rangaiah, G. P., Evaluation of genetic algorithms and Simulated Annealing for phase equilibrium and stability problems, Fluid Phase Equilibria, 187-188, 83 (2001).         [ Links ]

Renon, H. and Praunitz, J. M., Local compositions in thermodynamic excess functions for liquid mixtures, AIChE J., 14, 135 (1968).         [ Links ]

Sayegh, S. G. and Ratcliff, G. A., Excess Gibbs energies of binary systems of isopentanol and n-pentanol with hexane isomers at 25 ºC: Measurement and prediction by analytical group solution model, J. Chem. Eng. Data, 21, 1, 71 (1976).         [ Links ]

Smiley, M. W. and Chun, C. J., An algorithm for finding all solutions of a nonlinear system, Computational and Applied Math., 137, 293 (2001).         [ Links ]

Souza, A. T., Cardozo-Filho, L. and Wolff, F., Application of interval analysis for Gibbs and Helmholtz free energy global minimization in phase stability analysis, Braz. J. Chem. Eng., 23, 117 (2006).         [ Links ]

Souza, A. T., Corazza, M. L., Cardozo-Filho, L., Guirardello, R. and Meireles, M. A. A., Phase equilibrium measurements for the system clove (Eugenia caryophyllus) oil + CO2, J. Chem. Eng. Data, 49, 352 (2004).         [ Links ]

Sun, A. C. and Seider, W. D., Homotopy-continuation method for stability analysis in the global minimization of the Gibbs free energy, Fluid Phase Equilibria, 103, 203 (1995).         [ Links ]

Teh, Y. S. and Rangaiah, G. P., Tabu search for global optimization of continuous functions with application to phase equilibrium calculations, Computers & Chem. Eng., 27, 1665 (2003).         [ Links ]

Tessier, S. R., Brennecke, J. F. and Stadtherr, M. A., Reliable phase stability analysis for excess Gibbs energy models, Chem. Eng. Science, 55, 1785 (2000).         [ Links ]

Zhu, Y. and Xu, Z., A reliable prediction of the global phase stability for liquid-liquid equilibrium through the Simulated Annealing algorithm: Application to NRTL and UNIQUAC equations, Fluid Phase Equilibria., 154, 55 (1999).         [ Links ]

Zhu, Y., Xu, Z. and Wen, H., Global stability analysis and phase equilibrium calculations at high pressures using the enhanced Simulated Annealing algorithm, Chem. Eng. Science, 55, 3451 (2000).         [ Links ]



(Received: February 19, 2007 ; Accepted: March 5, 2008)



* To whom correspondence should be addressed

Creative Commons License All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License