Modified Bounded Homotopies in the Solving of ... - ACS Publications

ARTICLE pubs.acs.org/IECR

Modified Bounded Homotopies in the Solving of Phase Stability Problems for Liquid
Liquid Phase-Splitting Calculations Jani Kangas,* Ilkka Malinen, and Juha Tanskanen Department of Process and Environmental Engineering, University of Oulu, P.O. Box 4300, FI-90014 University of Oulu Finland ABSTRACT: In this study, the modified bounded homotopies presented by Malinen and Tanskanen (Malinen, I.; Tanskanen, J. Modified bounded homotopies to enable a narrow bounding zone. Chem. Eng. Sci. 2008, 63, 3419) are investigated in order to solve phase stability analysis problems in liquid
liquid equilibrium cases for phase-splitting calculations. The tangent-plane distance criterion is used to analyze the phase stability. The emphasis is on approaching the first root on the homotopy path. According to the observations, the bounding of the homotopy path with respect to the problem variables aids in the solving of a phase stability analysis problem. The main shortcomings of the modified bounded Newton, affine, and fixed-point homotopies are the starting point isolas, the existence of unfeasible solutions, and the convergence to only certain roots, respectively. The attraction domain of the global minimum was observed to be the largest with the fixed-point homotopy. Sequential usage of the fixed-point homotopy and starting points near pure components results in a robust algorithm for phase stability analysis and gives feasible initial estimates for the solving of phase splitting.

1. INTRODUCTION An accurate mathematical model of a separation process unit demands detailed thermodynamic data obtained from experiments. These experimental data are crucial for determining the parameters for a thermodynamic model. In addition to an accurate thermodynamic model, a robust analysis procedure for determination of the phase stability is essential in process simulation. Inaccuracies in the phase equilibrium prediction stage may result in process design failures.1 One possible source of inaccuracies is the failure of the solving method used in phase stability analysis. Evaluation of the phase stability is difficult chiefly in a separation process that involves two or more liquid phases. This fact originates from the properties of the local composition models, such as the NRTL and UNIQUAC models, frequently used to describe the phase equilibrium in the system. The mathematical forms of the models have in common the possibility of having multiple roots. The multiplicity property is generally demanded to describe the behavior of the fluid accurately and flexibly. In particular, the flexibility of the model may result in a false number and composition of phases. This makes a thorough analysis of the computed results obligatory. To outline the importance of accurate phase equilibrium solution, Bekiaris et al.,2 among others, have studied the effects of inaccuracies in phase equilibrium models for the simulation and design of heterogeneous azeotropic distillation systems. They have observed that the inaccuracies may cause, for example, multiple steady states in the systems, and if there are inaccuracies in the representation, it is very likely that the multiplicities will be affected qualitatively and quantitatively. In addition, if we are searching for azeotropes in a mixture prone to forming multiple liquid phases, we have to apply a rigorous stability test to make our search thorough and complete.3 It is also worth noting that the local composition models have in common the possibility of modeling, i.e., of predicting, the behavior of multicomponent mixtures with the binary interaction parameters obtained from pure component and binary measurement r 2011 American Chemical Society

data. However, in certain cases, they tend to predict more liquid phases than actually exist in the real system [see, e.g., the discussion on spurious solutions in liquid
liquid equilibrium (LLE) in Mitsos et al.4]. McDonald and Floudas5 and Xu et al.,6 among others, have presented examples of unfortunate errors in the selection of appropriate thermodynamic model parameters and the solving of phase equilibrium compositions. Thus, accurate analysis of the phase stability using robust solving methods and careful examination of the measurement results are vital for the efficient design of separation processes with multiple liquid phases. During the last couple of decades, the research work performed in the field of the robust solving of phase equilibrium compositions has been extensive and varied, in terms of both problem formulation approaches and solving methods (see sections 2.1 and 2.2). Therefore, this paper first gives an overview of the topic and a brief presentation of the theory behind phase stability analysis. The modified bounded homotopies proposed by Malinen and Tanskanen7 are then presented as serviceable solving methods for phase stability analysis. One aim of this paper is to seek ways to improve the performance of the modified bounded homotopy methods as robust and multipurpose solving methods. The properties of modified bounded homotopies are investigated using several multicomponent LLE cases and a LLLE case, in which the well-known excess Gibbs energy models (NRTL and UNIQUAC) are used. In these cases, we examine the convergence of various homotopy methods from different starting points to the first root on the homotopy path in particular. The attraction domain study presented in this paper has not been utilized previously in the literature in conjunction with homotopy continuation methods for phase stability analysis problems. However, a few attraction domain studies have been performed for Received: September 16, 2010 Accepted: April 19, 2011 Revised: March 28, 2011 Published: April 19, 2011 7003

dx.doi.org/10.1021/ie101907h | Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research some general nonlinear equation (NLE) sets; see Wayburn and Seader.8 The attraction domain is defined in this paper as the domain of starting points in which the solving method converges to a certain root first. Finally, a proposal for solving of the phase stability problem for subsequent phase-splitting calculations is presented based on the attraction domain study that was carried out.

2. OVERVIEW OF THE SOLVING OF PHASE EQUILIBRIUM COMPOSITIONS A common multiphase problem consists of finding the correct number and types of phases and their corresponding equilibrium compositions such that the Gibbs free energy of the system is at the global minimum. A considerable number of methods and algorithms have been formulated to solve the phase equilibrium problem. Generally, the basis in the algorithms is to solve a sequence of subproblems. The algorithms can be separated into two groups: equation-solving and Gibbs free energy minimization approaches.9,10 The equation-solving approaches solve a set of NLEs, based on the component mass balances and phase equilibrium relations. The equations fulfill the necessary conditions of equality of the chemical potential of every component in each phase, but they do not ensure that the sufficient condition of phase equilibrium is satisfied.10,11 The second approach of minimizing the Gibbs free energy yields the sufficient condition of phase equilibrium.12 Thus, the Gibbs free-energy minimization approach is more preferable than the equation-solving approach in the robust evaluation of the phase equilibrium conditions. The core of the Gibbs freeenergy minimization approaches is to evaluate the multicomponent and multiphase equilibrium based on the global minimum of the Gibbs free energy of the system.11 There are fundamentally two routes that can be followed to obtain or verify that an equilibrium solution corresponds to the global minimum of the Gibbs free energy. The two optimization problems associated with these approaches are the direct minimization of the Gibbs free energy and the minimization of the tangent-plane distance function (TPDF). 2.1. Gibbs Free-Energy Minimization. A variety of optimization methods can be used to minimize the Gibbs free energy directly. However, the main issue to be tackled is that the number and types of phases present at equilibrium are not known beforehand. Therefore, the application of optimization methods to the solution of the phase equilibrium problem may bring about severe computational difficulties. The difficulties can be bypassed by fixing the number of phases and using nonlinear programming to minimize the Gibbs free energy of the system among the phases. In this approach, it is worth noting that the number of phases has to be sufficiently high at the start of the solution in order to reach the correct solution robustly. This, in turn, may be inefficient with respect to the solving time. The other alternative is to add consecutively additional phases in the computations and test the stability of the resulting phase configuration by minimizing the TPDF. Minimization of the TPDF has received a lot of attention in the literature in recent years (e.g., Lucia et al.,13 Sun and Seider,14 Tessier et al.,15 Khaleghi and Jalali,16 Zhu and Xu,17,18 and Bonilla-Petriciolet et al.19). In spite of this recent activity in the field, Baker et al.11 were the first to describe phase stability analysis in terms of the tangent-plane criterion with a practical numerical implementation. The seeking of minima in the tangent-plane distance is usually performed by solving a system of NLEs for the stationary points, as formulated by Michelsen.20 Nagarajan et al.21

ARTICLE

