Homotopy continuation method in multi-phase multi-reaction equilibrium systems

Homotopy continuation method in multi-phase multi-reaction equilibrium systems

Computers and Chemical Engineering 23 (1999) 1319 – 1331 www.elsevier.com/locate/compchemeng Homotopy continuation method in multi-phase multi-reacti...

162KB Sizes 4 Downloads 133 Views

Computers and Chemical Engineering 23 (1999) 1319 – 1331 www.elsevier.com/locate/compchemeng

Homotopy continuation method in multi-phase multi-reaction equilibrium systems F. Jalali a, J.D. Seader b,* b

a Department of Chemical Engineering, Uni6ersity of Tehran, Tehran, Iran Department of Chemical and Fuels Engineering, Uni6ersity of Utah, Salt Lake City, UT 84112, USA

Received 5 October 1998; received in revised form 17 September 1999; accepted 21 September 1999

Abstract In nonideal chemical systems in which chemical reactions occur in more than one phase, more than one equilibrium condition may be calculated, but only one of them will be stable. In this research, a nonlinear optimization is applied to locate all minima of the Gibbs free energy function, from which the global stable solution can be ascertained. To find multiple solutions to the chemical and phase equilibria (CPE) problems in nonideal systems homotopy continuation is applied, accounting for bifurcation phenomena. Solution branches are traced by pseudo-arclength continuation. © 1999 Elsevier Science Ltd. All rights reserved. Keywords: Homotopy continuation method; Multi-phase; Multi-reaction; Equilibrium

1. Introduction Over the past three decades, many computer-aided methods for process design and control have evolved for application to the chemical industry. Today, chemical processes are modeled with increasing complexity. Among the fundamental tasks in computer-aided design of chemical processes is that of solving the largescale nonlinear algebraic equation systems involved in steady-state modeling (Zitney & Stadtherr, 1988; Seader, Kuno, Lin, Johnson, Unsworth & Wiskin, 1990). There are two reasons for this interest. The first comes from the broadly recognized need to automate the task of numerical computing. The second reason is that a reliable method for solving the required nonlinear equation systems still needs to be developed. Most available software sometimes experiences difficulty in converging iterative calculations. Problems related to chemical process design and control can often be formulated as nonlinear optimization problems. These problems usually have non-convexities that imply difficulties in determining a global solution. Chemical engineering examples of such prob* Corresponding author. E-mail address: [email protected] (J.D. Seader)

lems are mostly in reactor network synthesis, phase and/or chemical reaction equilibrium, heat exchanger network design, batch processes, and optimal design of distillation. In phase equilibrium with or without chemical reaction, a feed composition at a specified temperature and pressure is given. The number and state of phases existing at equilibrium, as well as the composition and quantity of each phase, is determined by one of the several approaches to the calculation of equilibrium. Optimization approaches for the calculation of phase equilibrium began with White, Johnson and Dantzig (1958) and has been implemented with different techniques since then. 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). Non-convexities appear due to the form of the objective function and lead to meta-stable equilibria or points whose multiphase solution converges to a trivial solution. The global optimal solution of the mathematical formulation corresponds to the true equilibrium solution. The mathematical properties of the nonlinear objective function depend on the mathematical form of the equation of state (EOS) or activity-coefficient correlation chosen to represent each of the phases that may exist at equilibrium.

0098-1354/99/$ - see front matter © 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 9 9 ) 0 0 2 9 4 - X

1320

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331

