Available online at www.sciencedirect.com
Computers and Chemical Engineering 32 (2008) 2333–2345
Global solution approaches in equilibrium and stability analysis using homotopy continuation in the complex domain Farhang Jalali a,∗ , J.D. Seader b , Saeed Khaleghi a a
Chemical Engineering Department, University of Tehran, 16 Azar Street, Tehran, Iran b Chemical Engineering Department, University of Utah, Salt Lake City, UT, USA
Received 27 April 2007; received in revised form 1 December 2007; accepted 3 December 2007 Available online 15 December 2007
Abstract Recently, significant attention has been shown to physical and chemical equilibrium and stability analysis in the real and complex domains. In this work, a new procedure involving the continuation method in the complex domain using bifurcation theory is propounded. Based on this method, homotopy branches in real and complex space are connected to each other through bifurcation branches. Thus, by just one initial guess, multiple solution branches are found. When calculations are only made in the real domain, multiple solutions are not always found from an arbitrary initial guess. Examples are presented to show the application of the method to nonlinear sets of equations in phase equilibrium, chemical and phase equilibrium, and stability analysis. These types of problems are believed to contain significant nonlinearities in process simulations. The results can be applied to flowsheet calculations. © 2008 Published by Elsevier Ltd. Keywords: Global solution; Equilibrium; Stability analysis; Homotopy continuation; Complex domain
1. Introduction Over the past three decades, many computer-aided methods for process design and control in the chemical industry have evolved. Today, chemical processes are simulated with an increasing number of rigorous, complex models. Among the fundamental tasks in computer-aided design of chemical processes is that of solving the large-scale, nonlinear-algebraic-equation systems involved in steady-state flowsheet modeling, optimization, and stability analysis. Two reasons may be given for the interest in this task. The first comes from a broadly recognized need to automate numerical modeling. The second reason is that a reliable method for solving the nonlinear-equation-system models still needs to be improved. Several difficulties in converging iterative calculations have been experienced using available software and that software typically does not search for multiple solutions, among which is the desired solution. Programming problems form a major subset of the field of mathematical programming. In particular, problems related to
∗
Corresponding author. Tel.: +9821 88632975; fax: +9821 88632976. E-mail address:
[email protected] (F. Jalali).
0098-1354/$ – see front matter © 2008 Published by Elsevier Ltd. doi:10.1016/j.compchemeng.2007.12.001
chemical process design and control can often be formulated as nonlinear optimization problems. These problems usually have nonconvexities that imply difficulties in determining a global solution. Chemical engineering examples of such problems are mostly in reactor network synthesis, phase equilibrium, simultaneous phase and chemical reaction equilibrium, stability analysis, heat exchanger network design, batch processing, and optimal design of distillation operations. Much of the design literature is related to solving the material and energy balances, phase equilibrium, transport equations, stability of equilibrium, and chemical kinetic expressions for these types of problems. The objective of the work reported here was to study the behavior of the continuation method using bifurcation theory in the complex domain for multivariable phase equilibrium, chemical and phase equilibrium, and stability analysis. It is shown that the non-convergent behavior, sometimes observed in real space, can be transferred through the complex domain to the real domain to find more solutions, among which may be the best solution. Here, global solutions are approached for the equilibrium and stability problem for the case when nonideal phases can be modeled using different activity-coefficient equations. This is an important result because if a global approach shows that
2334
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345
1.1. Mathematics of the computation Nomenclature A
y0 yi z zi
coefficient in the van Laar activity coefficient model coefficient of element k in component i activity of component i in phase j coefficient in the van Laar activity coefficient model total atoms of element k in the system number of components in the system number of elements function, objective function function Gibbs function (free energy) partial molar Gibbs function or chemical potential standard Gibbs function of component i in phase j function of x homotopy function equality constraint function complex number total number of moles in the system number of moles of component i in phase j number of moles of component i in the feed number of phases in the system gas constant real space inequality function temperature homotopy parameter mole fraction, general variable initial value of x mole number defined as exp (−K)yi imaginary part of complex variable, mole fraction vector in trial phase initial value of y mole fraction of component i complex variable, mole fraction vector mole fraction of i in the original system
Greek G γ ij ∂ ε λ μ
change in Gibbs free energy activity coefficient of component i in phase j partial deferential sign amount of the second phase vector of Kuhn–Tucker multiplier vector of chemical potential
Aik aij B bk C E f F G ¯ ij G Goij g(x) H h i NT nij nif P R Rn S T t x x0 Yi y
a solution satisfying the condition of equal chemical potentials is stable with respect to all possible perturbations, then global equilibrium and stability can be definitively asserted. Examples are provided that demonstrate the success in finding global solutions of the equilibrium and stability problems, regardless of the assumed starting point.
For steady-state analysis, the resulting nonlinear equation system, f (x- ) = 0, is solved in process simulators, mostly using Newton-based methods. Here, f is a vector function of vector x. Nonlinear equations, however, have many roots in general, and it is highly desirable in many chemical engineering problems to find all real solutions. Methods for finding additional roots by continuation are discussed by Allgower and Georg (1980). By far, the most widely used and successful of these approaches are those for finding all roots, real and complex of systems of polynomial equations by a global homotopy method (e.g. Watson, Billups, & Morgan, 1987). Here, the emphasis is on general systems of nonlinear algebraic equations. 1.2. Continuation method The continuation method is an important advance in the development of methods for solving general systems of nonlinear algebraic and transcendental equations. Potentially, it is a global method that has the capability of finding all the roots to a nonlinear system, from which physically meaningful and optimal solutions can be selected. This aspect has been discussed by Chow, Mallet-Paret, and Yorke (1978), Garcia and Zangwiil (1979), Watson (1986), Morgan (1987), Li (1987), Watson et al. (1987), Wayburn and Seader (1987), Seader et al. (1990), Sun and Seider (1995), Watson, Sosonika, Melville, Morgan, and Walker (1997), Jalali and Seader (1999), Jalali and Seader (2000), Bausa and Marquardt (2000), Wise, Sommese, and Watson (2000), Watson (2000), Gritton, Seader and Lin (2001), Li (2003), and Wu et al. (2005). 1.3. Definition of homotopy function A system of nonlinear algebraic equations can be written as: f (x- ) = 0 (1) where f:Rn → Rn is a continuously differentiable function, and x ∈ Rn is the variable vector. A homotopy function can provide a smooth transition between an approximation to the solution (often linear or nearly linear) and the true solution by gradually introducing the nonlinearities. In many formulations, the homotopy function is defined as: H(x, t) = tf (x) + (1 − t)g(x)
(2)
where f(x) is the system of equations to be solved and g(x) is a simple system of equations for which a solution is known or easily found. Here, t is a scalar homotopy parameter, which is gradually varied from 0 to 1 as the path is tracked from the staring point to a solution. Here, that range is extended to obtain all roots to the system of equations. The importance of the continuation technique comes from the disadvantage of Newton’s method, which can generate a singularity, or near-singularity of the Jacobian matrix and, thus, prevent convergence. However, recent work by Lucia and Feng
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345
2335
(2002, 2003) has shown that global terrain methods can overcome these limitations of Newton’s method and can find both physically meaningful solutions and singular points. In highly nonlinear systems of equations, a large or even global convergence domain is often desirable. This is especially true when the algorithm is part of a large simulation system where a good initial estimate may be difficult to provide and a failure to converge may result. Wayburn and Seader (1987) discuss the conditions under which homotopy continuation is guaranteed to converge, but cannot guarantee to find all of the solutions. In order to expand the domain of the continuation, calculation in the complex domain of homotopy continuation is receiving considerable attention as a means to find all solutions. That technique is applied here. 1.4. Bifurcation theory Bifurcation theory deals with the analysis of branch points of nonlinear equations in Banach space. Branch space is a real or complex vector space in which a vector has a non-negative length (or norm), and in which every Cauchy sequence converges to a point of the space. A key aspect of bifurcation theory is the implicit function theorem (Golubitsky & Schaeffer, 1985; Fujii & Ramm, 1997; Saha, Banerjee, & Chowdhury, 2001; Bari, 2002; L´opez-G´omez & Mora-Corral, 2004; and Gatermann & Hosten, 2005). By exploiting this theorem, one can determine if a set of equations has a unique solution. 1.5. Complex computation Cayley (1879) and Brolin (1966) studied the behavior of iterative maps of complex variables in mathematical physics. Many important contributions to the theory of nonlinear systems, such as bifurcation theory, have stemmed from their studies. Many of the important features of computations in the complex domain have found applications in fluid mechanics, chemically reacting systems, process control, and other areas of science and engineering. However, only recently (e.g. Allgower & Georg, 1992; and Andrei & Chin, 2004) has attention been directed to calculations in the complex domain for more general application to methods for solving general systems of nonlinear algebraic and transcendental equations. Recent studies that deal with multivariable maps, of interest to chemical engineering, in the complex domain include Lucia and Xu (1992), Lucia and Taylor (1992), Lucia, Guo, and Wang (1993), Lucia and Wang (1995), Taylor et al. (1996), Lucia (2000), and Gritton et al. (2001). To illustrate some of the basic details of the complex computation using homotopy continuation, consider the following cubic equation, presented by Choi and Book (1989): f (x) = x − x3 = −x(x − 1)(x + 1) = 0 Using a fixed-point homotopy, with g(x) = x − xo in Eq. (2), the homotopy path can not reach any of the three solutions in real space with x = 2 as the starting point. However, by extending the computation into complex space, all three solutions to the cubic equation are found as follows:The complex form of the
Fig. 1. Homotopy path for Choi-book Example using complex space.
equation is: f (z) = z − z3 = −z(z − 1)(z + 1) = 0,
where, z = x + iy
This equation in terms of z can now be solved by using the fixed-point homotopy method discussed by Gritton et al. (2001), which tracks the path in the complex domain, or, as implemented here, the equation can be split into two equations; one for the real part of the equation and the other for the imaginary part, as follows: F (1) = x(−x2 + 3y2 + 1) = 0 F (2) = y(y2 − 3x2 + 1) = 0 Using the fixed-point homotopy function, H(1) = tF (1) + (1 − t)(x − xo ) H(2) = tF (2) + (1 − t)(y − yo ) The computation is started from. x = 2, y = 0, The result is shown in Fig. 1, where it can be seen that transferring from real to complex space and vise versa takes place by a bifurcation branch shown as a bridge between the upper-left path and the lower-right path. The result is that all three roots located on the lower-right path are reached from the starting point located on the upper left path. 1.6. Application to phase equilibrium Many chemical process design problems, such as the design of distillation columns or alternative separation systems, involve phase equilibrium, where, using the criterion of equilibrium for a system under constraints of fixed pressure and temperature, the Gibbs function has to reach its minimum value. If it is assumed that each of the components, i, is present in each phase, the problem can be stated as follows. Given a set of feed compositions for a non-reacting multicomponent system at a given temperature and pressure, the number and states of phases existing at equilibrium as well as the composition and quantity of each phase is determined. Meth-
2336
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345
ods for calculating the equilibrium composition of a mixture of chemical species can be divided into two groups. 1- Equilibrium-constant methods 2- Free-energy-minimization methods A continuing discussion of the advantages and disadvantages of the two methods has ensued since 1960. Before the 1970s, chemical engineers relied mainly on Method 1, while engineers in the rocket industry relied on Method 2. Zeleznik and Gordon (1968) discussed the necessity of separate and organized technical discussions concerning these two methods. Today, in most computer programs that involve chemical equilibrium, minimization methods are the basis for phase-equilibrium calculations (e.g. Gordon & McBride, 1971, and Gautam & Seider, 1979a, 1979b, 1979c). For solving the equations of phase equilibria, there are several approaches used for Method 1. In one approach, Chavez, Seader, and Wayburn (1986) applied continuation from multiple starting points to find all real positive solutions to three examples of interlinked distillation systems where physical equilibrium was assumed at each stage. Lin, Seader, and Wayburn (1987) extended that work by using global homotopy continuation to find all real positive solutions from just one starting point by continuing the path until all roots were found. Heidemann (1983) reviewed advances in the computation of high-pressure phase equilibria using a successive-substitution algorithm. Also, Null (1970), Prausnitz, Eckert, Orye, and Oconnell (1980), and Anderson and Prausnitz (1980) discussed the practical implementation of a successive-substitution algorithm. Optimization approaches for Method 2 of the phaseequilibrium problem began with White, Johnson, and Dantzig (1958) and have been implemented with different methods since then, as mentioned above. The nonlinear objective function is the Gibbs free energy, which is subjected to a set of linear and nonlinear constraints (mass balances and non-negativity of moles). The mathematical properties of the nonlinear objective function depend completely on the mathematical form of the equation of state (EOS) or fugacity correlations chosen to represent each of the phases that may exist at equilibrium. These dependencies can be simple for the case of an ideal state or highly complicated for systems in which the NRTL (Non-Random Two-Liquid model) or UNIQUAC (UNIversal QUAsi-Chemical model) or an EOS near the critical region are used to represent nonideal liquid–liquid interactions. Nonconvexities can appear due to the form of the objective function and can lead to meta-stable equilibria or points whose multiphase solutions converge to a trivial solution. The global optimal solution of the mathematical formulation corresponds to the true physical equilibrium solution Therefore, equilibrium compositions are given at the global minimum of the Gibbs free energy surface. Poor initial guesses, however, for phase distribution and composition can lead to local minima. Constrained minima occur when too few phases are assumed, or, when the correct number of phases is assumed, but composition guesses are poor. Also local minima can occur when too many phases are assumed.
The probability of reaching local or constrained minima increases with the degree of nonideality and uncertainty in initial guesses for compositions. This occurs particularly when processes involving nonideal mixtures are analyzed using iterative procedures that produce large changes in compositions, thus increasing uncertainties in guessed compositions. Different calculations have been applied based on each of the criteria above and different researchers have used them for stability analysis of the phases at equilibrium (e.g. Shah, 1980; Ferraris & Morbidelli, 1981; Nagarajan, Cullick, & Griewank (1991a,b); Zhu, Wen, & Xu, 2000; Rangaiah, 2001; Balogh, Csendes, & Stateva, 2003; Burgos-Sol´orzano, Brennecke, & Stadtherr, 2004; and Xu, Haynes, & Stadtherr, 2005). Using the criterion of equilibrium for a system under constraints of fixed pressure and temperature, the Gibbs function reaches its minimum value. For a multiphase, multicomponent system, the Gibbs function is given by: G=
P C
¯ ij nij G
(3)
j=1 i=1
where nij is the number of moles of component i in phase j and ¯ ij is the partial molar Gibbs function or chemical potential, G which may be expressed as: ¯ ij = Goij + RT ln aij G
(4)
where Goij is the molar Gibbs function for component i in its standard state and aij is the activity of component i. The following constraints are present for phase equilibria: 1- Conservation of mass: nif =
P
nij ,
∀i
(5)
j=1
where the subscript, f, refers to the total for all phases. 2- Non-negativity: nij ≥ 0,
∀i, j
(6)
It has been assumed that each component, i, is present in each phase. The formulation can be written in terms of Lagrange multipliers to transform the problem to an unconstrained form. Activities, activity coefficients, fugacities, and fugacity coefficients are expressed in a variety of ways, depending upon the type of species and the type of mixture. The molar Gibbs free energy of mixing for a system of C components can be stated in a nonlinear programming problem (NLP) as: Minimize
G=
P C
nij Gij
i=1 j=1
=
P C i=1 j=1
nij (Goij + RT ln(γij xij )
(7)
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345
S.T.
h- (n- ) = nif =
P
nij ,
∀i
2337
(8)
j=1
and S- (n- ) = nij ≥ 0,
∀i, j
(9)
where γij = γij (n- j , T ) is the activity coefficient, which is a nonlinear function of composition and temperature, but is assumed not to be a function of pressure. The available thermodynamic models (empirical or semiempirical) for the Gibbs free energy may not be convex functions. Thus, the global minimum may or may not be found during the optimization process. It has been proved that a unique solution exists for a single ideal-solution phase system (Shapiro & Shapley, 1965). The nonuniqueness for nonideal systems has been discussed by many authors (Othmer, 1976; Caram & Scriven, 1976; Heidemann, 1978; Vandongen & Doherty, 1983; Harvie, Greenberg, & Weare, 1987; Castier, Rasmussen, & Fredenslund, 1989; Greiner, 1991; and Zhu et al., 2000). 1.7. Test problems for phase equilibrium In this section, two examples are presented of the continuation algorithm as applied to phase equilibrium in nonideal systems to show the ability of the proposed method to find all solutions, including the global equilibrium solution. Example 1. The first example is a modification of one used by Eubank and Barrufet (1988), where the van Laar activitycoefficient model is used, rather than the Margules equation for the Gibbs free energy of a binary mixture: G Ax1 x2 = + x1 ln x1 + x2 ln x2 RT Ax1 /B + x2
(11)
Values of 1.75 and 3.5 are used for A and B, respectively. Differentiation of the above equation and equating the result to zero yields the nonlinear equation: ∂(G/RT ) A(A/Bx22 − 1/x12 ) x1 =0 (12) =− + ln 2 ∂x1 x2 (A/Bx2 + 1/x1 ) In which the variables, x1 and x2 , are considered to be complex. In this work, Eq. (12) was combined with an equation for the sum of the mole fractions equaling 1 to reduce the system to one equation in x1 . Calculations were first carried out in real space only, using an initial guess for x1 of 0.8. The result is shown in Fig. 2, where only one solution is found. The problem was then computed in the complex domain with two equations in variables, x1 and y1 , with a starting point of (0.8, 0.0). As shown in Fig. 3, three solutions were found in the real domain at x1 equal to 0.2037 and 0.9646 at the minimum and 0.7103 at the maximum of the Gibbs free energy. These solutions are close to those of Eubank and Barrufet, who used the Margules equation for activity coefficients. Fig. 3 shows the real part of a complex domain that bridges the two real-domain paths, each of which contains a root for x1 in Eq. (12). The occurrence of complex-domain bridges are well known in homotopy continuation. In chapter 18 of Garcia and Zangwill (1981), a quadratic equation example is presented that involves
Fig. 2. Homotopy path for Example 1 using real space.
a real path that starts from x = 1 at t = 0, enters the complex domain at t = 0.168, re-enters the real domain at t = 0.952, and terminates at a real root of 3.0 at t = 1. More practical applications involving calculations with complex-domain bridges in order to find real roots for phase equilibria are presented by Lucia and Taylor (1992) and Gow, Guo, Liu, and Lucia (1997). However, a distinct difference exists between the type of bridge shown here and by Garcia and Zangwill and that of Lucia and Taylor. With homotopy continuation, the bridge occurs because of the change in the homotopy parameter, which changes the function(s) so that it (they) transition(s) from having only real roots to including complex roots. The type of bridge shown for the first time by Lucia and Taylor involves only the original function(s) and does not depend on parameterization. Example 2. In this second example, another binary phase equilibrium problem (based on Example 7.1(b) in the book by Walas, 1985) was solved using the van Laar activity-coefficient equation. In this example, when the problem is solved in the real domain with an initial value of 0.5 for x1 , only one solution
Fig. 3. Homotopy path for Example 1 using complex space.
2338
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345
the liquid phase is modeled using the NRTL or UNIQUAC equations. The Gibbs tangent-plane criterion for mixtures where the Gibbs free energy is a continuous first-order function is proved by Baker, Luks, and Pierce (1982). They demonstrated the problems associated with conventional stability criteria when several solutions exist that satisfy the condition of equal chemical potentials. Smith, Missen, and Smith (1993) obtained the necessary and sufficient conditions for stability of systems where chemical reaction may or may not be taking place, using the Karush–Kuhn–Tucker (KKT) optimality conditions as the basis of their consideration. 1.9. Test problem for simultaneous chemical and phase equilibrium Fig. 4. Homotopy path for Example 2 using real space.
is reached by the fixed-point homotopy branch, as shown in Fig. 4. When the problem is repeated in real and complex space with the same initial value (x1 = 0.5, y1 = 0.0), all three solutions (x1 = 0.1448, 0.5, and 0.8552) are found, as shown in Fig. 5. 1.8. Application to simultaneous chemical and phase equilibria Simultaneous phase and chemical equilibrium, which is of crucial importance in several process separation applications, is a difficult problem in chemical engineering. The true equilibrium state for constant temperature and pressure corresponds to the global minimum of the Gibbs free energy function. Because of the nonlinearity of the thermodynamic models and the possibility of multiple solutions, finding the global minimum is not guaranteed. However, McDonald and Floudas (1995b) present a global optimization approach that will find a global solution for the minimization of the Gibbs free energy function, when
Example 3. To illustrate the application of the homotopycontinuation method to simultaneous chemical and phase equilibrium, a problem has been solved involving a chemical reaction in the vapor phase in the presence of a liquid phase, as discussed by Xiao, Zhu, Yuan, & Chien, 1989). The feed consists of an equimolar solution of ethanol and acetic acid, which undergoes the esterification reaction: Ethanol + Aceticacid → Ethylacetate + Water Data required for the chemical-equilibrium constant are given by Sanderson and Chien (1973). At a system temperature of 358 K and a pressure of one atmosphere, which was used here, a liquid phase will co-exist with a vapor phase. Using the criterion of chemical and phase equilibrium for a system under constraints of fixed pressure and temperature, the Gibbs function has to attain its minimum value. For a multiphase, multicomponent system, the Gibbs function is given by: G=
P C
¯ ij nij G
(13)
j=1 i=1
where nij is the number of moles of component i in phase j and ¯ ij is the partial molar Gibbs function or chemical potential, G which may be expressed as: ¯ ij = Goij + RT ln aij G
(14)
where Goij is the molar Gibbs function for component i in its standard state and aij is the activity of component i. The following constraints apply for simultaneous chemical and phase equilibria: 1- Conservation of mass: P C
Aik nij = bk
k = 1, . . . , E
(15)
j=1 i=1
2- Non-negativity: Fig. 5. Homotopy path for Example 2 using complex space.
nij ≥ 0
∀i, j
(16)
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345
where Aik is the coefficient of element k in component i, E is the number of elements in the system, and bk is the total amount of atoms of element k in the system. It is assumed that each of the components i is present in each phase. Of course, in the case of pure condensed phases, this is not true and modifications to the equations are necessary (Eriksson & Rosen, 1973). Smith (1980) has shown that a number of equivalent formulations can be made which are useful in different situations. A second formulation can be written in terms of Lagrange multipliers to transform the problem to an unconstrained form as follows. The molar Gibbs free energy of mixing for a system of C components can be stated in a nonlinear programming problem (NLP) as:
Minimize G =
P C
2339
Table 1 Binary interaction parameters for Wilson equation System
A12
A21
EtOH–HAC EtOH–EtAC EtOH–H2 O HAC–EtAC HAC–H2 O EtAC–H2 O
101.6588 288.2011 198.1757 1749.9343 2.0311 26981.1421
−130.6527 28.879 466.1059 −464.1592 403.1564 1195.67
¯ ij nij G
i=1 j=1
=
P C
[nij Goij + RT ln(γij xij )]
(17)
i=1 j=1
¯ n) = S.T. h(¯
P C
Aik nij − bk = 0
(18)
j=1 i=1
and nij λij = 0
∀i, j
(19)
where γij = γij (n- ij , T ) is the activity coefficient, which is a nonlinear function of composition and temperature, and λij are Kuhn–Tucker multipliers. The global minimum of the Gibbs free energy during conventional optimization processes may or may not be found since the available thermodynamic models (empirical or semi-empirical) may not be convex functions. The non-uniqueness for nonideal systems has been discussed by many authors (e.g. Othmer, 1976; Caram & Scriven, 1976; Heidemann, 1978). In this work, the homotopy-continuation method for solving an optimization problem is used to develop a computerized procedure for predicting phase equilibria and for finding the global minimum by searching through real and complex space. In this example, liquid-phase activity coefficients, γ, were correlated by the Wilson equation, while the vapor phase was assumed to be ideal. The homotopy functions were based on Eqs. (17)–(19), with the functions converted from real to complex variables. A fixed-point homotopy function was applied to the set of equations, which were tracked with the AUTO program
Fig. 6. Homotopy path for Example 3 using real space.
of Deodel (1986), which applies continuation and bifurcation calculations. The binary interaction parameters for the Wilson equation are given in Table 1. Three methods were applied to solve Example 3: Newton–Raphson, homotopy continuation in real space, and homotopy continuation in real and complex space. The Newton–Raphson method did not converge. The results for the homotopy-continuation methods are shown in Table 2. Homotopy continuation in real space reached only one solution, while homotopy continuation in real and complex space found two sets of solutions, one of which is the desired solution. Figs. 6 and 7 show the homotopy branches for real space and complex domain calculations, respectively for the mol% water in the vapor phase. As can be seen in Fig. 7, continuation started from 40 mol% water in the vapor phase, moved to the right until a limit point was reached and then moved to the left. However, the limit point was also a bifurcation point, from which a real part moved to
Table 2 Results of homotopy-continuation calculations in real and complex space Homotopy in real space
EtOH HAC EtAC H2O
Homotopy in complex space
Liquid (mol%)
Vapor (mol%)
Liquid, 1st set (mol%)
Vapor 1st set (mol%)
Liquid, 2nd set (mol%)
Vapor, 2nd set (mol%)
5.294 22.6 12.44 59.64
8.803 6.52 45.45 39.22
4.82 20.84 11.49 56.18
8.65 6.47 45.48 39.39
5.29 22.6 12.44 59.64
8.8 6.52 45.45 39.22
2340
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345
Fig. 7. Homotopy path for Example 3 using complex space.
the right until another limit and bifurcation point was reached. From there, two branches continued in real space, with two solutions determined when the homotopy parameter reached values of one. The second set of results in Table 2 is in close, but not exact agreement with the results of Xiao et al. (1989), who used a slightly lower temperature of 355 K and the UNIQUAC equation for estimating liquid-phase activity coefficients. Neither that study nor this one allowed for nonideality and dimerization of acetic acid in the vapor phase or possible formation of a second liquid phase, as in the work of Castier et al. (1989). We plan to consider this in future work. 1.10. Application to stability analysis Calculation of equilibrium using thermodynamic stability analysis has found considerable interest during the past several years, with some applications of the stability test for chemical and phase equilibria computations reported in the literature. Sander, Rasmussen and Fredenslund (1986) apply the test to solid–liquid equilibria of aqueous salt solutions. Li, Ngheim, and Heidemann (1986) apply the stability test to VLE calculations using a modified form of the Peng and Robinson (1976) equation of state (EOS) to consider chemical association in both phases. Cairns and Furzer (1990) consider the phase stability of multicomponent, three-phase azeotropic distillation. Sun and Seider (1995) apply a homotopy-continuation method to the stability analysis of phase equilibria. McDonald and Floudas (1995a) present a global-optimization method for phase stability problems involving activity coefficients. Jiang, Chapman and Smith (1996) classify all possible types of behavior in multiphase, multi-reaction systems and introduce the reaction tangent plane. McKinnon and Mongeau (1998) develop algorithms to utilize the reaction tangent plane for stability analysis, and Jalali and Seader (2000) apply homotopy continuation to stability of multiphase, multi-reacting systems. The Gibbs tangent-plane criterion has become a powerful tool in solving equilibrium and stability problems. The ability to determine if a solution is stable with respect to perturbations in
any or all of the phases are important and useful in the search for the true equilibrium solution. Many approaches have been developed to find the stationary points of the tangent-plane-distance function, but due to the complex, nonlinear nature of the models used to predict equilibrium, none has guaranteed to find all stationary points. McDonald and Floudas (1995a) applied formulations of the stability criteria for a special class of problems where non-ideal liquid phases were modeled by the NRTL and UNIQUAC activity-coefficient equations. They showed how the global minimum of the tangent-plane-distance function can be obtained for that class of problems. Michelsen (1982a, 1982b) used the tangent-plane criterion efficiently for the first time. A two-stage approach was proposed whereby the stability problem was used to generate initial points to be used in the search for a global minimum of the Gibbs free energy. Experience revealed that if the postulated solution is indeed thermodynamically unstable, the solutions obtained from the stability problem usually provide a good starting point for the new search. Sun and Seider (1995) also applied a twostage approach with homotopy continuation to find the global solutions. Nagarajan et al. (1991a,b) reformulated the stability approach, replacing the mol numbers by mol densities, so that the Helmholtz free energy becomes the proper thermodynamic function to describe equilibrium. To avoid the singular Jacobian problem, Gupta, Bishnoi, and Kalogerkis (1991) combined the phase equilibrium and stability problem. Their approach is guaranteed to converge to a stationary point. Even though, the stability problem has been effectively used to improve the chances of finding the true equilibrium solution corresponding to a global minimum of the Gibbs free energy, there is no guarantee. It is usually assumed that if a postulated solution is found to be thermodynamically unstable, then a phase must be added or deleted. However, the true equilibrium solution may in fact be a solution with the same number of phases as the postulated one, but with a different distribution of the components within those phases. 1.10.1. Michelsen criteria The method of Michelsen (1982a, b) utilizes the tangentplane criterion illustrated by Baker et al. (1982). Based on this criterion, the equilibrium tie-line between the coexisting liquid phases must be tangent to the curve of the Gibbs free energy versus mole fraction, as in Fig. 8, at both α and β, in order to have equality of the chemical potential of species in the two phases. Here, a brief description of the method used by Michelsen is presented. The Gibbs free energy of a mixture of C components at a given temperature and pressure with mole numbers ni is: G=
C
n i μi
(20)
i=1
where μi is the chemical potential of component i. It is assumed that the mixture is split into two phases, I and II, with mole numbers (NT − ε) and ε, respectively, where ε is a very small amount of the second phase. Assuming the mole fraction of the
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345
2341
Fig. 9. Molar Gibbs free energy of mixing for a binary mixture showing the tangent at mole fraction z and tangent distance F at y.
tionary condition as: μi (y) − μi = K, Fig. 8. The dimensionless change in Gibbs free energy at constant temperature and pressure.
second phase is yi , the change in Gibbs free energy is: ΔG = GI + GII − G = G(NT − ε) + G(ε) − G
(21)
A Taylor series expansion of Gi yields: C C ∂G = G − ε yi μ i G(N − ε) = G(N) − ε yi ∂ni N i=1
(22)
Therefore,
i=1
(23)
i=1
where μi is the chemical potential for ni and μi (y) is the chemical potential for yi in ε. Stability of the original mixture requires that its Gibbs free energy be at the global minimum. Therefore, a necessary condition for stability is that F (y) =
C
yi [μi (y) − μi ] ≥ 0
(25)
where K is independent of the component i. Based on Eq. (23) where the Gibbs energy should be a minimum at equilibrium, the original system is stable if K is positive or equal to zero. Using the NRTL equation for chemical potential the above equation becomes: F (y) = yi [ln yi + ln γi − ln zi − ln γi (z)] RT C
(26)
i=1
i=1
C C ΔG = G(ε) − ε yi μi = ε [yi μi (y) − μi ]
∀i
(24)
i=1
for all trial compositions y. Geometrically, F(y) is the vertical distance from the tangent hyperplane to the molar Gibbs energy surface at composition z (mole fraction of original mixture) to the energy surface at composition y, as shown in Fig. 9 for a binary mixture. Because all minima of F(y) are located in the region (yi ≥ 0,
yi = 1), F(y) is positive if it is positive at stationary points. Differentiation with respect to the mole fractions gives the sta-
where γ i is the activity coefficient of component i, y is the mole fraction vector of the trial phase, and z is the mole fraction vector of the original system. The stationary condition (where the derivative F(y)/RT with respect to yi is equal to zero) is as follows: ln yi + ln γi − ln zi − ln γi (z) = K
i = 1, . . . , C
(27)
Introducing new variable, Yi = exp(−K)yi , in order to eliminate K from the computation, the above equation becomes: ln Yi + ln γi (y) − ln zi − ln γi (z) = 0
i = 1, . . . , C
(28)
where Michelsen interprets the new variable Yi as a mole number. The stability is verified if at all stationary points K ≥ 0, which means that Yi ≤ 1, and conversely, the phase is unstable if a stationary point is located where Yi > 1. Also, Michelsen devised an algorithm for locating the minimum vertical distance from the tangent hyperplane to the surface of the Gibbs free energy of mixing using either a direct iteration scheme or a minimization technique. The method also provides an initial guess of the phase-split for an unstable phase, and either a free-energy-minimization method or the direct-iteration method of Henley and Rosen (1969) can be used to determine the equilibrium compositions. This initial guess insures that the molar Gibbs free energy will have a lower value than that of the original single-phase system.
2342
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345
Michelson’s stability methods also require an initial guess. He recommends using as many initial guesses as there are components. Each initial guess for y is taken to be a pure component and is converged simultaneously for a maximum of four direct iterations. If one of the initial guesses is dropped from consideration, the initial guess with the negative F(y) is selected for final convergence. Here, Michelson’s method is combined with the homotopycontinuation method for stability analysis. When homotopy continuation is restricted to real space for the stability test, just some of roots are found from an initial guess. Alternatively, homotopy branches in real and complex space are connected to each other through bifurcation branches. Therefore, multiple solutions are found by just one initial guess, as illustrated in the following test problem. Fig. 10. Homotopy branch in real space for Example 4.
1.11. Test problem for stability analysis Example 4. In this problem, a mixture of water (1) and trimethylamine (2) in equal molar amounts and at atmospheric pressure is brought to phase equilibrium with two liquid phases. The NRTL equation is used to model the activity coefficients, with the following binary interaction coefficients: g12 = 2072.9748, g21 = 2058.3746, a12 = a21 = 0.4144. Initial mole numbers are set equal to 0.5. To check the stability, the following system of equations is solved: F (Yi ) = ln(Yi ) + ln γi (y) − ln(zi ) − ln γi (z) = 0
(29)
A fixed-point homotopy is used here: tF (Y ) + (1 − t)(Y − Y o ) = 0
(30)
For a binary system, the equations that must be solved are: ln(Y1 ) + ln γ1 (y) − ln(z1 ) − ln γ1 (z) = 0 ln(Y2 ) + ln γ2 (y) − ln(z2 ) − ln γ2 (z) = 0 Y1 y1 = Y1 + Y 2 Y2 y2 = Y1 + Y 2
Fig. 11. Homotopy branch in complex space for Example 4.
(31)
In the above system, zi is specified and ln γ i (z) is calculated. The unknown variables are Y1 and Y2 . In this work, the equation was first calculated in real space with the initial guess Y1 = 1 and Y2 = 0.5. The result is shown in Fig. 10. As can be seen, only one solution (u1 ) is found. The other roots are on another branch, with limit points (LP) on the first branch so that we cannot reach the second branch by using continuation only in real space because there is a complex region between the two branches. Therefore, the calculations were repeated in both the real and complex domain by considering Y as a complex variable Z. Now, as shown in Fig. 11, the continuation method is capable of tracing the first branch in complex space from branching point (BP) to reach the other branch in real space. The computation is done with the same initial guess Y1 real = 1, Y1 complex = 0, Y2 real = 0.5, and Y2 complex = 0. All solutions are traced using homotopy branches in both domains. A complex Gibbs energy
Table 3 Results of stability analysis for Example 4 Roots
Y1
Y2
Yi
Stability
1 2 3 4 5
0.959050 0.500000 0.259498 0.052932 0.734975
0.050563 0.500000 0.737483 0.955488 0.262072
1.009613 1.000000 0.996981 1.008420 0.997047
Unstable Trivial Stable Unstable Stable
is never used in the calculations. The five roots of Eq. (31) and the results of stability analysis are in Table 3.A plot of the Gibbs free energy change is shown in Fig. 12. 2. Discussion In this research, calculations of multiple solutions for phase equilibria, simultaneous phase and chemical equilibria, and phase stability analysis are presented using a homotopy continuation method in the complex plane with bifurcation branches
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345
Fig. 12. Change in Gibbs energy for Example 4.
as a connecting tool between two real spaces. For the above problems, a fixed-point homotopy was applied because it is sensitive to the initial guess. The fixed-point homotopy can have a scaling problem, and in the above examples, scaling is an important factor in the computations. However, it did not prevent us from obtaining the multiple solutions. Complexification cannot be applied if the equations cannot be converted into complex variables. In that case, there are some modifications available to overcome the difficulties. In order to prevent a singularity during the computations, it is possible to increase the precision of the iteration, so the calculations can pass through any singular points. Acknowledgment Part of this research has been supported by a grant from the University of Tehran, Faculty of Engineering. References Allgower, E., & Georg, K. (1980). Simplicial and continuation methods for approximating fixed points and solutions to systems of equations. SIAM Review, 22, 28. Allgower, E. L., & Georg, K. (1992). Continuation and path following. Acta Numerica, 1–64. Anderson, T. F., & Prausnitz, J. M. (1980). Computational methods for highpressure phase equilibria and other fluid-phase properties using a partition function. Part I. pure fluids. Industrial & Engineering Chemistry Product Research and Development, 19, 1. Andrei, S., & Chin, W. (2004). Solving a class of higher-order equations over a group structure. Journal of Symbolic Computation, 37, 329. Baker, L. E., Luks, K. D., & Pierce, A. C. (1982). Gibbs energy analysis of phase equilibria. Society of Petroleum Engineers Journal, 22, 731. Balogh, J., Csendes, T., & Stateva, R. P. (2003). Application of a stochastic method to the solution of the phase stability problem: cubic equations of state. Fluid Phase Equilibria, 212, 257. Bari, R. (2002). Local bifurcation theory for multiparameter nonlinear eigenvalue problems. Nonlinear Analysis, 48, 1077. Bausa, J., & Marquardt, W. (2000). Quick and reliable phase stability test in VLLE flash calculations by homotopy continuation. Computers & Chemical Engineering, 24, 2447.
2343
Brolin, H. (1966). Invarient sets under Iteration of rational functions. Arkiv for Matematik, 6, 103. Burgos-Sol´orzano, G. I., Brennecke, J. F., & Stadtherr, M. A. (2004). Validated computing approach for high-pressure chemical and multiphase equilibrium. Fluid Phase Equilibria, 219, 245. Cairns, B. P., & Furzer, I. A. (1990). Multicomponent three-phase azeotropic distillation. 2. Phase-stability and phase-splitting algorithms. Industrial & Engineering Chemistry, Research, 29, 1364. Caram, H. S., & Scriven, L. E. (1976). Non-unique reaction equilibria in nonideal systems. Chemical Engineering and Science, 31, 163. Castier, M., Rasmussen, P., & Fredenslund, A. (1989). Calculation of simultaneous chemical and phase equilibria in nonideal systems. Chemical Engineering and Science, 2, 37. Cayley, A. (1879). The Newton–Fourier imaginary problems. American Journal of Mathematics II, 97. Chavez, C. R., Seader, J. D., & Wayburn, T. L. (1986). Multiple steady-state solutions for interlinked separation systems. Industrial & Engineering Chemistry Fundamentals, 25, 566. Choi, S. H., & Book, N. L. (1989). An evaluation of the starting point criteria for the Global fixed-point homotopy continuation method. Department of Chemical Engineering, University of Missouri-Rolla. Chow, S. N., Mallet-Paret, J., & Yorke, J. A. (1978). Finding zeroes of maps: Homotopy methods that are constructive with probability one. Mathematics of Computation, 32, 887. Deodel, E. (1986). AUTO 86 User manual, software for continuation and bifurcation problems in ordinary differential equations. In Applied Mathematics. Pasadena, California: California Institute of Technology. Eubank, P. T., & Barrufet, M. A. (1988). A simple algorithm for calculation of phase separation. Chemical Engineering and Education, 36–41. Eriksson, G., & Rosen, E. (1973). Thermodynamic studies of high temperature equilibria VIII. General equations for the calculation of equilibria in multiphase systems. Chemica Scripta, 4, 193. Ferraris, G. B., & Morbidelli, M. (1981). Distillation models for two partially immiscible liquids. American Institute of Chemical Engineering Journal, 27, 881. Fujii, F., & Ramm, E. (1997). Computational bifurcation theory: Path-tracing, pinpointing and path-switching. Engineering Structures, 19, 385. Garcia, C. B., & Zangwiil, W. I. (1979). Determining all solutions to certain systems of nonlinear equations. Mathematical of Operations Research, 4, 1. Garcia, C. B., & Zangwill, W. I. (1981). Pathways to solutions, fixed points, and equilibria. New Jersey: Prentice-Hall, Englecliff Cliffs. Gatermann, K., & Hosten, S. (2005). Computational algebra for bifurcation theory. Journal of Symbolic Computation, 40, 1180. Gautam, R., & Seider, W. D. (1979a). Calculation of phase and chemical equilibria, Part I, local and constrained minima in Gibbs free energy. AIChE Journal, 25, 991. Gautam, R., & Seider, W. D. (1979b). Calculation of phase and chemical equilibria, Part II: phase-splitting. AIChE Journal, 25, 999. Gautam, R., & Seider, W. D. (1979c). Calculation of phase and chemical equilibria; Part III: electrolytic solutions. AIChE Journal, 25, 1007. Golubitsky, M., & Schaeffer, D. G. (1985). Singularities and groups in bifurcation theory, Vol. 1. NY: Springer-Verlag. Gordon, S., McBride, B.J. (1971). Computer Program for Calculation of Complex Chemical Equilibrium Compositions, Rocket Performance, Incident and Reflected Shocks, and Chapman-Jouguet Detonations. NASA SP-273. Gow, A. S., Guo, X., Liu, D., & Lucia, A. (1997). Simulation of refrigerant phase equilibria. Industrial Engineering and Chemical Research, 36, 2841. Greiner, H. (1991). An efficient implementation of Newton’s method for complex nonideal chemical equilibria. Computers & Chemical Engineering, 15, 115. Gritton, K. S., Seader, J. D., & Lin, W. (2001). Global homotopy continuation procedures for seeking all roots of a nonlinear equation. Computers & Chemical Engineering, 25, 1003. Gupta, A. K., Bishnoi, P. R., & Kalogerkis, N. (1991). A method for the simultaneous phase equilibria and stability calculations for multiphase reacting and non-reacting systems. Fluid Phase Equilibria, 65, 63.
2344
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345
Harvie, C. E., Greenberg, J. P., & Weare, J. H. (1987). A chemical equilibrium algorithm for highly non-ideal multiphase systems: Free energy minimization. Geochimica et Cosmochimica Acta, 51, 1045. Heidemann, R. A. (1978). Non-uniqueness in phase and reaction equilibrium computations. Chemical Engineering Science, 33, 1517. Heidemann, R. A. (1983). Computation of high pressure phase equilibria. Fluid Phase Equilibria, 14, 55. Henley, E. J., & Rosen, E. M. (1969). Material and energy balance computation. New York: Wiley. Jalali, F., & Seader, J. D. (1999). Homotopy continuation method in multi-phase multi-reaction equilibrium systems. Computers & Chemical Engineering, 23, 1319. Jalali, F., & Seader, J. D. (2000). Use of homotopy-continuation method in stability analysis of multiphase, reacting systems. Computers & Chemical Engineering, 24, 1997. Jiang, Y., Chapman, G. R., & Smith, W. R. (1996). On the geometry of chemical reaction and phase equilibria. Fluid Phase Equilibria, 118, 77. Li, Y. K., Ngheim, L. X., & Heidemann, R. A. (1986). Investigation of phase behavior predictions with chemical association equation of state. Report CMG.R12.11, University of Calgary, Canada. Li, T. Y. (1987). Solving polynomial systems. Mathematical Intelligence, 9(3), 33. Li, T. Y. (2003). Numerical solution of polynomial systems by homotopy continuation methods. Handbook of Numerical Analysis, 11, 209. Lin, W. J., Seader, J. D., & Wayburn, T. L. (1987). Computing multiple solutions to systems of interlinked separation columns. AIChE Journal, 33(6), 886. L´opez-G´omez, J., & Mora-Corral, C. (2004). Minimal complexity of semibounded components in bifurcation theory. Nonlinear Analysis, 58, 749. Lucia, A., & Xu, J. (1992). Global minima in root finding. In C. A. Floudas & P. M. Pardalos (Eds.), Recent Advances in Global Optimization. Princeton, NJ: Princeton University Press. Lucia, A., & Taylor, R. (1992). Complex iterative solutions to process model equations? Computers and Chemical Engineering, 16S, S387. Lucia, A., Guo, X., & Wang, X. (1993). Process simulation in the complex domain. AIChE Journal, 39, 461. Lucia, A., & Wang, X. (1995). Complex domain process dynamics. Industrial Engineering & Chemistry Research, 34, 202. Lucia, A. (2000). Complex domain chemical process simulation in theory and in practice. Industrial & Engineering Chemistry Research, 39, 1713. Lucia, A., & Feng, Y. (2002). Global terrain methods. Computers & Chemical Engineering, 26, 529. Lucia, A., & Feng, Y. (2003). Multivariable terrain methods. AIChE Journal, 49, 2553. McDonald, C. M., & Floudas, C. A. (1995a). Global optimization for the phase stability problem. American Institute of Chemical Engineering Journal, 41, 1798–1814. McDonald, C. M., & Floudas, C. A. (1995b). Global optimization for the phase and chemical equilibrium problem: Application to the NRTL equation. Computers and Chemical Engineering, 19, 1111. McKinnon, K., & Mongeau, M. (1998). A generic global optimization algorithm for the chemical and phase equilibrium problem. Journal of Global Optimization, 12, 325. Michelsen, M. L. (1982a). The isothermal flash problem: Part I: stability. Fluid Phase Equilibria, 9, 1. Michelsen, M. L. (1982b). The isothermal flash problem. Part II. phase-split calculation. Fluid Phase Equilibria, 9, 21. Morgan, A. (1987). Solving polynomial systems using continuation for engineering and scientific problems. NJ: Prentice-Hall. Nagarajan, N. R., Cullick, A. S., & Griewank, A. S. C. (1991a). New strategy for phase equilibrium and critical point calculations by thermodynamic energy analysis. Part I. stability analysis and flash. Fluid Phase Equilibria, 62, 191. Nagarajan, N. R., Cullick, A. S., & Griewank, A. S. C. (1991b). New strategy for phase equilibrium and critical point calculations by thermodynamic energy
analysis. Part II. Critical point calculations. Fluid Phase Equilibria, 62, 211. Null, H. R. (1970). Phase equilibrium in process design. NY: Wiley-InterScience. Othmer, H. G. (1976). Nonuniqueness in equilibria in closed reacting systems. Chemical Engineering Science, 31, 993. Peng, D. Y., & Robinson, D. B. (1976). A new two constant equation of state. Industrial & Engineering Chemistry Fundamentals, 15, 59. Prausnitz, J. M., Eckert, C. A., Orye, R. V., & Oconnell, J. P. (1980). Computer calculations for multicomponent vapor–liquid and liquid–liquid equilibria. NJ: Prentice-Hall. Rangaiah, G. P. (2001). Evaluation of genetic algorithms and simulated annealing for phase equilibrium and stability problems. Fluid Phase Equilibria, 187–188, 83. Saha, P., Banerjee, S., & Chowdhury, A. R. (2001). Bifurcation and chaos in a system simulating magnetic field configurations. Chaos, Solitons & Fractals, 12, 2247. Sander, Rasmussen, B. P., & Fredenslund, A. (1986). Calculation of solid–liquid equilibria in aqueous solutions of nitrate salts using an extended UNIQUAC equation. Chemical Engineering Science, 41, 1197. Sanderson, R. V., & Chien, H. H. Y. (1973). Simultaneous chemical and phase equilibrium calculation. AIChE Journal, 12(1), 81–85. Seader, J. D., Kuno, M., Lin, W. J., Johnson, S. A., Unsworth, K., & Wiskin, J. W. (1990). Mapped continuation methods for computing all solutions to general systems of nonlinear equations. Computers & Chemical Engineering, 14(1), 71. Shapiro, N. Z., & Shapley, L. S. (1965). Mass action laws and the Gibbs free energy function. Journal Social & Industrial Applied Magazine, 13, 353. Shah, V. B. (1980). Multicomponent Distillation with Two Liquid Phases, Ph.D. Thesis, University of Toledo. Smith, W. R. (1980). The calculation of chemical equilibria in complex systems. Industrial Engineering & Chemistry Fundamentals, 19, 1. Smith, J. V., Missen, R. W., & Smith, W. R. (1993). General optimality criteria for multiphase multireaction chemical equilibria in complex systems. AIChE Journal, 39, 707. Sun, A. C., & Seider, W. D. (1995). Homotopy—continuation method for stability analysis in the global minimization of the Gibbs free energy. Fluid Phase Equilibria, 103, 213. Taylor, R., Achuthan, K., & Lucia, A. (1996). Complex domain distillation calculations. Computers and Chemical Engineering, 20, 93. Vandongen, D. B., & Doherty, M. F. (1983). Calculation and stability of multiple equilibrium points for nonideal mixtures. Fluid Phase Equilibria, 14, 335. Walas, S. M. (1985). Phase equilibria in chemical engineering. Stoneham, MA: Butterworth. Watson, L. T. (1986). Engineering applications of the Chow-York algorithm. Applied Mathematics and Computation, 9, 111. Watson, L. T., Billups, S. C., & Morgan, A. P. (1987). Algorithm 652, HOMPACK: A suite of codes for globally convergent homotopy algorithm. ACM Transactions of Mathematical Software, 13, 281. Watson, L. T., Sosonika, M., Melville, R. C., Morgan, A. P., & Walker, H. F. (1997). Algorithm 777: HOMPACK90, a suite of Fortran 90 codes for globally convergent homotopy algorithms. ACM Transactions on Mathematical Software, 23, 514. Watson, L. T. (2000). Theory of globally convergent probablility-one homotopies for nonlinear programming. SIAM Journal of Optimization, 11, 761. Wayburn, T. L., & Seader, J. D. (1987). Homotopy continuation methods for computer-aided process design. Computers & Chemical Engineering, 11, 7. White, W. B., Johnson, S. M., & Dantzig, G. B. (1958). Chemical equilibrium in complex mixtures. Journal of Chemistry and Physics, 28, 751. Wise, S.M., Sommese, A.J., & Watson, L.T. (2000). Algorithm 801: POLSYS PLP: A partitioned linear product homotopy code for solving polynomial systems of equations. ACM Transactions of Mathmatical Software, 26, 176. Wu, T. M. (2005). A study of convergence on the Newton-homotopy continuation method. Applied Mathematics and Computation, 168, 1169.
F. Jalali et al. / Computers and Chemical Engineering 32 (2008) 2333–2345 Xiao, W., Zhu, K., Yuan, W., & Chien, H. H. (1989). An algorithm for simultaneous chemical and phase equilibrium calculation. AIChE Journal, 35(11), 1813–1820. Xu, G., Haynes, W. D., & Stadtherr, M. A. (2005). Reliable phase stability analysis for asymmetric models. Fluid Phase Equilibria, 235, 152.
2345
Zeleznik, F. J., & Gordon, S. (1968). Calculation of complex chemical equilibria. Industrial & Engineering Chemistry, 60(6), 27. Zhu, Y., Wen, H., & Xu, Z. (2000). Global stability analysis and phase equilibrium calculations at high pressures using the enhanced simulated annealing algorithm. Chemical Engineering Science, 55, 3451.