presented a reformulation of Michelsen’s method of phase stability analysis, which employed component molar densities as the primary variables in place of mole numbers. They stated that the reformulation improves the numerical stability of the Newton iterative scheme. Recently, Mitsos and Barton22 proposed an alternative interpretation of phase stability analysis via Lagrangian duality, in which differentiability of the Gibbs free-energy surface is not required. An alternative to minimization of the TPDF is to maximize the area (or volume) between the Gibbs free-energy surface and the tangent plane. This approach has been proposed by Eubank et al.23 The “area method” of Eubank et al.23 relies on minimization of the Gibbs free energy of the system by integration, rather than differentiation, of the Gibbs energy curve to provide the maximum area rather than its derivative, the tangent line. Later Balogh et al.24 generalized the area method from a two-component, one- and two-phase system analysis to hypersurfaces that are characteristic for the Gibbs free energy of multicomponent, multiphase systems. The generalization of the area method has not yet been implemented in the literature, and as Balogh et al.24 state, the generalization “will stimulate fresh research”. As a whole, the field of phase equilibrium analysis is an extensive one. Readers are urged to familiarize themselves with a thorough review of Wakeham and Stateva9 on the numerical solution of isothermal, isobaric phase equilibrium problems. The solution of the phase stability analysis problems requires robust solution methods and is thus assessed in depth in the following section. 2.2. Solution Methods for Phase Stability Analysis. The phase stability analysis problem has been tackled using different solution methods. These methods can be divided into local and global methods. Local methods include reformulations of Newton’s method. Local methods have been found to be insufficient to tackle the phase stability analysis problem with a good degree of certainty. Thus, the majority of studies use global methods. There are two classes of global methods: stochastic methods and global deterministic optimization methods. Global methods guarantee at least theoretically that the optimum is obtained irrespective of the starting point. 2.2.1. Stochastic Methods. The interest in stochastic methods in the literature is wide. Balogh et al.25 used a reformulation of the stochastic solving method presented by Boender et al.26 to predict phase stability. They applied the method to cubic equations of state to study vapor
liquid equilibrium (VLE) cases. Zhu and Xu18 presented a simulated annealing algorithm to evaluate global phase stability in LLE cases. The algorithm has been refined by Zhu et al.27 into an enhanced simulated annealing algorithm to make the algorithm capable of controlling the “cooling schedules” themselves. More recently, Bonilla-Petriciolet et al.19 compared the performances of different stochastic global optimization methods in the case of both nonreactive and reactive mixtures. They observed that the most reliable method is the simulated annealing method. A similar evaluation of solution methods was performed by Rangaiah,10 who compared the performance of the simulated annealing method with respect to that of genetic algorithms in phase equilibrium calculation cases, including phase stability analysis cases. Contrary to the observations of Bonilla-Petriciolet et al.,19 Rangaiah10 found that the genetic algorithm is more efficient and reliable in phase equilibrium cases than the simulated annealing method. The genetic algorithm is able to locate the global minimum of the TPDF reliably in all cases.10 Srinivas and Rangaiah28 assessed the performance of the tabu search compared to the differential evolution method. They used 7004

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research a local minimizer in the vicinity of the global minimum to refine the solution. Both methods are highly reliable according to Srinivas and Rangaiah,28 but the computational efficiency of the differential evolution is better than that in the tabu search method. 2.2.2. Deterministic Global Optimization Methods. Deterministic global optimization methods have been used successfully to solve phase stability analysis problems. The methods include the interval Newton methods, branch-and-bound methods, homotopy continuation methods, the tunneling method, and the terrain method. According to Hua et al.,29 the interval Newton methods, when combined with a generalized bisection approach, provide the power to find with confidence all of the solutions of a system of NLEs. In addition, the interval Newton methods are able to find with total reliability the global minimum of a nonlinear function, provided only that the upper and lower bounds are available for all variables.29 Hua et al.29,30 applied the method to cubic equations of state for VLE cases. More recently, Tessier et al.15 used the method to predict the phase stability with Gibbs excess energy models. The study of Tessier et al.15 examined the appearance of multiple liquid phases instead of the VLE. Xu et al.6 continued the research with asymmetric models; i.e., the vapor phase is modeled with an equation of state and the liquid phase with a Gibbs excess energy model. They defined the new concept of “pseudotangent plane distance” in order to avoid difficulties with the nondifferentiable objective function, i.e., minimization of separate tangent plane distances for both the vapor and liquid phases. Additionally, Gecegormez and Demirel31 reported success in the use of an interval Newton method to predict the phase stability of liquid phases. The cases studied by Gecegormez and Demirel31 include several binary and ternary mixtures with the NRTL model. The branch-and-bound methods are also seen as a feasible alternative for solving phase-stability-related problems. Zhu and Xu17 introduced, in addition to a simulated annealing algorithm, a new branch-and-bound algorithm for phase stability analysis of excess Gibbs energy models. Wasylkiewicz et al.32 applied the ideas of differential geometry to their approach of finding stationary points of the TPDF in several LLE cases. They also used the knowledge of the topology of the surface (“topological criterion”) to specify a situation in which all of the stationary points have been found. Exploration of the function domain was performed by traveling the ridges and valleys of the surface systematically. McDonald and Floudas33 presented a global optimization approach to phase and chemical equilibrium problems. In the approach, they transformed the phase stability problem into a biconvex problem as the NRTL model was used. Thereupon, they solved the NRTL-related problem with the GOP algorithm of Floudas and Visweswaran.34 Furthermore, McDonald and Floudas5 converted the form of the UNIQUAC model later on into the difference of two convex functions, where the concave portion is separable. The phase stability problem with the UNIQUAC model was solved with a branch-and-bound algorithm. Similar to the UNIQUAC model, McDonald and Floudas35 showed that the UNIFAC, modified Wilson, and ASOG models can be transformed to a form for which the branch-and-bound algorithm is applicable. On the basis of previous approaches, McDonald and Floudas36 generalized the approaches into an algorithm that is able to solve both the phase and chemical equilibrium problems for a large group of Gibbs free-energy models. Nichita et al.37 applied the tunneling method for solving phase stability problems. This method has two stages. In the first stage, a local bounded optimization method is used to minimize the

ARTICLE

objective function. In the second stage, either global optimality is ascertained or a feasible initial estimate for a new minimization is generated. Nichita and Gomez38 later refined their formulation of the tunneling method. Lucia et al.13 applied the terrain methods to find all of the stationary points of the TPDF in the phase stability problem. They exploited the knowledge of stationary points, singular points, changes in convexity, and terrain paths obtained during the phase stability calculations to reliably initialize multiphase equilibrium calculations for any investigated mixture composition. 2.3. Homotopy Continuation Methods in Phase Equilibrium Examinations. Homotopy continuation methods are generally used in the field of phase behavior analysis to solve nonlinear algebraic equation sets. The literature regarding homotopy continuation methods in phase equilibrium studies can be divided into the use of either problem-dependent or problem-independent homotopies. 2.3.1. Problem-Dependent Homotopies. The emphasis in the literature is to take into account the application properties. Thus, the homotopy is reformed frequently into a problem-dependent form. The problem dependency of a homotopy is defined here to include one or several of the options below: • Usage of a problem variable in place of an artificial homotopy parameter. • Change in the properties of the problem along the homotopy path; e.g., the nonideality of the mixture increases from Raoult’s law behavior to actual nonideal behavior as the value of the homotopy parameter is increased. • Usage of the homotopy parameter only for certain equations of the problem. Among others, Fidkowski et al.39 computed homogeneous azeotropes in multicomponent mixtures using a Raoult’s lawbased starting point for component mole fractions in vapor and liquid phases. The solution procedure is problem-dependent because the increase in the homotopy parameter value results in an increase of the system nonideality. Later Eckert and Kubicek40 generalized the approach of Fidkowski et al.39 to the calculation of heterogeneous mixtures. Wasylkiewicz et al.3 applied the Fidkowski et al.39 approach to find also heterogeneous azeotropes. They used the rigorous stability analysis presented in their earlier paper32 to evaluate the LLE. Tolsma and Barton12,41 studied the same application of calculating heteroazeotropes, yet with different solution procedures. More recently, Wasylkiewicz et al.42 analyzed the pressure sensitivity of azeotropes with the Fidkowski et al.39 approach. Wang et al.43 calculated the critical loci of binary mixtures with a homotopy continuation method. Aslam and Sunol44 evaluated the sensitivity of the azeotropic states to activity coefficient model parameters and system variables using either one of the parameters or one of the variables as the homotopy parameter. The study of Aslam and Sunol44 is a followup of their earlier research45 of computing binary homogeneous azeotropes at high pressures through an equation of state. 2.3.2. Problem-Independent Homotopies. Problem-independent homotopies have also been used in place of problemdependent homotopies. In this approach, the phase stability analysis problem was investigated with the help of an artificial homotopy parameter. However, some heuristics regarding the problem properties were utilized usually to ensure that a feasible starting point was chosen (see, e.g., Sun and Seider14 and Bausa and Marquardt46). 7005

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research Bausa and Marquardt46 presented a hybrid method in which a priori knowledge of the phase diagram properties was used to determine the heterogeneous (VLLE or LLE) regions. Both the phase stability stage and the flash calculations are performed with a homotopy continuation method. Steyer et al.47 proposed improvements to the algorithm presented by Bausa and Marquardt46 by reducing the number of variables. They also presented a ratebased approach to liquid
liquid phase-splitting calculations. According to Steyer et al.,47 the computational efficiency increases significantly in the cases studied. Sun and Seider14 studied both symmetric and asymmetric phase equilibrium models and restricted evaluation to the Newton homotopy. Sun and Seider14 obtained all of the stationary points of the TPDF in the cases that they studied. The algorithm included a choice of starting points based on the suggestions of Michelsen.20 Jalali-Farahani and Seader48 applied a homotopy continuation method for multiphase, reacting systems. They used the Newton homotopy to calculate the chemical and phase equilibrium and the fixed-point homotopy to solve the phase stability problem. Jalali and Seader49 reported that the Newton homotopy was preferred because of its self-scaling property. They also assessed the fixed-point homotopy properties and stated that the fixedpoint homotopy may have a scaling problem, which is an important issue in the computations. They also claimed that, because of the high nonlinearity of the functions, calculation of the function derivatives is problematic with the (scale-invariant) affine homotopy. Khaleghi and Jalali16 relaxed calculations to include the complex space and used bifurcation theory to find more solutions through the tracking of branches originating from bifurcation points. They applied the Newton homotopy to solve the phase equilibrium and the fixed-point homotopy to test the phase stability, similarly to Jalali-Farahani and Seader.48 According to Khaleghi and Jalali,16 the fixed-point homotopy was chosen because of the sensitivity of the homotopy to the selection of an initial guess. The affine homotopy function was not included in the Khaleghi and Jalali16 study. More recently, Jalali et al.50 stated that relaxation to the complex space was not applicable unless the equations were converted to include complex variables.

3. PHASE STABILITY ANALYSIS FOR PHASE-SPLITTING CALCULATIONS 3.1. Problem Formulation. The stability of a phase can be determined by using the Gibbs tangent-plane criterion.20 The theory states that a phase with a composition z in temperature T and pressure P is unstable if its molar Gibbs mixing energy surface is below the tangent plane to the surface at any composition xi ∈ [0, 1] and ∑ni=1xi = 1. Thus, the distance of the tangent plane to the molar Gibbs mixing energy surface, D(x),
 n Dm DðxÞ ¼ mðxÞ
m0
ðxi
zi Þ ð1Þ i ¼ 1 Dxi 0

ARTICLE

Figure 1. Tangent-plane distance surface in mole number (mol) coordinates for case 3 in the Appendix. The mixture studied is z = [0.1, 0.8, 0.1]. The mole number of n3 is constant throughout the tangent-plane distance surface (n3 = 0.01 mol). The horizontal plane indicates the zero plane of D(x).

Gibbs energy of mixing is given by mðxÞ ¼

ð2Þ