These dependencies can be simple for the case of ideal solutions or highly complex for systems in which the NRTL (Non-Random Two-Liquid model), UNIQUAC (UNIversal QUAsi-Chemical model), or an EOS near the critical region is used to represent nonideal liquid–liquid interactions. Equilibrium compositions are given at the global minimum of the Gibbs free energy surface. For steady-state analysis, the resulting nonlinear equation system, i.e. f( (x¯,p¯ )=0

(1)

is solved, often using Newton-based methods. Here, f is a function of general variable, x, and parameters, p. The overbar indicates a vector. Methods for finding the zeros of Eq. (1) can be divided into local methods and global methods. Local methods only converge from a good initial guess. Global methods can converge from any initial guess. Homotopy continuation methods are global methods because they utilize a global property of the mapping that is preserved under homotopy, namely, topological degree. Nonlinear equations, however, can have multiple roots and it is highly desirable in many chemical engineering problems to find all real solutions. Methods for finding additional roots are discussed by Allgower and Georg (1980) and include the use of multiple starting points, deflation, global Newton methods, global homotopy methods, and a d-trick homotopy developed by them. By far the most successful of these approaches is limited to systems of polynomial equations (Chow, Mallet-Paret & Yorke, 1978; Garcia & Zangwill, 1979; Watson, 1986; Li, 1987; Morgan, 1987; Watson, Billups & Morgan, 1987). With the development of advanced computers there has been a great interest in developing algorithms that locate the global optimum. Extensive surveys and books on the approaches for global optimization are available (e.g. Archetti & Schoen, 1984; Torn & Zlinskas, 1987; Ratschek & Rokne, 1988; Hansen, Jaumard & Lu, 1989; Horst & Tuy, 1990; Floudas & Pardalos, 1990). They utilize and develop global convergence methods, algorithms for parameterization and stability analysis, methods for locating singularities and bifurcation points, and algorithms for tracking periodic solutions and locating secondary bifurcations. Their methods also utilize singularity and bifurcation theories, with expansions in the vicinity of singular points. However, no one has presented a general technique for finding the global solution. In this research, an attempt is made to reach a global solution for equilibrium between phases with and without reaction by combining new methods from applied mathematics. In this work, the number of phases is known as a priori. The stability of the

phases is discussed in a separate forthcoming paper. In the following sections, the mathematical bases for continuation, nonlinear optimization, and bifurcation analysis, are first discussed. Then applications are made to phase equilibria (PE) and CPE.

2. Continuation method The continuation method is a major advance in the solution of nonlinear systems of equations. It is a global method that has the potential for finding all the roots to a nonlinear system, from which physically meaningful solutions can be extracted. The continuation technique avoids the disadvantage of Newton’s method, which can generate a singularity, or near-singularity of the Jacobian matrix and, thus, prevent convergence. A method that allows for a greatly expanded convergence domain relative to Newton’s method is differential-arc-length homotopy continuation. Wayburn and Seader (1987) discuss the conditions under which homotopy continuation is guaranteed to converge. When non-linearity results in multiple solutions, local convergence algorithms encounter singularities, e.g. at limit or turning points, where the iterates can oscillate between two steady-state attractors. In such situations, homotopy-continuation methods are very effective in achieving global convergence. Because of these features, the application of homotopy continuation to engineering is receiving much attention. Seydel and Hlavacek (1987) review the role of continuation methods in engineering problems. Another review on continuation methods and their application to separation problems is given by Wayburn and Seader (1987). Paloschi (1995) has developed a bounded homotopy, which can restrict the continuation path to a prescribed domain. Vickery and Taylor (1986) solved the MESH equations for a distillation tower by first assuming ideal liquid phases and then using a homotopy parameter to gradually introduce the nonideal terms into the model. Sophisticated homotopy continuation methods for computing azeotropes are given by Fidkowski, Malone and Doherty (1993), Eckert and Kubicek (1997). Kovach and Seider (1987) solve the liquid– liquid equilibrium (LLE) equations using a Newton homotopy. DeGance (1993) uses a continuation method with an equation of state to compute phase equilibria. He uses the pressure or a mole fraction as the homotopy parameter and calculates bubble and dew points. Sun and Seider (1995) use homotopy continuation for stability analysis by global minimization of Gibbs free energy for phase equilibrium with a check on stability of the system by a Newton homotopy-continuation method.

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331

2.1. Definition of homotopy function A system of nonlinear algebraic equations can be written as: f( (x¯ ) =0

(2)

where f: R “ R is a continuously differentiable function, and x R n is a variable vector. Homotopy functions provide a smooth transition along a path between an approximation to the solution and the true solution by gradually introducing a complex non-linearity. In many formulations a linear homotopy function is used, defined by: n

n

H (x¯, t)= tf (x¯ )+(1 −t) g (x¯ )

(3)

where f is the system of equations to be solved, g is a simple system of equations for which a solution is known, and 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. Three choices for g are widely used: l. Newton homotopy: H (x¯, t)= tf (x¯ )+(1 − t) [ f (x¯ ) − f(x¯ 0)]

(4)

2. Fixed-point homotopy: H (x¯, t)=tf (x¯ )+ (1− t) (x¯ −x¯ 0)

(5)

3. Affine homotopy: H(x¯, t)=tf (x¯ )+(1−t) f’(x¯ 0) (x¯ − x¯ 0)

(6)

where x 0 is the initial value of x and f % is the derivative of f.

Nonlinear optimization methods are often used in chemical process design. Applications extend from optimal design of individual pieces of process equipment to overall flowsheet optimization. Minimization of Gibbs free energy to determine equilibria in nonideal mixtures is another common application. General representation of an optimization problem is: Optimize subject to and

f(x¯, p¯ ) h( (x¯, p¯ ) = 0 g¯ (x¯, p¯ )]0

1. Inability to find more than just one local minimum. 2. Failure to converge from poor initial guesses due to the presence of nonlinear equalities. Floudas, Aggarwal and Ciric (1988) discuss the application of the generalized Benders’ decomposition to non-convex problems. They use decomposition as a heuristic tool during the search, but their method does not guarantee finding the global optimum. Much of the work on global optimization involves the unconstrained case, Torn and Zlinskas (1987). McDonald and Floudas (1994) present a technique for global optimization. They reformulate the MINLP and perform a global optimization algorithm to achieve o-global solution. In another work, McDonald and Floudas (1995) apply global optimization for the stability analysis of phase equilibria using NRTL and UNIQUAC activity coefficient models. The constrained optimization problem can be transformed into an unconstrained one by forming the Lagrange function, L, and applying the stationary conditions 9L = 0. The result is a system of nonlinear equations (NLE): 9 (f( + p¯ Th( + l( Tg¯ )=0

(7)

h( = 0

(8)

gili = 0

Öi

(9)

where p and l are Lagrange and Kuhn–Tucker multipliers, respectively, with the additional requirements that gi ] 0 li ] 0

3. Nonlinear optimization

1321

(10) (for minimum)

(11)

Singularities such as limit points and bifurcation points in equality constraints usually cause complex solution spaces. Singularities may involve the following problems: 1. Local optima may occur when the design variables are represented as parameters. 2. The optimization problem may become complicated with highly nonlinear inequality constraints.

4. Bifurcation analysis equality constraints inequality constraints

A number of methods for solving optimization problems in recent years have focused on the development of successive quadratic programming (SQP) strategies e.g. Locke, Edahl and Westerberg (1983), Biegler (1990)and Lucia and Xu (1990). These methods are capable of solving for local minima efficiently, but the starting point should be feasible. Two shortcomings are:

Because homotopy continuation paths often include singular points (such as limit points) where the determinant of the Jacobian matrix is zero, it is important to trace the paths with sufficient accuracy while bypassing the singular points. There are different types of algorithms available for solving this problem. One algorithm is presented by Keller (1977) and implemented in AUTO by Doedel (1986), in which a successive continuation technique is applied to find stationary points.

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331

1322

In homotopy-continuation algorithms, the steadystate solution branches are followed and only limit points are classified. Higher-order singular points are normally skipped. These include real bifurcation points where two or more solution branches cross, and Hopf bifurcation points where a periodic branch of solutions is born. Stability analysis can be performed to locate real and Hopf bifurcation points along the solution branches as shown by Kubicek and Marek (1983). At limit points, an eigenvalue becomes zero and a change in stability often occurs. Also, real bifurcation points have similar properties at the intersection of two or more branches.

4.1. Pseudo-arclength continuation The pseudo-arclength continuation method presented by Keller (1977) is given by:

!

(12)

=−



n n Du y1 Dp y1

( f 1p)y p; 0

f(u y1, p y1) (u −u0) u; 0 +(p y1 −p0) p; 0 −Ds y 1

T

n

(13)

The next direction vector (u; 1, p; 1) is defined by the equations f 1u u; 1 +f 1p p; 1 = 0

C

(15)

j=1 i=1

where P is the number of phases, C is the number of components, nij is the number of moles of component i in phase j, and G( ij is the partial molar Gibbs function or chemical potential expressed by: G( ij = G 0ij + RT ln aij

(16)

where G 0ij is the molar Gibbs function at the standard state and aij is the activity of the component i in phase j. The following constraints are present for phase equilibria. Conservation of mass: nif = % nij

(17)

j=1

This method finds a solution (u1, p1) to f in a hyperplane that is at a distance Ds from (u0, p0) and is perpendicular to the direction vector (u; 0, p; 0). Newton’s method for solving these equations is: ( f 1u)y u; T0

P

G= % % nij G( ij

P

f(u 1, p1) = 0 T (u 1 −u0) u; 0 +(p1 − p0) p; n0 −Ds =0



Gibbs function must reach its minimum value. For a multi-phase, multicomponent system, the Gibbs function is given by:

(14)

u; T0 u; 1 +p; 0 p; 1 =1 and (u; 1, p; 1) can be found by back-substitution at the end of the Newton iterations. Also, the orientation of the branch is preserved if Ds is sufficiently small.

5. Application to phase equilibrium Two methods for calculating the equilibrium composition of a mixture of chemical species are commonly applied: 1. Equilibrium constant method 2. Free energy minimization methods. There has been a continuing discussion of the advantages and disadvantages of these two methods since at least 1960. Today, however, in most computer programs minimization methods are the basis for phase equilibrium calculations. Using the criterion of equilibrium for a system under constraints of fixed pressure and temperature, the

where the subscript f refers to the total for all phases. Non-negativity: nij ] 0

(18)

It is assumed that each component 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 as shown by Eriksson and Rosen (1973). Smith (1980) shows that a number of equivalent formulations can be made, which can prove useful in specific situations. To show the efficiency of the continuation method for phase equilibrium, two test problems are presented here. First, a liquid–liquid system is considered. Second, a system of four phases (three liquid phases and one vapor phase) is presented.

5.1. Test problem 1 This test problem, which is taken from Heidemann and Mandhane (1973), involves an LLE system of one mole of water and one mole of n-butyl acetate at 25°C and 1 atm, using the NRTL activity coefficient model. For this system the Gibbs function is: 2

2

2

2

G= % % nijG( ij = % % nij [G 0ij + RT ln (gij xij )] i=1 j=1

(19)

i=1 j=1

The parameters for the NRTL model for estimating the activity coefficients, gij, are given by Heidemann and Mandhane (1973) as: t12 = 3.00498 t21 = 4.69071

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331

a12 =a21 =0.391965 Eq. (19) can be considered a constrained problem with Eqs. (17) and (18) as constraints. This constrained problem can be transformed to an unconstrained one by forming a Lagrangian function, L, as follows: 2

2

I=1

j=1

2

2

L = DG( + % pi % (nij −nif ) − % % nij lif

(20)

i=1 j=1

DG( =

In order to expand the range of feasibility of the constraint (Eq. (23)) the Mangasarian theorem (Vasudevan, Watson & Lutze, 1989) is used. This theorem replaces the complementary slackness condition (Eq. (23)) with nonlinear equations whose solution is feasible and has the proper sign for the KT multipliers. If a function u is defined so that u(0)= u%(0)=0 then: u{ gi (x)− li }− u{gi (x)}− u(li )= 0

where DG = % % bj xij ln(xij gij ) NTRT i = 1j = 1 m

2

2

1323

Öi

Function u can be either z z or z z 2. The results for mole fractions are:

2

DG m = G− % nif G 0i

Mij = (xij − lij )2 − xij xij − lij lij

i=1

where NT is the total number of the moles in the system, and bj is Nj /NT. Applying the stationary conditions 9L =0 to the above equation, results in a nonlinear system of equations (NLE) known as the Kuhn–Tucker (KT) necessary conditions (Reklaitis, Ravindran & Ragsdell, 1983): 9n¯L =9n¯ 9G( (n¯ )+ p¯ T 9n¯h( (n¯ ) − l( T9n¯S( (n¯ ) =0

(21)

where n is the mole number and 2

h( (n¯ ) = % nij − nif = 0

(22)

i =1,…C

(24)

or Mij = (1− xij − lij )2 − (1− xij ) 1− xij − lij lij

(25)

where for all of the above equations, i= 1, 2 and j=1, 2. To solve the above system of equations for LLE, Newton homotopy functions are used, where mole fractions replace mole numbers. tF(x¯, l( , p¯ , b( )+(1−t) (F(x¯, l( , p¯ , b( )− F(x¯ 0,l( 0,b( 0))=0 (26)

j=1

S( (n¯ )l( =n¯l( = 0

(23)

where

th( (x¯, b( )+ (1−t) [h( (x¯, b( )− h( (x¯0, b( 0)]= 0

nij ] 0 lij ] 0

tM( (x¯, l( , b( )+ (1− t) [M( (x¯, l( , b( )− M( (x¯ 0, l( 0, b( 0)]=0 (27) (28)

where F(x¯, l( , p¯ )=9x¯L

for minimum

5.2. Initialization Table 1 Mole fractions of water and n-butyl acetate in two phases of test problem 1 Set no.

1 2 3 4

Phase I

Phase II

x11

x21

x21

x22

0.00456436 0.00459335 0.00455711 0.05382077

0.9954356 0.9954067 0.9954429 0.9461792

0.9351421 0.8000118 0.5919872 0.5000000

0.06485794 0.1999882 0.4080128 0.5000000

Table 2 Gibbs free energy of the LLE system for test problem 1 Set no.

DG

1 2 3 4

−0.019606281 −0.017252349 −0.020198312 −0.017577543

Initial values for mole fractions are chosen so that they are consistent with the equations. It is better to choose either minimum or maximum possible values, e.g. x11 = 0.00001, x21 = 0.99999; choosing values in between sometimes leads to local solutions. The initial values of l are arbitrary positive values (like 0.01) and the initial values of p can be 1.0. In this problem, i=1 is for water and i= 2 is for n-butyl acetate. The results are shown in Tables 1 and 2. In Table 1, x11, x12 are the mole fractions of water in the first and second phases, respectively, while x21, x22 are the mole fractions of n-butyl acetate. Four solutions were found. As can be seen in Table 2, the Gibbs free energy of the system for the third solution set is lower than the other three. Hence, the global minimum of the system is − 0.020198312, and the other three solutions are local minima. Fig. 1 shows the homotopy path of the mole fraction of water in liquid phase 2 versus the homotopy parameter, t. The four solutions at t=1 are indicated by points.

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331

1324

Here, again the NRTL activity coefficient model is used with the following parameters: t12 = 4.9579 t13 = 10.9998 t21 = 0.01038 t23 = 1.4471 t31 = 6.2160 t321.9369 where the first component is water, the second component is aniline, and the third component is n-hexane. The value of a is taken equal to 0.2 for all three binary pairs.

5.4. Initialization

Fig. 1. Homotopy path for the mole fraction of water in liquid phase 2, test problem 1. Table 3 Mole numbers of water, aniline, and hexane for test problem 2, first solution set Component

Liquid phase I

Liquid phase II

Liquid phase III

Vapor phase

Water

0.9293322

0.0336938

3.348

0.03663921

−4

Aniline Hexane

0.00256858 1.268×10−6

0.4675863 0.0128917

× 10 0.0213838 0.9730883

0.5084613 0.01401868

Table 4 Mole numbers of water, aniline, and hexane for test problem 2, second solution set Component Liquid phase I Water Aniline Hexane

Liquid phase II

Liquid phase III

Vapor Phase −8

0.04324229 5.422×10 0.625801 0.3309566 0.00280074 0.00148118 0.9957177 4.125×10−7 1.55×10−6 8.164×10−7 0.9999972 4.76×10−7

The computation is started in the same way as in Test problem 1 but in terms of mole numbers instead of mole fractions. Even if the initial values for the mole numbers do not satisfy the constraints, e.g. n1j = 1.0, n2j = 0.00001, and n3j = 0.00001, the computation can proceed and find the solutions. The values of Lagrange multipliers are completely arbitrary. The results are summarized in Tables 3–5, where three sets of solutions are found along the homotopy path. The Gibbs functions for the first, second, and third sets of solutions are − 0.22793727, − 0.072608479, and − 0.072608477. Therefore, the global minimum for this system is the first solution set. The other two sets are local minima, almost identical to each other in the Gibbs function. Fig. 2 shows the homotopy branch for moles of aniline in the vapor phase versus the homotopy parameter. As can be seen in the figure, turning points are not smooth and in some places they are very sharp.

6. Application to simultaneous chemical and phase equilibria The subject of phase equilibria with chemical reaction is of great interest in chemical engineering. Many processes use operations where simultaneous chemical reactions and mass transfer occur in the same vessel. Examples are reactive distillation and chemical absorption.

5.3. Test problem 2 In this test problem, the equilibrium of a system of four phases (three liquid phases and one vapor phase) is investigated. The system involves three components, water, n-hexane, and aniline, with one mole of each component at the conditions of 25°C and 1 atm. The procedure for solving the system is similar to the first problem. However, since (L/(n is linear with respect to l, it is possible, as shown by Jalali (1992), to reduce the number of unknowns and therefore the number of equations without loss of generality.

Table 5 Mole numbers of water, aniline, and hexane test problem 2, third solution set Component

Liquid phase I

Liquid phase II

Liquid phase III

Vapor phase

Water

2.51×10−8

0.9567577

0.0399365

0.00330582

Aniline

2.373×10−7

0.0042819

0.9195964

0.07612144

Hexane

2.702×10−7

2.375×10−6

0.9235487

0.0764486

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331

1325

method and the SQP method is given by Lantagne, Marcos and Cayrol (1988). The Gibbs function also provides a means of identifying the equilibrium state when chemical reactions occur. The general procedure for solving multiphase chemical equilibria is to formulate the chemical equilibrium in one phase and relate the compositions of the phases by phase equilibrium relations. For vapor–liquid reactions, it is preferable to formulate the chemical reaction in the vapor phase and relate the compositions of the phases with the equilibrium flash equation. For liquid–liquid phases, it is more efficient to formulate the chemical equilibrium in the phase that has the greater concentrations of reactants, if possible. The total Gibbs function for the chemical and phase equilibria system is the same as Eq. (15) and the partial molar Gibbs function can be expressed as Eq. (16). The following constraints apply at CPE. 1. Conservation of mass: P

C

% % Aik nij = bk k= 1,…E Fig. 2. Homotopy path for the moles of aniline in vapour phase, test problem 2.

The calculation of simultaneous chemical and phase equilibria is an important aspect of process design. There are many examples in industrial applications that show the need for robust and efficient algorithms for CPE, such as phase equilibria in reactive distillation (Barbosa & Doherty, 1988), and the utilization of chemical theory (Prausnitz, Lichtenthaler & de Azevedo, 1986) to model phase equilibria, in which the components in a mixture react and form new chemical species. In CPE, the number and identity of the phases present at equilibrium are unknown in advance. Many researchers have followed the strategies used for phase equilibria such as assuming a large number of phases and eliminating the phases in excess during the course of the calculation (Castillo & Grossmann, 1981; Uchida, 1987), or starting from one phase and adding phase(s) one at a time during the calculation (Seider, Gautam & White, 1979; Michelsen, 1982). Smith, Missen and Smith (1993) propose a general optimality criteria for multi-phase, multi-reaction equilibria. Some CPE problems are linearly constrained, and therefore, direct minimization can also be performed (Harvie, Greenberg & Weare, 1987; Castier, Rasmussen & Fredenslund, 1989). The whole class of linear-constrained linear and nonlinear minimization methods is discussed by Gill, Murray and Wright (1981). The algorithm most often used in these cases is sequential quadratic programming (SQP), one implementation of which is in the OPT program by Biegler and Cuthrell (1985). A comparison study using a penalty function

(29)

j=1 i=1

2. Non-negativity: nij ] 0

(30)

where Aik is the coefficient of element k in component i, k is the subscript for the number of elements (E) 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. Assuming there are C components in the system, including reactants and products, and P phases, one can rewrite the objective function (Gibbs function) in terms of fugacity as: p

C

P

C



Minimise G= % % nijG( ij = % % nij G 0ij + RT ln J = 1i = 1

j = 1i = 1



f( ij f 0ij (31)

where the ratio of f( ij, the molar fugacity of component i in phase j, to f 0ij, the standard-state fugacity of component i phase j, replaces aij in Eq. (16). The constraints now are as follows: P

C

subject to h( (n¯ )= % % Aik nij − bk = 0

(32)

j = 1i = 1

and nij lij = 0 Thermodynamic models (e.g. fugacity forms) of the Gibbs function may not be convex functions. Thus, the global minimum may or may not be found during the optimization process. Non-uniqueness for nonideal systems has been discussed by many authors (e.g. Othmer, 1976; Caram & Scriven, 1976; Heidemann, 1978). Here, as in the previous cases, a Newton continuation technique for solving an optimization problem is used to

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331

1326

Table 6 Feed composition for test problem 3

Table 8 The values of Aik and bk for methanol synthesis

Components

Mole number

CO CO2 H2 H2O CH3OH CH4

15.0 8.0 74.0 0.0 0.0 3.0

Components

Aik C

CO CO2 H2 H2O CH3OH CH4

develop a procedure for predicting phase equilibria with chemical reaction and for finding the global minimum, even when the model is non-convex. To show the application of the homotopy-continuation method to CPE, two test problems are presented using the Soave–Redlich – Kwong equation of state. First, a system of vapor – liquid phases with reaction is considered; and second, a reacting system in three phases (vapor–liquid–liquid) is investigated.

bk

O

CH4

1 1 0 0 1 0

0 0 2 2 4 0

1 2 0 1 1 0

0 0 0 0 0 1

23

148

31

3

six components present at equilibrium and the number of the elements is four (C, O, H, CH4). Table 8 shows the values of Aik and bk for this problem. The Lagrangian of the system is as follows:



4

2



6

2

6

L= G+ % pk % % Aiknij − bk − % % nijlij k=1

6.1. Test problem 3

H

j = 1i = 1

(33)

j = 1i = 1

and the Kuhn–Tucker necessary conditions (9L =0) are

In this problem, the equilibrium calculations for the vapor–liquid synthesis of methanol are presented. Graaf, Sijtsema, Stamhuis and Joosten (1986) and Chang, Rousseau and Kilpatric (1986) studied this system using equations of state to model the non-ideality of the phases. In particular, Chang et al. (1986) reported that their method failed to converge at 473.15 K and 300 atm with the feed composition shown in Table 6. The calculation method presented here for this problem confirms that using the SRK equation with interaction parameters shown in Table 7 derived from Castier et al. (1989), predicts the equilibrium compositions in two phases at 473.15 K and 300 atm. The two main reactions in the methanol synthesis are:

4 (L = G 0ij + ln f( ij + % Aikpk − lij = 0 (nij k=1 2

(34)

6

h(nij )= % % Aiknij − bk = 0

(35)

j = 1i = 1

M(nij )= (nij − lij )2 − nij nij − lij lij = 0

(36)

where i= 1,…6, and j= 1, 2 for the above equations. Since Eq. (34) is linear with respect to lij, it is possible to reduce the number of variables by taking lij as a function of nij and pk in Eq. (34). Therefore, the number of equations reduces to 16 without loss of generality. Considering this, the system of equations is solved using a Newton homotopy.

CO +2H2 lCH3OH CO2 +H2 lCO +H2O

6.2. Initialization

Also, methane is added to the system as an inert, and since it is not part of the reacting system in the calculation, it is considered as an element. Therefore, there are

The computation was initiated with the following mole numbers of the components in the system, Lagrange multipliers, and Kuhn–Tucker multipliers.

Table 7 Interaction parameters (kij ) for the SRK equation of state in the methanol synthesis Component

CO

CO2

H2

H2O

CH3OH

CH4

C18H38

CO CO2 H2 H2O CH3OH CH4 C18H38

0.0 −0.032 0.0804 0.0603 −0.288 0.0322 0.0

0.0 −0.343 0.0737 0.0148 0.0933 0.0

0.0 −0.25 0.0 −0.022 0.0

0.0 −0.079 0.5 0.4

0.0 −0.25 0.0

0.0 0.0

0.0

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331

n11 =n12 =n41 =15,

n32 =74,

n22 =8,

n62 =3,

n31 =n21 =n42 =n52 = n61 =n51 =1.0 ×10 − 15, pi =10

Table 11 The values of Aik and bk for methanol synthesis for test problem 4 Components

i=1, . . . ,4

It was not necessary for the initial values of the mole numbers to be consistent with the constraints. However, it is better to choose them close to possible maximum or minimum values of the components in the system at equilibrium. Table 9 Methanol synthesis results for test problem 3 Components

Liquid phase (moles)

Vapor phase (moles)

CO CO2 H2 H2O CH3OH CH4

0.0002904 0.006381 2.70937 6.74251 17.6095 0.689805

0.001640 0.01366 17.3546 1.23744 5.36857 2.31020

1327

Aik C

CO CO2 H2 H2O CH3OH CH4 C18H38 bk

H

O

CH4

C18H38

1 1 0 0 1 0 0

0 0 2 2 4 0 0

1 2 0 1 1 0 0

0 0 0 0 0 1 0

0 0 0 0 0 0 1

23

208

61

3

10

The computation starts at t= 0 and the continuation branch continues until it passes through t= 1, which is a solution to the system of equations. The results are presented in Table 9. The value of the Gibbs function is equal to − 0.110291× 104. With different initial values only one solution was found. Fig. 3 shows the solution branch for moles of water in the liquid phase. The system of equations in this problem was scaled by some factors to get better changes near t= 0. In CPE problems, which are formulated in terms of moles, scaling is more important than in the phase equilibrium problem. In PE problems, the total number of moles is constant so the calculation can be done in terms of mole fractions that are in the range 0 and 1.

6.3. Test problem 4

Fig. 3. homotopy path for the moles of water on liquid phase, test problem 3. Table 10 Feed composition for test problem 4 Components

Mole number

CO CO2 H2 H2O CH3OH CH4 C18H38

15 8 74 30 0.0 3 10

This test problem is an extension of Test Problem 3 to an additional liquid phase at equilibrium for 473.15 K and 100 atm. Ko, Lee and Kulik (1987) studied vapor–liquid methanol synthesis in the presence of inert n-octadecane. In this test problem, the equilibrium calculation is performed in three phases using the SRK equation of state with interaction coefficients, kij, given by Castier et al. (1989). The feed composition for this problem is shown in Table 10. Also, in Table 11, the values of Aik and bk of the components in the system are tabulated. The same system of equations as presented in Test problem 3 is used here except for the following differences: number of phases= 3 number of components=7 number of elements=5 An additional liquid phase is introduced by addition of n-octadecane to the system. Here, n-octadecane is treated as an element. Therefore, the number of the elements is 5.

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331

1328

6.4. Initialization The initial values of the mole numbers of the components in the system and the Lagrange multipliers are as follows: n13 =15, n33 = 74, n23 =8, n41 =30, n63 =3, n73 =10 n31 =n21 =n42 =n52 = n61 =n51 =1.0 ×10 − 15, pi =10

i=1, . . . ,5

n11 =n12 =n22 =n32 = n62 =n43 =n53 =n71 =n72 =1.0× 10 − 15 The results for the three solutions that were found are summarized in Tables 12 – 14. The values of the Gibbs function for the first, second, and third sets of solutions are −0.183188×104, −0.182020 × 104, and − 0.1820197×104. Therefore, the global minimum to the problem is the first solution set. The other two sets, which have almost the same Gibbs function are local

minima. Fig. 4 shows the homotopy branch for the mole numbers of CO in the first liquid phase. Again in this problem, the homotopy functions were scaled. 7. Discussion For all of the problems, three kinds of homotopy functions were tested. The Newton homotopy was preferred over the others for the following reasons: 1. The equations being solved are highly nonlinear, and each equation does not depend mainly on any one variable. Therefore, the change in one variable cannot represent the major change in that equation. 2. The fixed-point homotopy can have a scaling problem, and in the above examples, scaling is an important factor in the computations. The Newton homotopy is self-scaling. 3. The affine homotopy had the problem of derivatives of functions. Since the functions are highly non-linear, their derivatives caused difficulties in the computations.

Table 12 Methanol synthesis results for test problem 4, first solution set Components

Liquid phase I (moles)

Liquid phase II (moles)

Vapor phase (moles)

CO CO2 H2 H2O CH3OH CH4 C18H38

0.00021528 0.074609 0.182091 24.2852 6.94718 0.011871 0.8736×10−13

0.0053354 0.594598 2.55286 2.45935 6.06425 0.423826 9.92959

0.049752 2.50445 26.8966 8.08183 6.75961 2.5643 0.070406

Table 13 Methanol synthesis results for test problem 4, second solution set Components

Liquid phase I (moles)

Liquid phase II (moles)

Vapor phase moles)

CO CO2 H2 H2O CH3OH CH4 C18H38

8.15337×10−5 0.0159119 0.388572 28.7958 8.18090 0.0254610 1.29875×10−13

0.00758335 0.368856 20.7068 8.78470 14.3707 2.96330 9.96223

2.88135×10−5 0.00139862 0.0785158 0.0333098 0.0544909 0.0112362 0.0377747

Table 14 Methanol synthesis results for test problem 4, third solution set Components

Liquid phase I (moles)

Liquid phase II (moles)

Vapor phase (moles)

CO CO2 H2 H2O CH3OH CH4 C18H38

8.1566×10−5 0.0159119 0.388573 28.7958 8.18090 0.0254610 2.708×l0−12

0.0074918 0.364405 20.4569 8.67868 14.1973 2.92754 9.84199

0.00012027 0.00585027 0.328421 0.139330 0.227928 0.0469996 0.158006

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331

1329

tion points and stability, and to maintain the branches in the real domain. Appendix A. Nomenclature

NT ni nij nif p po P R Rn s S T T u t x xi u; 0 x0 z

activity of component i in phase j coefficient of element k in component i total atoms of element k in the system number of components in the system number of elements function, objective function derivative of f with respect to u 0 derivative of f with respect to x 0 the partial molar fugacity of component i in phase j the standard-state fugacity of component i phase j Gibbs function (free energy) partial molar Gibbs function or chemical potential standard Gibbs function of component i in phase j function of x inequality constraint function homotopy function equality constraint function Lagrange function Mangasarian function of component i in phase j total number of moles in the system mole number of component i number of moles of component i in phase j number of moles of component i in the feed parameter of nonlinear function bifurcation parameter at starting point number of phases in the system gas constant real space arc-length parameter inequality function transpose temperature general variable homotopy parameter mole fraction, general variable mole fraction of component i direction of u at starting point initial value of x variable in Mangasarian function

Greek a DG Ds gij p

non-randomness parameter for NRTL model change in Gibbs free energy arc-length increment activity coefficient of component i in phase j Kuhn–Tucker multiplier

aij Aik bk C E f f 0u f 0x f( ij f 0ij G G( ij G 0ij Fig. 4. Homotopy path for the moles of CO in liquid phase 1, test problem 4.

In order to prevent a singularity during the computation, bifurcation analysis was used to bypass the singular points. Also, it sometimes happens that negative moles or mole fractions are produced during the iteration, causing the computation to stop. To prevent such a situation, the produced negative value was replaced by the previous value divided by ten. This procedure always allowed the computation to proceed. The prediction of a state of combined chemical and phase equilibrium, using the equations above, can be quite complicated because of the large number of nonlinear equations that must be solved simultaneously. Frequently, the equations involved can be simplified or reduced in number by recognizing that some species are only very slightly soluble in certain phases, and that some reactions may go virtually to completion or do not measurably proceed at all in some phases. However, even with such simplifications, the calculations are likely to be quite tedious. The non-ideality of a multiphase reacting system makes it possible to have multiple solutions, like the three sets of solutions in Test problem 4. Castier et al. (1989) solved this problem and achieved only one solution to the system. To find the global minimum to the Gibbs function as a criterion for equilibrium of a reacting system, it is necessary to consider all of the possible solutions to the system. Further details of the four test problems are given by Jalali (1992). AUTO software was used in this research to solve all of the test problems in order to determine the bifurca-

g(x) g H h L Mij

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331

1330

u t y y ( Other –

Lagrange multiplier Mangasarian function interaction parameter in NRTL model iteration number partial deferential sign (over bar) indicates a vector

References Allgower, E., & Georg, K. (1980). Simplicial and continuation methods for approximating fixed points and solutions to systems of equations. SIAM Re6iew, 22, 28. Archetti, F., & Schoen, F. (1984). A survey on the global minimization problem: general theory and computational approaches. Annals of Operations Research, 1, 87. Barbosa, D., & Doherty, M. F. (1988). The influence of equilibrium chemical reactions on vapor–liquid phase diagrams. Chemical Engineering Science, 43, 529. Biegler, L. T., & Cuthrell, J. E. (1985). Improved infeasible path optimization for sequential modular simulators. II: The optimization algorithm. Computers &Chemical Engineering, 9 (3), 257. Biegler, L. T. (1990). Chemical process simulation. Chemical Engineering Progress, 85 (10), 50. Caram, H. S., & Scriven, L. E. (1976). Non-unique reaction equilibria in nonideal systems. Chemical Engineering Science, 31, 163. Castier, M., Rasmussen, P., & Fredenslund, A. (1989). Calculation of simultaneous chemical and phase equilibria in nonideal systems. Chemical Engineering Science, 2, 37. Castillo, J., & Grossmann, I. E. (1981). Computation of phase and chemical equilibria. Computers & Chemical Engineering, 5, 99. Chang, T., Rousseau, R. W., & Kilpatric, R. (1986). Methanol synthesis reactions: calculations of equilibrium constant using equations of state. Industrial and Engineering Chemistry, Process Design and De6elopment, 25, 477. 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. DeGance, A. E. (1993). Ab initio equation of state phase equilibria computations via path continuation. Fluid Phase Equilibria, 89, 303. Doedel, E. (1986). AUTO 86 User Manual, Software for continuation and bifurcation problems in ordinary differential equations. Eckert, E., & Kubicek, M. (1997). Computing heterogeneous azeotropes in multicomponent mixtures. Computers & Chemical Engineering, 21, 347. 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. Fidkowski, Z. T., Malone, M. F., & Doherty, M. F. (1993). Computing azeotropes in multicomponent mixtures. Computers & Chemical Engineering, 17, 1141. Floudas, C. A., & Pardalos, P. M. (1990). A collection of test problems for constrained global optimization algorithms. Lecture notes in computer science, vol. 455. New York: Springer-Verlag. Floudas, C. A., Aggarwal, A., & Ciric, A. (1988). Global optimum search for nonconvex NLP and MINLP problems. Computers & Chemical Engineering, 13 (10), 1117. Garcia, C. B., & Zangwill, W. I. (1979). Determining all solutions

to certain systems of nonlinear equations. Mathematical Operations Research, 4, 1. Gill, P. E., Murray, W., & Wright, M. H. (1981). Practical Optimization. New York: Academic Press. Graaf, G. H., Sijtsema, J. J. M., Stamhuis, E. J., & Joosten, G. E. H. (1986). Chemical equilibria in methanol synthesis. Chemical Engineering Science, 41, 2883. Hansen, P., Jaumard, B., & Lu, S. H. (1989). Global optimization of univariant Lipschitz functions: II. New algorithms and computational comparisons. RUTCOR Research Report, Rutgers University, p. 23. Harvie, C. E., Greenberg, J. P., & Weare, J. H. (1987). A chemical equilibrium algorithm for highly nonideal multiphase systems: free energy minimization. Geochimica Cosmochimica Acta, 51, 1045. Heidemann, R. A. (1978). Non-uniqueness in phase and reaction equilibrium computations. Chemical Engineering Science, 33, 1517. Heidemann, R. A., & Mandhane, J. M. (1973). Some properties of the NRTL equation in liquid – liquid equilibrium data. Chemical Engineering Science, 28, 1213. Horst, R., & Tuy, H. (1990). Global optimization. Deterministic approaches. New York: Springer-Verlag. Jalali, F. (1992). Application of continuation methods to chemical and phase equilibria and stability analysis. Ph.D. Thesis, University of Utah. Keller, H. B. (1977). Numerical solution of bifurcation and nonlinear eigenvalue problems. In P. H. Rabinowitz, Applications of bifurcation theory. London: Academic Press. Ko, M. K., Lee, S., & Kulik, C. J. (1987). Multicomponent physical equilibrium of liquid-phase methanol synthesis. Energy & Fuels, 1 (2), 211. Kovach, J. W. III, & Seider, W. D. (1987). Heterogeneous azeotropic distillation: homotopy-continuation methods. Computers & Chemical Engineering, 11 (6), 593. Kubicek, M., & Marek, M. (1983). Computational methods in bifurcation theory and dissipati6e structures. New York: Springer-Verlag. Lantagne, G., Marcos, B., & Cayrol, B. (1988). Computation of complex equilibrium by nonlinear optimization. Computers & Chemical Engineering, 12 (6), 589. Li, T. Y. (1987). Solving polynomial systems. Mathematical Intelligence, 9 (3), 33. Locke, M. H., Edahl, R., & Westerberg, A. W. (1983). An improved successive quadratic programming optimization algorithm for engineering design problems. American Institute of Chemical Engineering Journal, 29, 5. Lucia, A., & Xu, J. (1990). Chemical process optimization using Newton-like methods. Computers & Chemical Engineering, 14 (2), 119. McDonald, C. M., & Floudas, C. A. (1995). Global optimization for the phase stability problem. American Institute of Chemical Engineering Journal, 41 (7), 1798 – 1814. McDonald, C. M., & Floudas, C. A. (1994). Decomposition based and branch and bound global optimization approach for the phase equilibrium problem. Journal of Global Optimization, 5, 205. Michelsen, M. L. (1982). The isothermal flash problem. Part II. Phase-split calculation. Fluid Phase Equilibria, 9, 21. Morgan, A. (1987). Sol6ing polynomial systems using continuation for engineering and scientific problems. NJ: Prentice-Hall. Othmer, H. G. (1976). Nonuniqueness of equilibria in closed reacting systems. Chemical Engineering Science, 31, 993. Paloschi, J. R. (1995). Bounded homotopies to solve systems of algebraic nonlinear equations. Computers & Chemical Engineering, 19 (12), 1243.

F. Jalali, J.D. Seader / Computers and Chemical Engineering 23 (1999) 1319–1331 Prausnitz, J. M., Lichtenthaler, R. N., & de Azevedo, E. G. (1986). Molecular thermodynamics of fluid-phase equilibria (2nd ed.). Prentice-Hall, Englewood Cliffs, NJ. Ratschek, H., & Rokne, J. (1988). New computer methods for global optimization. New York: Halsted Press. Reklaitis, G. V., Ravindran, A., & Ragsdell, K. M. (1983). Engineering optimization, methods and applications. New York: Wiley. 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. Seider, W. D., Gautam, R., White III, C. W. (1979). Computation of phase and chemical equilibrium: a review. Computer application to chemical engineering, ACS Symposium Series No. 124, 115. Seydel, J. D., & Hlavacek, V. (1987). Role of continuation in engineering analysis. Chemical Engineering Science, 42, 1281. Smith, W. R. (1980). The computation 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 equilibrium. American Institute of Chemical Engineering Journal, 39 (4), 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.

.

1331

Torn, A., & Zlinskas, A. (1987). Global optimization. Lecture notes in computer science No. 350. New York: Springer-Verlag. Uchida, M. (1987). MPEC2: a code for multi-phase chemical equilibria. Computers & Chemistry, 11 (1), 19. Vasudevan, G., Watson, L. T., & Lutze, F. H. (1989). A homotopy approach for solving constrained optimization problems. Proceedings of the 1989 American Cont. Conference (p. 780), Pittsburgh. Vickery, J., & Taylor, R. (1986). Path-following approaches to the solution of multicomponent, multistage separation process problems. American Institute of Chemical Engineering Journal, 32 (4), 547. 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 Trans. Mathematical Software, 13, 281. Wayburn, T. L., & Seader, J. D. (1987). Homotopy continuation methods for computer-aided process design. Computers & Chemical Engineering, 11 (1), 7. White, W. B., Johnson, S. M., & Dantzig, G. B. (1958). Chemical equilibrium in complex mixtures. Journal of Chemical Physics, 28, 751. Zitney, S. E., & Stadtherr, M. A. (1988). Computational experiments in equation-based chemical process flowsheeting. Computers & Chemical Engineering, 12 (12), 1171.