because excess Gibbs energy models are used. The first term on the right in eq 2 is convex, and the second term, i.e., the contribution of the mixture nonidealities, is nonconvex. The main difficulties in the minimization of D(x) arise, in general, from the nonconvex term. The Gibbs tangent-plane criterion is generally checked by minimizing D(x) with a mole fraction restriction. This minimization problem can be solved either directly or by solving an equivalent set of NLEs for the stationary points of D(x).15 An example of the D(x) surface can be viewed in Figure 1. The TPDF surface is divided into two parts, with the plane indicating the stability limit [D(x) = 0]. With certain test compositions, D(x) is negative, indicating that the mixture is unstable and will split into two or more phases depending on the subsequent flash calculations. However, in order to have feasible estimates for the individual phase compositions, at least the global minimum of D should be sought. In this paper, the global minimum of TPDF and all of the other roots of the problem are solved by creating a set of NLEs for the stationary points of D(x). Thus, the problem is described with the following equations:15 



 

 Dm Dm Dm Dm ¼ 0 i ¼ 1, :::, n
1


Dxi Dxn Dxi Dxn 0

ð3Þ

∑

is negative.15 Therefore, a phase with composition z is unstable in these conditions and will split into two or more phases. In eq 1, x is the mole fraction vector, n is the number of components, and m(x) is the reduced Gibbs energy of mixing. The subscript zero stands for evaluation at x = z. The reduced

n

∑ xi ln xi þ GE ðxÞ i¼1

1


n

∑ xi ¼ 0 i¼1

ð4Þ

The set of equations (n  1) has a trivial root, which corresponds to the composition of the inspected mixture (x = z). The equation set may also have multiple nontrivial feasible roots, i.e., mixture compositions that fulfill eqs 3 and 4 and the condition xi ∈ [0, 1]. Thus, the problem is a multiplicity problem. 7006

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research

ARTICLE

The number of roots depends on the case studied, i.e., the components, conditions, and the location of the mixture composition with respect to phase splitting, and the thermodynamic model used to describe the phase equilibrium. Finding only one root out of many or a negative TPDF value is not sufficient to ensure both accurate phase stability analysis and feasible estimates for subsequent phase-splitting calculations.38 The problem may have several roots, giving negative D(x) values, but only the smallest value corresponds to the global minimum. If the system is found to be unstable, the global minimum gives an initial estimate for the subsequent flash calculation.20 The global minimum and the other located stationary points can be used together to initiate the phase-splitting calculations efficiently (see, e.g., Sofyan et al.51). The consistency of the roots found, i.e., the stationary points, with respect to the TPDF surface topology, can be verified with the topological criterion presented by Wasylkiewicz et al.32 The topological criterion is defined for an n component mixture (d = n
1) as32 even odd Nmax þ ð
1Þd Nmin þ Nsad
Nsad ¼ ð
1Þd

ð5Þ

Thus, for a three-component mixture, the criterion reduces to the form32 Nmax þ Nmin
Nsad ¼ 1

ð6Þ

However, it should be noted that, in the topology criterion, it is not the absolute number of stationary points that is defined but rather the relative number of different stationary points. 3.2. Problem Solution. The phase stability analysis problem presented in eqs 3 and 4 is a set of NLEs that can be expressed as f ðxÞ ¼ 0

ð7Þ

In general, solving a system of NLEs is difficult for various reasons: • Some of the functions can be undefined for certain variable values, which creates a discontinuity (see, e.g., Shacham and Brauner52) or a boundary across the solution path. • The equation set may have solutions that are not physically feasible. • The functions may be highly nonlinear and/or badly scaled. • The feasible selection of the starting (or initial) point may be difficult. In addition, in process flowsheeting cases, the number of equations in the NLE set may be high and the Jacobian of the equation set may be sparse; i.e., the Jacobian matrix has a large portion of elements having a value of zero. This, in turn, may cause the solving of the NLE set to be inefficient unless the Jacobian matrix structure is taken into account in the numerical Jacobian matrix approximation. NLE sets are frequently solved iteratively with Newton
Raphson-based solution methods. The Newton and quasi-Newton methods are widely used local methods because they are simple and, theoretically at least, superlinearly convergent. These methods depend a great deal on the initial point and will not converge unless the initial point is reasonably close to the solution. In addition, these methods achieve only one solution from one initial point. The TPDF has multiple solutions, and selection of the initial point is not self-evident. The reasons for the difficult initial point selection also depend on other factors. In novel process schemes, the experimental data of the component behavior may be relatively limited or the user of the process simulation software may have

little experience with the software. Some heuristics have been proposed to tackle the difficulties; see, e.g., Michelsen.20 However, even with the use of numerous starting points along with the heuristics, it is not guaranteed that all of the minima of TPDF will be found.32 Thus, there is no theoretical guarantee of finding either the global minimum or any local minimum, proving the mixture instability for the TPDF and thus obtaining feasible initial estimates for subsequent phase-splitting calculations with local methods. The homotopy continuation methods enlarge generally the attraction domain with respect to Newton’s method and its reformulations. Thus, homotopy methods approach a solution or a multiple of them more robustly. This property has been demonstrated previously for solving of the phase stability analysis problem by Sun and Seider14 and Khaleghi and Jalali.16 3.2.1. Problem-Independent Homotopy Continuation Methods. Problem-independent homotopies are formed as a combination of the original problem function, f(x), and an auxiliary function, g(x). g(x) and f(x) are embedded with the help of a homotopy parameter, θ, to form a homotopy function h(x,θ). The first advantage of embedding is that the initial guess for h(x,θ) at θ = 0, i.e., the solution of g(x), may be obtained with a system of equations for which a solution is known or is easily obtained.8 The second advantage is that a continuous blend between the starting point and the solution of the original problem f(x) is created. The homotopy path is tracked with the help of the homotopy parameter from the known solution x0 at θ = 0 to the unknown solution x* at θ = 1, based on the homotopy equation hðx, θÞ ¼ θf ðxÞ þ ð1
θÞgðxÞ ¼ 0

ð8Þ

Depending on selection of the function g(x), the fixed-point (eq 9), Newton (eq 10), and scale-invariant affine (eq 11) homotopies are obtained:53 hðx, θÞ ¼ θf ðxÞ þ ð1
θÞðx
x 0 Þ

ð9Þ

hðx, θÞ ¼ θf ðxÞ þ ð1
θÞ½f ðxÞ
f ðx 0 Þ

ð10Þ

hðx, θÞ ¼ θf ðxÞ þ ð1
θÞf 0 ðx 0 Þ ðx
x 0 Þ

ð11Þ

The fixed-point and scale-invariant affine homotopies have the advantage that only one root satisfies h(x,θ) = 0 at θ = 0. However, the path can consist of branches that are only connected at the opposite infinities or branches in the complex domain. The Newton homotopy, on the other hand, has the advantage that a closed path for x on θ may exist within a finite domain. The main drawback is that the homotopy path may have multiple roots at θ = 0; i.e., more than one homotopy path branch may exist.53 In addition, the possibility of a starting point multiplicity may imply that the path will return to θ = 0 before reaching the first root. This is inefficient in relation to the computation efficiency. The homotopy continuation methods are particularly useful in solving problems in which sufficiently good starting guesses for the unknown variable values cannot be given. The homotopy continuation methods presented by, e.g., Allgower and Georg,54 Wayburn and Seader,8 Seader et al.,53 and Sun and Seider14 have in theory the property of global convergence, but they have their own shortcomings. One challenge is that the homotopy path may run partly or completely outside the feasible problem domain. Outside the 7007

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research

ARTICLE

domain, a problem may rise if the physicochemical properties of the system under study have to be evaluated and the procedure for doing so does not accept values outside the feasible domain. A typical example of this is calculation of the thermodynamic properties using thermodynamic relations built into process simulators. Thus, the computer subroutines evaluating the phase equilibrium constants and enthalpies may generate fatal errors.8 Another challenge encountered frequently in chemical engineering is boundary striking in the numerical solution of NLE sets. Boundary striking may occur if the solution path runs in the vicinity of the domain boundary without crossing it. Even though the correct solution path does not cross the domain boundary, the numerical predictor
corrector path-tracking algorithm may predict variable values outside the feasible domain. The boundary crossing occurs usually in homotopy continuation methods because of an unacceptably long predictor step. In some cases, the problems can be tackled with the use of absolute value functions8 and different mapping techniques (see, e.g., Seader et al.53). A promising proposal for overcoming the problems of the unboundedness and boundary striking of the homotopy path is the class of bounded homotopy methods presented by Paloschi.55 Bounded homotopy methods can be expressed in the form hb ðx, θÞ ¼ πðxÞ hðx, θÞ þ vðxÞ
vðx b Þ ¼ 0

ð12Þ

56

Paloschi extended the capabilities of bounded homotopies to handle also sparse problems. The difference compared to eq 12 is that the single penalty function value π(x) has been replaced with a penalty matrix P(x). The proposed methods yield a bounded homotopy path that preserves a Jacobian with the same sparsity pattern as the original problem. Thus, the possibility of exploiting sparse matrix techniques is maintained. This is important from an efficiency point of view. The main shortcoming of the bounded homotopies presented by Paloschi55,56 is that the tracking of the homotopy path inside the narrow bounding zone is numerically challenging. For example, a predicted step may fall outside the problem domain. In addition, the capability of the bounded homotopy methods to keep the homotopy path strictly bounded with respect to the problem variables is questionable in the form proposed by Paloschi. These weaknesses can be overcome with the modified bounded homotopies proposed by Malinen and Tanskanen.7 3.2.2. Modified Bounded Homotopies. The modified bounded homotopies presented by Malinen and Tanskanen7 are described in the form hmod ðxinf , θÞ ¼ Πðx inf Þ hðx inf , θÞ þ v ðxinf Þ
v ðxb, inf Þ b

¼0

x

x

ð13Þ The most significant difference between eqs 12 and 13 is that, in the modified bounded homotopies, the homotopy path is tracked with mapped variables. The idea is to track the homotopy path in an infinite space instead of the finite space bounded by the problem domain boundaries. Although the homotopy path is tracked in an infinite space, the variables are still mapped back into the finite space before substituting them into the original problem equations, f(x). Logarithmic mapping from a finite space into an infinite space is carried out for every variable xi of the problem on the basis of min the defined maximum (bmax i ) and minimum (bi ) values, as

presented in the following: ! 2ðxi
bmin inf i Þ , xi ¼ log max bi
bmin i

xinf i

1 þ bmin when xi < ðbmax i Þ ð14Þ 2 i

! 0:5ðbmax
bmin i i Þ ¼ log , bmax
xi i

1 þ bmin when xi g ðbmax i Þ 2 i ð15Þ

(bmax i )

(bmin i )

The maximum and minimum values can be realized as domain boundary values, either as real physical restrictions or as artificial values. By tracking of the homotopy path in the mapped variable space, path tracking becomes more flexible near the domain boundaries. Thus, the width of the bounding zone can be decreased in terms of the unmapped values. The implementation details on modified bounded homotopies can be found from the original papers of Malinen and Tanskanen7,57 and the work of Paloschi.55,56 3.2.3. Modified Bounded Homotopies in Phase Stability Analysis. The homotopy path may run outside the physically meaningful mole fraction domain with the traditional problemindependent homotopy continuation methods as the equation set formed by eqs 3 and 4 is solved. This is illustrated in Figure 2a, where the mole fraction of lauryl alcohol reaches almost the mole fraction value of 1.05 at θ ≈ 0.3. In this case, the mole fractions of the other components stayed inside the defined domain. In general, crossing the mole fraction boundary could cause the model solving to stop because of errors in the thermodynamic property calculation routines. However, in this case, the thermodynamic model equations are in a form that surpasses the usage of high nonphysical mole fraction values. It is worth noting that tracking of the path to the negative mole fraction domain will necessitate the usage of complex arithmetic. By utilization of the modified bounded Newton homotopy, the path can be restricted to run inside the feasible mole fraction domain [0, 1]. In addition, the width of the bounding zone can be set very narrow. In the case of Figure 2, a bounding zone width of 1  10
7 has been utilized. The advantage of a narrow bounding zone can be seen in cases in which one or multiple components are almost completely immiscible with other components and/or if some components are present in the mixture in trace amounts. Thus, one or more solutions of the phase stability analysis problem are located near the sides of the composition triangle (or polygon). Especially in these situations, the bounding zone should be set narrow. Figure 2 also illustrates how selection of the domain boundary min values, i.e., the maximum (bmax i ) and minimum (bi ) values, affects the course of the homotopy path in the mapped variable space. In the case of the unbounded homotopy path, the mole fraction domain has been specified as [0, 1.05], and in the case of the bounded homotopy path, as [0, 1]. In the mapped variable space, the change in the domain boundary values can be noticed as completely different paths. However, in the unmapped variable space, the bounded and unbounded Newton homotopy paths are equal except when the bounding routine is in force. It should be noted that in the case of the modified bounded fixed-point and affine homotopies the selection of domain and bmin has a crucial effect on the course boundary values bmax i i of the bounded homotopy path in the unmapped variable space, 7008

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research

ARTICLE

Figure 2. Bounded and unbounded Newton homotopy paths from the starting point to the first solution in an (a) unmapped and (b) mapped variable space for case 1 presented in the Appendix. The mixture composition studied is z = [0.2030, 0.7071, 0.0899]. The starting point in the unmapped variable space is x0 = [0.7455, 0.7932, 0.8746].

also when the bounding routine is not in force. When the bounding routine is not in force, eq 13 is reduced to the form inf inf inf hmod b ðx , θÞ ¼ hðx , θÞ ¼ θf ðxÞ þ ð1
θÞgðx Þ ¼0

ð16Þ

In eq 16, the difference between evaluation of the original and auxiliary equation sets is emphasized. As was already mentioned, the mapped variables are mapped back into the unmapped form before substitution into the original problem equations f(x). However, in fixed-point and affine homotopies, the function g(xinf) includes a term calculated only with the mapped variables (xinf
x0,inf). Thus, the course of the homotopy path is changed in the unmapped space through selection of the domain boundary values with fixed-point and affine homotopies, even inside the feasible domain.

4. STARTING POINT AND HOMOTOPY FUNCTION SELECTION From the convergence point of view, the favorable selection of the starting point is an important issue. Michelsen20 stated that one good choice for the starting point is a pure phase, which helps to detect the liquid immiscibility in highly immiscible cases. Sun and Seider14 used several different sets of starting points, ranging from the pure trial phase to an arithmetic mean of the overall composition. In the following sections, the convergence of modified bounded homotopies with the homotopy functions defined in eqs 9
11 is investigated. In addition, the attraction domain behaviors of the homotopies and the fsolve routine of Matlab, i.e., a local solver, are compared. The local solver uses the trust-region dogleg method to solve the NLE set. The emphasis in the study is to consider the effects of the starting point selection on the first root obtained and how the homotopy path is qualitatively changed in the presented cases. The reason for studying the performance of the methods with respect to the first root is that it would be highly beneficial to detect the phase-splitting phenomenon both effectively and with reasonable certainty. In particular, the effectiveness of the method is dependent on the rate of the path tracking. Thus, the contribution of the tracking algorithm to the overall solver

performance should be managed. A significant option in this respect is to track the homotopy path until the first root is obtained. If the tangent-plane distance at the first root is negative, the conclusion is that the mixture will split into two or more liquid phases. With this information, we can continue the phasesplitting examination by introducing representative equations for the new liquid phases in the system and solving the equation set subsequently. 4.1. Starting Point Selection. The implications of the change of the starting point are illustrated in Figures 3 and 4. These figures represent the first root obtained from each studied starting point with different modified bounded homotopy methods. The root has been approached by tracking the homotopy path starting in the positive homotopy parameter direction. The roots obtained, i.e., stationary points, have been classified as maxima, minima, and saddle points with evaluation of the Hessian matrix eigenvalues for TPDF. In addition, the topology criterion of Wasylkiewicz et al.,32 i.e., eq 6, has been used to verify the consistency of the results. On the basis of these results, some interesting observations can be made. 4.1.1. Vertex Point as a Starting Point. As can be noticed in Figures 3a
c and 4a
c, the vertices of the composition space triangle are not necessarily in the same area with respect to the first root reached; i.e., they are not in the same attraction domain. This implies that the solution of the problem can be performed effectively by initializing the calculation from the respective vertices. It is also noticeable that the global minimum corresponding root is found in both cases with the proposed initialization approach. However, the robustness of such a procedure has not been generally investigated in this study. In addition, the attraction domain of the global minimum is relatively large in both cases. Thus, the certainty of finding the root is high. It is noticeable that, with the local solver (fsolve routine of Matlab), the attraction domains cannot be clearly defined. Thus, starting from different vertex points sequentially will not necessarily result in different roots with local solving methods. Another point to be recognized relating to the issue includes the number of variables, i.e., the number of vertex points. The more vertex points that are found in the mixture studied, the more times the solver has to be initialized to have a reasonable certainty of finding feasible solutions, i.e., the local and global 7009

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research

ARTICLE

Figure 3. First root obtained from various starting points with modified bounded (a) affine, (b) fixed-point, and (c) Newton homotopies and (d) the fsolve routine of Matlab for case 1 presented in the Appendix. The components studied are ethylene glycol (1), lauryl alcohol (2), and nitromethane (3). The molar composition of the mixture, i.e., the trivial root (local minimum), studied is z = [0.270 78, 0.473 02, 0.256 20]. The composition of z is in the vicinity of one of the LLLE phase compositions (see the Appendix). The molar composition at each root is included in the legend.

minima of TPDF, to be used as initial estimates in subsequent phase-splitting calculations. In engineering terms, the path tracking may become the limiting factor in the efficiency (the calculation time) of the problem solving algorithm. 4.1.2. Starting Points near the Trivial Root. The knowledge of the trivial root location of the phase stability analysis problem is at the same time a bonus and a drawback. The trivial root is known before the start of the solution. Therefore, obtaining the trivial root as the first root along the homotopy path is inefficient with respect to the solution. Thus, the means to avoid approaching the trivial root should be found in the solving algorithms applied to phase stability analysis problems. The first alternative to prevent reaching the trivial root in a general phase stability analysis case is to choose the starting point based on a careful examination of the results in each case. The attraction domain of the trivial root is highly dependent on the case studied. In case 1 (see Figures 3a
c), the usage of a starting point near the trivial root, i.e., performing small perturbations to the starting point with regard to the trivial root, causes the solver

to reach the trivial root. However, in case 2 (see Figures 4a
c), the solver reaches the trivial root only from certain starting points in the vicinity of the trivial root. In both cases, the trivial root is a local minimum of the problem. The main difference between cases 1 and 2 is that the roots of the problem in case 2 (shown in Figure 7) are located in the vicinity of the trivial root. In contrast, in case 1, the roots are not located solely near the trivial root. Thus, it can be concluded that starting points in the vicinity of the trivial root do not necessarily result in the trivial root, but the behavior is case-specific. It should be noted that, theoretically speaking, starting the solution from the trivial root (or any other feasible root) forces the homotopy-based solver to converge to the same root. However, the attraction domain of a root may be very small. Thus, the solver may converge to another root even in a case in which very small perturbations are performed on the starting point. 4.2. Homotopy Function Selection. 4.2.1. Trivial Root as the First Root. A change of the homotopy function is another alternative 7010

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research

ARTICLE

Figure 4. First root obtained from various starting points with modified bounded (a) affine, (b) fixed-point, and c) Newton homotopies and (d) the fsolve routine of Matlab for case 2 presented in the Appendix. The components studied are n-propanol (1), n-butanol (2), and water (3). The molar composition of the mixture, i.e., the trivial root (local minimum), studied is z = [0.148, 0.052, 0.8]. The composition of z is in the vicinity of the plait point (see the Appendix). “No Root” indicates starting points from which no root whatsoever is obtained. The molar composition at each root is included in the legend.

to affecting the first root approached. The high probability of obtaining the trivial root as the first root is evident in Figures 3a
c and 4a
c. Selecting the fixed-point homotopy causes the area of the trivial root to be the largest in both cases. Therefore, usage of the fixed-point homotopy should be avoided in this respect. It is important to note for later analysis of the results in this paper that the trivial root is a local minimum, not a saddle point, in both cases 1 and 2. It is also worth noting that the shape of the trivial root attraction domain is (qualitatively) similar to that of the affine and Newton homotopies. In addition to changing the solver options (i.e., the starting point and homotopy function), the reformulation of the problem equations could also be a feasible alternative for preventing the approaching of the trivial root. This, in turn, could imply that the trivial root would be extracted from the feasible roots through mathematical modifications.

4.2.2. Unfeasible Roots as the First Root. Both Figures 3 and 4 show that, in addition, to the original problem roots, unfeasible roots (Figures 3a and 4a,c) may be approached as the first root. This behavior is shown in Figure 5 for the affine homotopy in case 2. The homotopy path travels into the bounding zone and stays inside the bounding zone until θ = 1. The appearance of unfeasible roots originates from the bounding effect of the modified bounded homotopies. If the path were not bounded, the path would travel outside the feasible domain. Thus, the unfeasible roots located in all of the cases of this study (see also Tables 1 and 2) indicate a situation in which the unbounded homotopy path, after running into the unfeasible problem domain, will not turn back into the feasible problem domain. In phase stability analysis problems, the feasible domain is typically defined as the mole fraction domain, xi ∈ [0, 1]. It should be noted that starting the solution toward the minus direction with respect to θ for the affine homotopy will not help 7011

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research

ARTICLE

in finding feasible roots, as illustrated by Figure 5. The reason is that only one homotopy branch exists at θ = 0 with the affine

homotopy. The reason for this is that g(x) is satisfied by only one root at θ = 0. However, it should be noted that this is true only inside the feasible domain, not including the bounding zone. In some cases, the path or another homotopy branch may cross θ = 0 inside the bounding zone. All of the cases in this paper also indicate that unfeasible roots are obtained only with affine and Newton homotopies but not with the fixed-point homotopy. As a whole, the probability of reaching an unfeasible root seems to be highest with the affine homotopy. Thus, the affine homotopy should be avoided in order to minimize the possibility of obtaining unfeasible roots. However, it should be noted that unfeasible roots do not fulfill the original problem equations, i.e., eqs 3 and 4. For example, in the case presented in Figure 5, it is observed that f(x*) ≈ [29.24, 30.11, 1
1.221  10
12]. Thus, the requirement of f(x*) = 0 is evidently conflicted. As a whole, unfeasible roots can easily be separated from feasible roots. 4.2.3. No Root. In Figure 4c, no root has been achieved, which, in turn, is certainly adverse for the solver performance. This behavior has been observed with the Newton homotopy applied for case 2, when the homotopy path approaches the homotopy parameter value of θ =
20. This indicates a starting point multiplicity in which the homotopy path continues until the minus infinity of θ is approached.

Figure 5. Homotopy path to the plus (strong curve) and minus directions (light curve) in case 2 with the affine homotopy. The studied mixture composition is z = [0.148, 0.052, 0.8]. The starting point at θ = 0 is x0 = [0.05, 0.1, 0.85] and the unfeasible root is located at θ = 1, x* = [1
1.017  10
12, 1
6.2  10
14, 1.42  10
13].

Table 1. Proportion of First Root Types with Various Modified Bounded Homotopies and the fsolve Routine of Matlab for Case 3 Presented in the Appendix with and without Starting Point Restrictiona proportion (%) Newton homotopy root

∑ix0i

=1

∑ix0i

6¼ 1

affine homotopy ∑ix0i

=1

∑ix0i

6¼ 1

isolab

13.7

4.5

trivial (local minimum)

13.4

12.5

14.5

12.2

saddle point

42.4

49.8

39.1

45.0

global minimum

13.5

9.2

15.1

11.2

unfeasible

0.1

0.6

31.3

31.6

no rootc,d

16.9

23.4

fixed-point homotopy ∑ix0i

=1

39.4 60.4

∑ix0i

6¼ 1

38.9 61.1

fsolve routine of Matlab ∑ix0i = 1

∑ix0i 6¼ 1

48.1

49.5

51.7

50.5

0.2

0.2e

a

1000 randomly (with the rand function of Matlab) selected starting points have been studied per option and function. The mixture under investigation was z = [0.299 89, 0.200 06, 0.500 05]. b No root found; a starting point isola is detected. c No root; the path crosses θ =
20. d With fsolve, x enters the negative mole fraction domain and the number of function evaluations reaches 2000. e Path tracking failure.

Table 2. Proportion of First Root Types with Various Modified Bounded Homotopies and fsolve Routine of Matlab for Case 4 Presented in the Appendix with and without Starting Point Restrictiona proportion (%) Newton homotopy

affine homotopy

fixed-point homotopy

fsolve routine of Matlab

root

∑ix0i = 1

∑ix0i 6¼ 1

∑ix0i = 1

∑ix0i 6¼ 1

∑ix0i = 1

∑ix0i 6¼ 1

∑ix0i = 1

∑ix0i 6¼ 1

trivial (saddle point) local minimum

60.4 21.3

74.2 13.9

60.1 21.8

67.1 13.9

48.2

70.5

75.4 17.0

81.0 12.4

global minimum

51.8

29.5

7.6

5.8

isolab

13.3

8.7

13.3

8.9

unfeasible

0.4

0.6

4.8

10.1

no rootc,d

4.6

2.6

0.8

a

1000 randomly (with the rand function of Matlab) selected starting points have been studied per option and function. The mixture under investigation was z = [0.2, 0.2, 0.2, 0.2, 0.2]. b No root found; a starting point isola is detected. c No root; the path crosses θ =
20. d With fsolve, x enters the negative mole fraction domain and the number of function evaluations reaches 8000. 7012

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research

Figure 6. Starting point isola in case 3 with the Newton homotopy. The starting point (*) is x0 = [0.666 134 708, 0.096 134 988, 0.237 730 303]. The starting point multiplicity (O) is located at θ = 0, x = [0.797 066 0, 0.029 672 9, 0.173 261 1]. The studied mixture composition is z = [0.299 89, 0.200 06, 0.500 05].

It is important to notice that no root is obtained also when a starting point isola is formed. This has not been the case in Figure 4c but has been found in case 3 (see Table 1). The starting point isola is defined here as a starting point multiplicity in which the homotopy path returns back to the starting point without crossing the θ = 1 plane, as shown in Figure 6. Both starting point multiplicity types are typical for the Newton homotopy and can be seen as problematic for usage of the Newton homotopy. However, the effect of the weaknesses of the Newton homotopy can be assuaged with certain provisions depicted below. 4.2.4. Tracking the Path in the Negative Direction. The problematic situation of obtaining either no root at all or an unfeasible root with the Newton homotopy can be alleviated in certain situations by tracking the path in the negative direction from the starting point with respect to the homotopy parameter. This issue has been studied in case 2 for the Newton homotopy, and the results are illustrated in Figure 7. From all of the starting points selected inside the shaded areas in Figure 7, the global minimum is approached as the first root by starting the path tracking in the negative θ direction. Thus, in the studied case 2, obtaining an unfeasible or no root is not detrimental to the performance of a Newton homotopy, but all of the starting points are connected to at least one root of the original problem. Clearly, the inherent property of the Newton homotopy of having starting point multiplicities may be used to our advantage. The selected starting point may be connected through the starting point multiplicities to feasible roots. However, in the case of the starting point isolas (section 4.2.3), it must be stressed that changing the direction of the homotopy path tracking will not help in approaching the feasible roots. However, in case 2, illustrated by Figures 4c and 7, the starting point isolas have not been detected. 4.2.5. Comparing the Attraction Domains of Various Homotopies. The attraction domain of the unfeasible root should be as small as possible to increase the efficiency of the calculation procedure. On the other hand, the attraction domain of the root corresponding to the global minimum should be as large as possible.

ARTICLE

Figure 7. Magnification of the “No Root” (dark shaded area) and “Unfeasible Root” (lightly shaded area) attraction domains of Figure 4c. Starting toward the minus direction from the points located inside the shaded areas results in approaching the global minimum. The composition studied is z = [0.148, 0.052, 0.8], marked with 0.

The probability of achieving the global minimum as the first root can be increased by changing the homotopy function. Figures 3 and 4 indicate that the fixed-point homotopy exhibits a different type of attraction domain than Newton and affine homotopies. The fixed-point function contracts the attraction domains of some of the roots, i.e., saddle points 1 and 2 in Figure 3b and the saddle point in Figure 4b. Thus, the probability of approaching the saddle points first from the specified starting point is diminished. However, the design of a solving procedure beforehand, such that the desired root is obtained first, seems to be challenging. The fixed-point homotopy does not locate saddle points and unfeasible roots as the first root but locates at least one minimum corresponding root of the original problem, as shown in Figures 3b and 4b. The observation that certain roots are not obtained as the first root with the fixed-point homotopy irrespective of the starting point can be considered either harmful or beneficial to the performance of the solving procedure. In all cases studied, the saddle points have positive or zero TPDF values, which would indicate that the mixture is stable at these conditions. The minima, in turn, have zero or negative TPDF values. Thus, in all of the cases, the mixture splits into two or more phases. Therefore, the observation of obtaining only the minima of TPDF is not harmful but beneficial from the perspective of phase stability analysis. Nonetheless, it should be recognized that while a general NLE set is solved, the lack of convergence to certain roots is not generally desirable. The remedy for obtaining the roots could be that, after approaching one root, path tracking is continued with the aim of also approaching the other roots. Second, if usage of the Newton homotopy results in attaining an unfeasible root or a starting point isola, changing to the fixed-point homotopy could aid in obtaining at least one feasible root for the problem. 4.3. Summation Restriction of the Starting Point. As the first guess, the starting point dependency of the homotopy methods could be used to our advantage with the provision that a heuristic rule be formed for the dependency. In previous sections, the dependency was studied with the attraction domains of each 7013

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research root. Now, let us consider the dependency from the perspective of the properties of the problem, i.e., the NLE set. The phase stability analysis problem includes the mole fraction summation equation as eq 4. Therefore, one promising selection could be to use either starting points with the mole fraction summation of 1 (naturally selected) or points having the mole fraction summation unequal to 1 (unnaturally selected). The mole fraction summation is utilized generally in deterministic process models, for example, the steady-state description of a distillation column tray. The possible gains could be high in the field of chemical engineering problem solving if the utilization of an unnaturally defined starting point would have a considerable and predictable effect on the attraction domains. The comparison between natural, i.e., at the starting point 6 1 and x0i ∈ [0, 1], ∑ni=1x0i = 1 and x0i ∈ [0, 1], and unnatural, ∑ni=1x0i ¼ starting point selection was investigated by tracking the homotopy path from the starting point in the positive direction with respect to the homotopy parameter θ. The first root obtained was stored. The results are summarized in Tables 1 and 2. Both Tables 1 and 2 indicate similar behavior. Choosing between the starting point types (natural or unnatural) does not significantly change the probability of approaching different roots. The probabilities are affected only slightly without a clear trend. However, a comparison between the fsolve routine of Matlab and homotopies shows that the probability of finding the global minimum is enhanced considerably with usage of the fixed-point homotopy. With Newton and affine homotopies, the increase is minor but, nonetheless, relevant. Another significant change is observed in Table 2 for the fixedpoint homotopy. The unnatural starting point selection increases the probability of reaching the local minimum with respect to the global minimum. The difference is detrimental to the performance of the solving method, but the generality of the qualitative behavior is currently undefined. Some other observations can also be made. The different homotopies converge toward different roots as the first root. The fixed-point homotopy approaches only the local or global minimum corresponding roots and lacks altogether the ability to reach other roots, which confirms the previous findings made based on Figures 3 and 4 for cases 1 and 2. On the other hand, the Newton homotopy has the drawback of a high proportion of starting point isola formation in case 3 indicated in Table 1. The affine homotopy reaches unfeasible roots at a high rate, which, in turn, reflects the probability of the homotopy path crossing the feasible domain boundaries with unbounded homotopies in both cases 3 and 4. Tables 1 and 2 are also representative of the similar behavior of the Newton and affine homotopy functions and the dissimilarity of the fixed-point function with respect to the others. The fixedpoint homotopy function has the benefit of converging only to feasible roots.

5. CONCLUSIONS Modified bounded homotopies assist in the determination of the phase stability by bounding the homotopy path with respect to the problem variables. Thus, fatal errors caused by unfeasible variable values in thermodynamic subroutines can be avoided. The mapping of the variables utilized in modified bounded homotopies helps in the tracking of the path near the mole fraction boundaries, which is both useful and typical for immiscibility cases. Even though modified bounded homotopies

ARTICLE

may induce unfeasible roots, fortunately unfeasible roots are easily separated from feasible roots. The choice of the homotopy function apparently has a significant impact on the first root approached from a starting point. The affine and Newton homotopy functions are relatively similar with regard to this behavior. Usage of the fixed-point homotopy function may lead to situations in which some of the roots are unattainable as the first root regardless of the starting point. However, in the cases studied, the unattainable roots were classified as saddle points, not as the global minimum or local minima. In addition, the TPDF values at all of the saddle points were positive or zero. Thus, the invaluable information regarding the phase stability of the considered mixture was not lost. Furthermore, it was observed that the attraction domain of the global minimum was considerably larger with the fixed-point homotopy than with the other homotopies and the fsolve (local) solver routine of Matlab. Hence, the findings imply that the fixedpoint homotopy should be favored as the first homotopy function option in the solving of the global minimum of TPDF. The main drawback of the affine homotopy function is that unfeasible roots are frequently obtained as the first root. The major problem of the Newton homotopy function is the formation of starting point isolas. Starting the solution in the negative direction with respect to the homotopy parameter increases the robustness of a Newton homotopy based solver. This is especially true in a situation where no root (excluding the starting point isolas) or an unfeasible root has been obtained with the Newton homotopy. The appearance of a starting point isola can be circumvented by choosing a new starting point or another homotopy function. The starting point has a clear effect on the performance. Generally, the vertices of the composition triangle or polygon are situated in different convergence areas. Thus, starting the calculations sequentially from pure (or near-pure) compositions results in convergence to different roots. The option is, therefore, a plausible but not completely certain way to increase the robustness in the phase stability determination. The random choice of the starting point, irrespective of the mole fraction summation limitation, has been found to be an unfeasible way to increase the robustness. According to the cases studied here, verification of the phase stability/instability of a liquid mixture into several liquid phases along with feasible initial estimates for the subsequent phasesplitting calculation is achieved with a high degree of certainty by (i) using a modified bounded fixed-point homotopy, in particular, if the trivial root is a saddle point, (ii) starting the solution sequentially from each vertex point, and (iii) tracking the homotopy paths until the first root is obtained. However, an important open question to be addressed by the user of the proposed heuristics is in terms of efficiency: Do we save time by tracking just one homotopy path branch and expecting to obtain all of the roots on the path, or do we track the homotopy path branches starting from all of the vertices consecutively until the first root is found and end there? This issue will hopefully be tackled in future work. It is worth noting that the proposed approach is not an efficient route to verifying the phase stability solely. The reason for this is that eq 4 is not necessarily fulfilled along the solution path. Thus, the value of TPDF cannot be evaluated along the path but only at the obtained roots. The efficiency of the suggested approach originates rather from the supply of feasible initial estimates for subsequent flash calculations. 7014

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research

ARTICLE

Table 3. UNIQUAC Binary and Pure Component Parameters for Case 1, the Ethylene Glycol (1)
Lauryl Alcohol (2)
Nitromethane (3) Mixture, τij = exp(
Δuij/RT)a

Table 8. LLE Compositions and Phase Fractions for n-Propanol (1)
n-Butanol (2)
Water (3) Composition z = [0.148, 0.052, 0.8]

Δuij/R i

j=1

1

a

j=2

0

ri

qi = qi0

2.4088

2.248

j=3

247.2

54.701

2

69.69

0

305.52

8.8495

7.372

3

467.88

133.19

0

2.0086

1.868

The parameters are taken from McDonald and Floudas.5 P = 1 atm.

phase fraction

L1

0.171 50

0.116 73

0.037 07

0.846 20

L2

0.828 50

0.154 47

0.055 09

0.790 44

Table 9. Binary Parameters of the NRTL Model for Case 3, the Toluene (1)
Water (2)
Aniline (3) Mixturea components, i
j

Table 4. Feasible Roots and Corresponding Values of D for Ethylene Glycol (1)
Lauryl Alcohol (2)
Nitromethane (3) Composition z = [0.270 78, 0.473 020, 0.256 20]a root, [x1, x2, x3] 0.270 780

0.473 020

0.256 20

feed/trivial,

a

0.156 741

0.496 339

τij

j

τij

Rij = Rji

C7H8
H2O

1

2

4.930 35

7.770 63

0.2485

C7H8
C6H7N

1

3

1.598 06

0.035 09

0.3000

H2O
C6H7N

2

3

4.184 62

1.279 32

0.3412

The temperature is 298 K, and the pressure is 1 atm. The parameters were taken from McDonald and Floudas.5

0

local minimum 0.346 920

i

a

D

type

L, [x1, x2, x3]

liquid phase

saddle point 1

0.034 544
6

0.619 885 0.023 341

0.005 619 0.001 726

0.374 497 0.974 933

local minimum global minimum


4.1164  10
0.058 761

0.369 416

0.015 329

0.615 256

saddle point 2

7.9576  10
3

Table 10. Feasible Roots and Corresponding Values of D for Toluene (1)
Water (2)
Aniline (3) Composition z = [0.299 89, 0.200 06, 0.500 05] root, [x1, x2, x3]

The temperature is 295 K. P = 1 atm.

0.299 89

0.200 06

D

type 0.500 05

trivial, local

0

minimum

Table 5. LLL Equilibrium Compositions for Ethylene Glycol (1)
lauryl Alcohol (2)
Nitromethane (3) Composition z = [0.270 78, 0.473 020, 0.256 20]a

a

L, [x1, x2, x3]

liquid phase

phase fraction

L1

0.961 43

0.278 99

0.491 91

0.229 10

L2

0.036 35

0.027 76

0.002 06

0.970 18

L3

0.002 22

0.692 80

0.003 99

0.303 21

components, i
j

i

j

τij

τij

Rij = Rji

C3H8O
C4H10O

1

2


0.612 59

0.716 40

0.30

C3H8O
H2O

1

3


0.071 49

2.742 5

0.30

C4H10O
H2O

2

3

0.900 47

3.513 07

0.48

6.693 7  10

Liquid phase Phase fraction

a

0.052

0.8

trivial, local minimum

0.114 336 0.035 993 0.849 671 global minimum 0.143 614 0.049 882 0.806 503 saddle point

D 0 4.5711  10
8

In addition, the solving of a phase stability analysis problem using modified bounded homotopy methods is worth studying in terms of approaching all of the roots robustly. This involves an indepth assessment of the effect of the auxiliary function, which

0.13515

0.34674

0.07584

0.57742

9.13  10
3 0.9949475 4.9612  10
3

The temperature is 295 K, and P = 1 atm.

i\j


9.8510  10
6

0.86485

L2

τij

a

0.148

L1

L, [x1,x2,x3]

Table 12. Binary and Pure Component Parameters of the UNIQUAC Model for Case 4, the Acetic Acid (1)
Benzene (2)
Furfural (3)
Cyclohexane (4)
Water (5) Mixturea

Table 7. Feasible Roots and Corresponding Values of D for n-Propanol (1)
n-Butanol (2)
Water (3) Composition z = [0.148, 0.052, 0.8] type

0.996 865 0.003 068 global minimum
0.294 54

Table 11. LLE Compositions and Phase Fractions for Toluene (1)
Water (2)
Aniline (3) Composition z = [0.299 89, 0.200 06, 0.500 05]a

The parameters are taken from McDonald and Floudas.5

root, [x1, x2, x3]

2.969 8  10
7

0.209 041 0.498 409 saddle point
5

The temperature is 295 K. P = 1 atm.

Table 6. Binary NRTL Parameters for Case 2, the n-Propanol (1)
n-Butanol (2)
Water (3) Mixturea

a

0.292 55

ri

qi = qi0

1.54662

2.2024

2.072

0.09441

3.1878

2.400

0.96249

0.60488

3.1680

2.484

0.26222

1.0

0.08839

4.0464

3.240

0.69066

0.19491

1.0

0.9200

1.400

1

2

3

4

1

1.0

1.26362

3.36860

0.85128

2

0.99972

1.0

1.02041

0.89333

3

0.31633

0.79027

1.0

4

0.49739

1.09619

5

2.44225

0.13507

5

The parameters were taken from Tessier et al.15

controls the path in the bounding zones. A means to diminish the attraction domain of the trivial root should be sought in future algorithms. An alternative is to utilize preknowledge of the trivial root location. Furthermore, usage of the concept of the homotopy parameter bounding, as presented by Malinen and Tanskanen,58 should be implemented and evaluated in phase stability analysis. 7015

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research

ARTICLE

Table 13. Feasible Roots and Corresponding Values of D for Acetic Acid (1)
Benzene (2)
Furfural (3)
Cyclohexane (4)
Water (5) Composition z = [0.2, 0.2, 0.2, 0.2, 0.2] root, [x1, x2, x3, x4, x5]

D

type

0.2

0.2

0.2

0.2

0.2

trivial, saddle point

0.227 485 0.030 397

0.003 345 0.302 415

0.047 503 0.034 395

0.001 584 0.622 212

0.720 083 0.010 582

global minimum local minimum

Table 14. LLE Compositions and Phase Fractions for Acetic Acid (1)
Benzene (2)
Furfural (3)
Cyclohexane (4)
Water (5) Composition z = [0.2, 0.2, 0.2, 0.2, 0.2] liquid phase

phase fraction

L1

0.466 59 0.094 857 0.352 995 0.133 926 0.394 035 0.024 187

L2

0.533 41 0.291 972 0.066 171 0.257 797 0.030 272 0.353 788

L, [x1, x2, x3, x4, x5]

’ APPENDIX The appendix presents the interaction parameters utilized for the thermodynamic models, the stationary points and respective TPDF function values, and the phase equilibrium compositions for the cases studied. Case 1: Ethylene Glycol (1), Lauryl Alcohol (2), and Nitromethane (3). This case was investigated by McDonald and

Floudas,5 Tessier et al.,15 and Chakravarty et al.,59 who supplied the parameters for the UNIQUAC activity coefficient model. The binary and pure component parameters are presented in Table 3. Tessier et al.15 presented the stationary points (i.e., roots) of the problem for a number of tested mixture compositions. Stateva et al.60 studied the LLLE behavior in the case but applied the original UNIFAC model. The feasible roots and corresponding values of D for composition z = [0.270 78, 0.473 020, 0.256 20] are presented in Table 4. The LLLE compositions in this case are presented in Table 5 taken from McDonald and Floudas.5 In this case (Table 4), d = 2, Nmax = 0, Nmin = 3, and Nsad = 2. Thus, according to the topology criterion presented by Wasylkiewicz et al.,32 the number of stationary points is consistent.

Case 2: n-Propanol (1), n-Butanol (2), and Water (3). The case was investigated, among others, by McDonald and Floudas,5 Walraven and Van Rompay,61 and Tessier et al.15 The binary parameters given by McDonald and Floudas5 are presented in Table 6. The calculated LLE diagram for the system is presented in Walraven and Van Rompay.61 The feasible roots and corresponding values of D for composition z = [0.148, 0.052, 0.8] are presented in Table 7. The studied mixture composition lies in the vicinity of the plaint point of the LLE area (see Walraven and Van Rompay61) and Table 8. In this case (Table 7), d = 2, Nmax = 0, Nmin = 2, and Nsad = 1. Thus, according to the topology criterion presented by Wasylkiewicz et al.,32 the number of stationary points is consistent. Case 3: Toluene (1), Water (2), and Aniline (3). The phase behavior is represented with the NRTL activity coefficient model. The parameters were acquired from McDonald and Floudas5 and are listed in Table 9. The LLE compositions in this case are presented in Table 11. The feasible roots and corresponding values of D for composition

0
0.177 65
0.113 54

z = [0.299 89, 0.200 06, 0.500 05] are presented in Table 10. In this case (Table 10), d = 2, Nmax = 0, Nmin = 2, and Nsad = 1. Thus, according to the topology criterion presented by Wasylkiewicz et al.,32 the number of stationary points is consistent. Case 4: Acetic Acid (1), Benzene (2), Furfural (3), Cyclohexane (4), and Water (5). The UNIQUAC activity coefficient

model is used with parameters obtained from Tessier et al.15 The parameters are listed in Table 12. The feasible roots and corresponding values of D for composition z = [0.2, 0.2, 0.2, 0.2, 0.2] are presented in Table 13. The LLE compositions in this case are presented in Table 14. In this case (Table 13), d = 4, Nmax = 0, Nmin = 2, and Nsad = 1. Thus, according to the topology criterion presented by Wasylkiewicz et al.,32 the number of stationary points is consistent.

’ AUTHOR INFORMATION Corresponding Author

*Tel.: þ358 8 553 2344. Fax: þ358 8 553 2304. E-mail: jani. kangas@oulu.fi.

’ NOMENCLATURE b=domain boundary D=tangent-plane distance d=dimension of the tangent-plane distance surface f=equation set f0 =Jacobian matrix of f G=reduced molar Gibbs free energy g=auxiliary function in the homotopy equation h=homotopy function m=reduced Gibbs free energy of mixing N=number of stationary points n=number of components, dimension of the vector, or number of moles, mol P=pressure, atm P=penalty matrix q=relative molecular surface area, qi = qi0 , used in the UNIQUAC model R=ideal gas constant, 8.314 J/(mol 3 K) r=relative molecular size used in the UNIQUAC model T=temperature, K u=interaction energy of molecules used in the UNIQUAC model, J/mol v=auxiliary function in the bounded homotopy equation x=variable vector z=composition of the studied mixture and trivial root Greek Letters

r=constant or binary interaction parameter in the NRTL model θ=homotopy parameter π=penalty function Π=penalty matrix 7016

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research τ=binary parameter evaluated from experimental data in the NRTL or UNIQUAC model Subscripts

b=bounded i=ith ij=interaction between the ith and jth molecules j=jth max=maximum min=minimum x=with respect to problem variables sad=saddle point 0=evaluation of the matrix at the tested mixture composition or starting point calculation Superscripts

b=bounded E=excess even=term in eq 5 included if the number of negative eigenvalues is even i=ith inf=infinity L=liquid max=maximum min=minimum mod=modified n=dimension or number of components odd=term in eq 5 included if the number of negative eigenvalues is odd 0=starting point *=solution point and root 0 =modified in the UNIQUAC model

’ REFERENCES (1) Saber, N.; Shaw, J. M. Rapid and robust phase behaviour stability analysis using global optimization. Fluid Phase Equilib. 2008, 264, 137. (2) Bekiaris, N.; G€uttinger, T. E.; Morari, M. Multiple Steady States in Distillation: Effect of VL(L)E Inaccuracies. AIChE J. 2000, 46, 955. (3) Wasylkiewicz, S. K.; Doherty, M. F.; Malone, M. F. Computing All Homogeneous and Heterogeneous Azeotropes in Multicomponent Mixtures. Ind. Eng. Chem. Res. 1999, 38, 4901. (4) Mitsos, A.; Bollas, G. M.; Barton, P. I. Bilevel optimization formulation for parameter estimation in liquid
liquid phase equilibrium problems. Chem. Eng. Sci. 2009, 64, 548. (5) McDonald, C. M.; Floudas, C. A. Global Optimization for the Phase Stability Problem. AIChE J. 1995, 41, 1798. (6) Xu, G.; Haynes, W. D.; Stadtherr, M. A. Reliable phase stability for asymmetric models. Fluid Phase Equilib. 2005, 235, 152. (7) Malinen, I.; Tanskanen, J. Modified bounded homotopies to enable a narrow bounding zone. Chem. Eng. Sci. 2008, 63, 3419. (8) Wayburn, T. L.; Seader, J. D. Homotopy continuation methods for computer-aided process design. Comput. Chem. Eng. 1987, 11, 7. (9) Wakeham, W. A.; Stateva, R. P. Numerical Solution of the Isothermal, Isobaric Phase Equilibrium Problem. Rev. Chem. Eng. 2004, 20, 1. (10) Rangaiah, G. P. Evaluation of genetic algorithms and simulated annealing for phase equilibrium and stability problems. Fluid Phase Equilib. 2001, 187
188, 83. (11) Baker, L. E.; Pierce, A. C.; Luks, K. D. Gibbs Energy Analysis of Phase Equilibria. Soc. Pet. Eng. J. 1982, 22, 731. (12) Tolsma, J. E.; Barton, P. I. Computation of heteroazeotropes. Part I: Theory. Chem. Eng. Sci. 2000, 55, 3817. (13) Lucia, A.; DiMaggio, P. A.; Bellows, M. L.; Octavio, L. M. The phase behavior of n-alkane systems. Comput. Chem. Eng. 2005, 29, 2363.

ARTICLE

(14) Sun, A. C.; Seider, W. D. Homotopy-continuation method for stability analysis in the global minimization of the Gibbs free energy. Fluid Phase Equilib. 1995, 103, 213. (15) Tessier, S. R.; Brennecke, J. F.; Stadtherr, M. A. Reliable phase stability analysis for excess Gibbs energy models. Chem. Eng. Sci. 2000, 55, 1785. (16) Khaleghi, S.; Jalali, F. Multiple Solutions in Stability Analysis Using Homotopy Continuation in Complex Space. Chem. Eng. Commun. 2007, 194, 1241. (17) Zhu, Y.; Xu, Z. Calculation of Liquid
Liquid Equilibrium Based on the Global Stability Analysis for Ternary Mixtures by Using a Novel Branch and Bound Algorithm: Application to UNIQUAC Equation. Ind. Eng. Chem. Res. 1999, 38, 3549. (18) Zhu, Y.; 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 Equilib. 1999, 154, 55. (19) Bonilla-Petriciolet, A.; Vazquez-Roman, R.; Iglesias-Silva, G. A.; Hall, K. R. Performance of Stochastic Global Optimization Methods in the Calculation of Phase Stability Analyses for Nonreactive and Reactive Mixtures. Ind. Eng. Chem. Res. 2006, 46, 4764. (20) Michelsen, M. L. The isothermal flash problem. Part I: Stability. Fluid Phase Equilib. 1982, 9, 1. (21) Nagarajan, N. R.; Cullick, A. S.; Griewank, A. New Strategy for Phase Equilibrium and Critical Point Calculations by Thermodynamic Energy Analysis. Part I. Stability Analysis and Flash. Fluid Phase Equilib. 1991, 62, 191. (22) Mitsos, A.; Barton, P. I. A Dual Extremum Principle in Thermodynamics. AIChE J. 2007, 53, 2131. (23) Eubank, P. T.; Elhassan, A. E.; Barrufet, M. A.; Whiting, W. B. Area Method for Prediction of Fluid-Phase Equilibria. Ind. Eng. Chem. Res. 1992, 31, 942. (24) Balogh, J.; Craven, R. J. B.; Stateva, R. P. The Area Method for Phase Stability Revisited: Further Developments. Formulation in Terms of the Convex Envelope of Thermodynamic Surfaces. Ind. Eng. Chem. Res. 2007, 46, 1611. (25) Balogh, J.; Csendes, T.; Stateva, R. P. Application of a stochastic method to the solution of the phase stability problem: cubic equations of state. Fluid Phase Equilib. 2003, 212, 257. (26) Boender, C. G. E.; Rinnooy Kan, A. H. G.; Timmer, G. T.; Stougie, L. A Stochastic Method for Global Optimization. Math. Program 1982, 22, 125. (27) Zhu, Y.; Wen, H.; Xu, Z. Global stability analysis and phase equilibrium calculations at high pressures using the enhanced simulated annealing algorithm. Chem. Eng. Sci. 2000, 55, 3451. (28) Srinivas, M.; Rangaiah, G. P. A study of differential evolution and tabu search for benchmark, phase equilibrium and phase stability problems. Comput. Chem. Eng. 2007, 31, 760. (29) Hua, J. Z.; Brennecke, J. F.; Stadtherr, M. A. Reliable Prediction of Phase Stability Using an Interval Newton Method. Fluid Phase Equilib. 1996, 116, 52. (30) Hua, J. Z.; Brennecke, J. F.; Stadtherr, M. A. Enhanced Interval Analysis for Phase Stability: Cubic Equation of State Models. Ind. Eng. Chem. Res. 1998, 37, 1519. (31) Gecegormez, H.; Demirel, Y. Phase stability analysis using interval Newton method with NRTL model. Fluid Phase Equilib. 2005, 237, 48. (32) Wasylkiewicz, S. K.; Sridhar, L. N.; Doherty, M. F.; Malone, M. F. Global Stability Analysis and Calculation of Liquid
Liquid Equilibrium in Multicomponent Mixtures. Ind. Eng. Chem. Res. 1996, 35, 1395. (33) McDonald, C. M.; Floudas, C. A. Global Optimization for the Phase and Chemical Equilibrium Problem: Application to the NRTL Equation. Comput. Chem. Eng. 1995, 19, 1111. (34) Floudas, C. A.; Visveswaran, V. A Global Optimization Algorithm (GOP) for Certain Classes of Nonconvex NLPs—I. Theory. Comput. Chem. Eng. 1990, 14, 1397. (35) McDonald, C. M.; Floudas, C. A. Global optimization and analysis for the Gibbs free energy function using the UNIFAC, Wilson and ASOG equations. Ind. Eng. Chem. Res. 1995, 34, 1674. 7017

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018

Industrial & Engineering Chemistry Research (36) McDonald, C. M.; Floudas, C. A. GLOPEQ: A New Computational Tool for the Phase and Chemical Equilibrium Problem. Comput. Chem. Eng. 1997, 21, 1. (37) Nichita, D. V.; Gomez, S.; Luna, E. Multiphase equilibria calculation by direct minimization of Gibbs free energy with a global optimization method. Comput. Chem. Eng. 2002, 26, 1703. (38) Nichita, D. V.; Gomez, S. Efficient location of multiple global minima for the phase stability problem. Chem. Eng. J. 2009, 152, 251. (39) Fidkowski, Z. T.; Malone, M. F.; Doherty, M. F. Computing Azeotropes in Multicomponent Mixtures. Comput. Chem. Eng. 1993, 17, 1141. (40) Eckert, E.; Kubicek, M. Computing Heterogeneous Azeotropes in Multicomponent Mixtures. Comput. Chem. Eng. 1997, 21, 347. (41) Tolsma, J. E.; Barton, P. I. Computation of heteroazeotropes. Part II: Efficient calculation of changes in phase equilibrium structure. Chem. Eng. Sci. 2000, 55, 3835. (42) Wasylkiewicz, S. K.; Kobylka, L. C.; Castillo, F. J. L. Pressure Sensitivity Analysis of Azeotropes. Ind. Eng. Chem. Res. 2003, 42, 207. (43) Wang, M. C.; Wong, D. S. H.; Chen, H.; Yan, W.; Guo, T.-M. Homotopy continuation method for calculating critical loci of binary mixtures. Chem. Eng. Sci. 1999, 54, 3873. (44) Aslam, N.; Sunol, A. K. Sensitivity of azeotropic states to activity coefficient model parameters and system variables. Fluid Phase Equilib. 2006, 240, 1. (45) Aslam, N.; Sunol, A. K. Reliable computation of binary homogeneous azeotropes of multicomponent mixtures at higher pressures through equation of states. Chem. Eng. Sci. 2004, 59, 599. (46) Bausa, J.; Marquardt, W. Quick and reliable phase stability test in VLLE flash calculations by homotopy continuation. Comput. Chem. Eng. 2000, 24, 2447. (47) Steyer, F.; Flockerzi, D.; Sundmacher, K. Equilibrium and ratebased approaches to liquid
liquid phase splitting calculations. Comput. Chem. Eng. 2005, 30, 277. (48) Jalali-Farahani, F.; Seader, J. D. Use of homotopy-continuation method in stability analysis of multiphase, reacting systems. Comput. Chem. Eng. 2000, 24, 1997. (49) Jalali, F.; Seader, J. D. Homotopy continuation method in multi-phase multi-reaction equilibrium systems. Comput. Chem. Eng. 1999, 23, 1319. (50) Jalali, F.; Seader, J. D.; Khaleghi, S. Global solution approaches in equilibrium and stability analysis using homotopy continuation in the complex domain. Comput. Chem. Eng. 2008, 32, 2333. (51) Sofyan, Y.; Ghajar, A. J.; Gasem, K. A. M. Multiphase equilibrium calculations using Gibbs minimization techniques. Ind. Eng. Chem. Res. 2003, 42, 3786. (52) Shacham, M.; Brauner, N. Numerical solution of non-linear algebraic equations with discontinuities. Comput. Chem. Eng. 2002, 26, 1449. (53) Seader, J. D.; Kuno, M.; Lin, W.-J.; Johnson, S. A.; Unsworth, K.; Wiskin, J. W. Mapped continuation methods for computing all solutions to general systems of nonlinear equations. Comput. Chem. Eng. 1990, 14, 71. (54) Allgower, E.; Georg, K. Simplicial and Continuation Methods for Approximating Fixed Points and Solutions to Systems of Equations. SIAM Rev. 1980, 22, 28. (55) Paloschi, J. R. Bounded Homotopies to Solve Systems of Algebraic Nonlinear Equations. Comput. Chem. Eng. 1995, 19, 1243. (56) Paloschi, J. R. Bounded Homotopies to Solve Systems of Sparse Algebraic Nonlinear Equations. Comput. Chem. Eng. 1997, 21, 531. (57) Malinen, I.; Tanskanen, J. A Rigorous Minimum Energy Calculation Method for a Fully Thermally Coupled Distillation System. Chem. Eng. Res. Des. 2007, 85A, 502. (58) Malinen, I.; Tanskanen, J. Homotopy parameter bounding in increasing the robustness of homotopy continuation methods in multiplicity studies. Comput. Chem. Eng. 2010, 34, 1761.

ARTICLE

(59) Chakravarty, T.; White, C. W., III; Seider, W. D. Computation of Phase Equilibrium: Optimization with Thermodynamic Inconsistency. AIChE J. 1985, 31, 316. (60) Stateva, R. P.; Cholakov, G. S.; Galushko, A. A.; Wakeham, W. A. A powerful algorithm for liquid
liquid
liquid equilibria predictions and calculations. Chem. Eng. Sci. 2000, 55, 2121. (61) Walraven, F. F. Y.; Van Rompay, P. V. An improved phasesplitting algorithm. Comput. Chem. Eng. 1988, 12, 777.

7018

dx.doi.org/10.1021/ie101907h |Ind. Eng. Chem. Res. 2011, 50, 7003–7018