Computation of heteroazeotropes. Part II: efficient calculation of changes in phase equilibrium structure

Chemical Engineering Science 55 (2000) 3835}3853

Computation of heteroazeotropes. Part II: e$cient calculation of changes in phase equilibrium structure John E. Tolsma, Paul I. Barton* Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Received 15 March 1999; accepted 13 December 1999

Abstract This paper, the second in a two-part series, discusses how a new approach for the computation of the heteroazeotropes present in a multicomponent mixture can be extended to compute e$ciently the changes in phase equilibrium structure with parameters such as pressure. One particular advantage of the approach is the capability to detect incipient azeotropes and heteroazeotrope (i.e., azeotropes and heteroazeotropes that do not appear under current conditions but may appear under di!erent conditions). The algorithm for computing the azeotropes and heteroazeotropes is analyzed in the "rst paper in this series: Computation of Heteroazeotropes. Part I: Theory. This paper includes a brief discussion of the implementation of this approach and some numerical examples as well as examples illustrating the e!ectiveness in determining the bifurcation pressures at which azeotropes and heteroazeotropes appear and disappear or switch between each other.  2000 Elsevier Science Ltd. All rights reserved. Keywords: Heteroazeotropy; Continuation methods; Bifurcation analysis; Interval arithmetic

1. Introduction Computation of the azeotropes and heteroazeotropes present in a multicomponent mixture is a necessary task when analyzing and designing separation systems. In addition to computing the azeotropes and heteroazeotropes at a given pressure (or temperature), it is often valuable to know how the azeotropic composition and temperature vary with pressure (or, equivalently, how the azeotropic composition and pressure vary with temperature). A systematic approach for analyzing such changes can be incorporated directly into design algorithms, dramatically increasing the space of alternative designs. Given the composition, temperature, and pressure of an azeotrope or heteroazeotrope, standard continuation methods can be applied to the necessary conditions for

* Corresponding author. Tel.: #1-617-253-6526; fax: #1-617-2539695. E-mail address: [email protected] (P.I. Barton)

azeotropy to determine how the state variables change with a given parameter (e.g., temperature or pressure). Applying a phase stability test during the continuation, it is possible to determine at what conditions an azeotrope becomes a heteroazeotrope or vice versa. However, there often arises situations where at the speci"ed temperature or pressure, the azeotrope or heteroazeotrope does not even exist. Furthermore, performing a phase stability test at each point during the continuation is extremely costly. This paper describes a systematic approach for identifying incipient azeotropes and heteroazeotropes and the computation of the bifurcation pressures at which they appear and disappear or switch between each other. Section 2 brie#y describes the algorithm analyzed in Tolsma and Barton (2000a) and includes a brief description of the implementation. Section 3 contains a description of how the approach can be extended to analyze the phase equilibrium structure under parameter variation. Section 4 includes numerical examples for both the computation of the azeotropes and heteroazeotropes present in a multicomponent mixture as well as examples illustrating the ability to analyze the phase equilibrium structure.

0009-2509/00/$ - see front matter  2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 0 9 - 2 5 0 9 ( 0 0 ) 0 0 0 3 3 - 6

3836

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

2. Overview of homogeneous and heterogeneous azeotrope 5nding algorithm The "rst paper in this series, Computation of Heteroazeotropes. Part I: Theory (Tolsma & Barton, 2000a), describes a new approach for the computation of heteroazeotropes. The approach, which is an extension of Fidkowski's approach for computing homogeneous azeotropes (Fidkowski, Malone & Doherty, 1993), is based on solving the necessary conditions for azeotropy and heteroazeotropy using homotopy continuation. The homogeneous homotopy map is FM (x, ¹,j)




"



x (1![jK (¹, P, x)#(1!j)PQ (¹)/P])    $ x (1![jK (¹, P, x)#(1!j)PQ (¹)/P]) L L L L x !1 G G

"0, (1)

where x31L is the liquid composition, ¹ the temperature, P the pressure, and j31 is the homotopy parameter. At a "xed pressure, system (1) is n#1 equations in terms of n#2 variables and de"nes a 1-manifold (referred to as a homogeneous branch) in (x, ¹,j)-space. Let c (m)3FM \(0) denote this homogeneous branch where m is some suitable parameterization (e.g., arclength). At j"0, the homotopy map is equivalent to Raoult's Law and since Raoult's Law cannot predict azeotropy, there are precisely n solutions to this system of equations, the pure components. At j"1, the homotopy map reduces to the necessary conditions for azeotropy. As described in Tolsma and Barton (2000a), j is initialized to zero and x is initialized to each of the pure components and the homotopy branch de"ned by (1) is tracked to j"1. Along some of these pure component branches, bifurcation points appear which correspond to intersections with binary branches (branches with two nonzero mole fraction elements). These branches cross j"1 at solutions satisfying the necessary conditions for a binary azeotrope. Similarly, along some of the binary branches, bifurcation points are found that correspond to intersections with ternary branches which result in solutions satisfying the necessary conditions for a ternary azeotrope at j"1. In general, k-ary azeotropes are obtained through bifurcations on (k!1)-ary branches. When this approach is applied to systems exhibiting heteroazeotropy, solutions are found at j"1 which are spurious homogeneous solutions corresponding to actual heteroazeotropes. The

 Alternatively, temperature can be "xed, in which case (1) de"nes a 1-manifold in (x, P, j)-space.

heterogeneous homotopy map is FM M(x, x', x'', ¹, s,j)"




x ![jK' (¹, P, x')#(1!j)PQ (¹)/P]x'     $



x ![jK' (¹, P, x')#(1!j)PQ (¹)/P]x' L L L L j[c' (¹, P, x')x' !c'' (¹, P, x'')x'' ]#(1!j)[x' !x'' ]       $ j[c' (¹, P, x')x' !c''(¹, P, x'')x'']#(1!j)[x' !x''] L L L L L L x!sx'!(1!s)x'' L x !1 G G L x' !1 G G "0,

(2)

where x31L is the overall liquid composition, x', x''31L are the composition of liquid phases I and II, respectively, ¹ the temperature, P the pressure, s31 is the liquidphase fraction (the fraction of the total number of moles of liquid in liquid phase I), and j31 is the homotopy parameter. Similar to the homogeneous case, at a "xed pressure, Eq. (2) de"nes a 1-manifold (referred to as the heterogeneous branch and denoted by c M(m)) in (x, x', x'', ¹, s,j)-space. The heterogeneous branch has the following property: when the s-component on the branch crosses zero or unity, a projection of the heterogeneous branch intersects the corresponding spurious homogeneous branch. As the spurious homogeneous branches are tracked, the points of intersection are identi"ed. The heterogeneous map is then used to track the heterogeneous branch from the intersection point to the heteroazeotrope at j"1. Since only solutions satisfying the necessary conditions for azeotropy and heteroazeotropy are obtained using the approach described above, a phase stability test must be performed on all homogeneous and heterogeneous solutions. Another, independent, mechanism for obtaining the heteroazeotropes is also described in Tolsma and Barton (2000a). Here, heteroazeotropes which have corresponding spurious homogeneous azeotropes that lie outside the physical composition space (i.e., have compositions outside the physical range of zero and unity) can be computed by using a third homotopy map discussed in Tolsma and Barton (2000a). As shown in Tolsma and Barton (2000a), bifurcation and intersection points of interest when computing the azeotropes and heteroazeotropes present in a multicomponent mixture will occur within 0(j(1. However, if the homotopy branches are tracked outside this range, bifurcation and intersection points are sometimes identi"ed that correspond to nonphysical branches, that is, branches that either do not cross j"1 (the only physically meaningful point on the curve) or cross j"1 at

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

a point where mole fraction vector elements or the phase fraction (in the heterogeneous case) are outside their physical range of zero and unity. For example, if a pure component homogeneous branch is tracked backwards from j"0 (backwards de"ned as the direction of decreasing j) and a bifurcation point is identi"ed, then the corresponding binary branch is not physical because it is unable to cross j"0 in order to reach j"1. (Recall, as a consequence of Raoult's Law, the only solutions to the homotopy map at j"0 are the pure components and thus, branches of dimension two or greater cannot in general cross j"0.) If a bifurcation or intersection point is identi"ed at some j'1 then the corresponding branch either does not cross j"1 or does so at a point where the composition or phase fraction is outside the physical bounds of zero and unity. These nonphysical branches themselves are not of any use, however, as will be shown in Section 3, the nonphysical bifurcation and intersection points (i.e., points occurring outside 0(j(1) are useful in examining the phase equilibrium structure of the mixture. 2.1. Implementation There are four main numerical calculations performed in the azeotrope and heteroazeotrope "nding algorithm: branch tracking, bifurcation point identi"cation and branch switching, intersection point identi"cation and computation, and phase stability testing. These four calculations are described brie#y below. A more detailed discussion is presented in Tolsma (1999). 2.1.1. Branch tracking The major calculation performed in the azeotrope and heteroazeotrope "nding algorithm is the tracking of the homogeneous and heterogeneous branches. Since the branches may exhibit turning points (particularly, the heterogeneous branches), it is necessary to perform a continuation method capable of handling such features. All branch tracking in the numerical examples presented in this paper was performed using the locally parameterized continuation code PITCON (Rheinboldt & Burkardt, 1983). 2.1.2. Bifurcation point identixcation and branch switching As shown in Tolsma and Barton (2000a), while moving along a k-ary homogeneous branch, a necessary condition for a bifurcation point corresponding to an intersection with a (k#1)-ary branch is a (m)"0 for some j3+k#1,2, n, H

(3)

where a "1![jK #(1!j)PQ /P]. If the bifurcation H H H point is identi"ed by monitoring the sign of the determi-

3837

nant of the Jacobian at each step on the branch during the continuation there is a risk that an even number of bifurcation points might be missed due to a cancellation of the sign change of the determinant or we may incorrectly conclude there is one bifurcation point when there is an odd number greater than one. The occurrence of a bifurcation along the homotopy branch is analogous to a state event in a di!erential/algebraic equation (DAE) system, that is, a point at which the functional form of the DAE changes during the course of a simulation as a consequence of an event whose time of occurrence is a function of the state of the system (or implicit discontinuity). Typically, the state events are identi"ed by zero crossings of an appropriate discontinuity function. Park and Barton describe an algorithm for state event location that not only guarantees the identi"cation of the state event, but also correctly identi"es the "rst zero crossing of the discontinuity function (Park & Barton, 1996). Condition (3) is analogous to the discontinuity function of a hybrid discrete/continuous DAE model. By constructing an interpolating polynomial for each a , H j"k#1,2, n, using l previously computed points, the algorithm described in Park and Barton (1996) can be used to identify e$ciently the bifurcation points along the path. Once the bifurcation point has been identi"ed and computed, a point on the (k#1)-ary branch is obtained by solving the following system of equations:


 x (1!KM (x, ¹,j))   $

x (1!KM (x, ¹,j)) I I "0, I>x !1 G G a (x, ¹,j) H x !e H

(4)

where KM "jK #(1!j)PQ /P, k#1)j)n, and e is H H H some su$ciently small positive constant. Once this point is computed, it is saved so that this branch may be tracked later. The branch tracking is then continued on the original k-ary branch to j"1. In the case of a k-ary heterogeneous branch, the bifurcation criteria is d (m)"[s(m)!1]c' (m)![s(m)!KM ' (m)]c''(m)"0 H H H H for some j3+k#1,2, n,,

(5)

where c' "jc' #(1!j), c''"jc''#(1!j), and H H H H KM ' "jK' #(1!j)PQ . As with the homogeneous case, H H H polynomials can be "tted to the d 's and the state event H location algorithm mentioned above can be used to locate the bifurcation points. Once the points are identi"ed, the following system of equations is solved for a point on

3838

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

the new (k#1)-ary branch: x !KM ' x'    $ x !KM ' x' L L L c' x' !c'' x''     $ "0, c' x' !c''x'' L L L L x!sx'!(1!s)x''

(6)

1. there is no root in X, or 2. there is a unique root in X, or 3. there may or may not be one or more roots in X.

L x !1 G G L x' !1 G G d (x, x', x'', ¹, s,j) H x !e H where k#1)j)n and e is some su$ciently small positive constant. This point is saved so that that this branch may be tracked later. It should be noted that heterogeneous branches are primarily obtained through intersections with spurious homogeneous branches. 2.1.3. Intersection point identixcation and computation At the intersection of an k-ary spurious homogeneous branch and a projection of a heterogeneous branch the following system of equations is satis"ed:


 x (1!KM )   $

x (1!KM ) I I c x !c'' x''     "0. $

(7)

c x !c''x'' I I I I I x''!1 G G I x !1 G G Eq. (7) was formed from Eq. (2) by simply setting s"1 and x"x' and removing the redundant equations, including those associated with zero-valued variables x , j"k#1,2, n. Elements 1 through k and element H 2k#2 of Eq. (7) are satis"ed on the spurious homogeneous branch. The intersection "nding problem amounts to moving along the spurious homogeneous branch and looking for solutions to Eq. (7). Since the cost of this root-"nding problem increases with the number of variables, locating roots of the following subset of equations has been found to be computationally more e$cient and just as robust:




""

c x !c'' x''     $ c x !c''x'' I I I I I x''!1 G G

The intersection "nding problem is divided into two steps: (1) intersection point identi"cation and (2) intersection point computation. The second step is e!ectively a branch switching calculation. The identi"cation problem is handled using a root-exclusion test. A root-exclusion test performed on a system of equations f (z)"0, f : DL1LP1L, provides the following information about X-D:

"0.

For a multidimensional problem, a root-exclusion test based on interval arithmetic is ideal. Interval arithmetic is a branch of mathematics concerned with operations with intervals or closed subsets of the real line (Moore, 1977). Let X ,[XJ , XS ]"+x31 " XJ )x)XS , de     note an interval. The set of all intervals in 1 is denoted by (1. Similarly, the set of all n-dimensional intervals in 1L is denoted by (1L (e.g., the set of all rectangles in 1 is denoted by (1). The familiar operations associated with real numbers can also be de"ned for intervals. For example, let
denote some binary operation de"ned for real numbers. The associated interval operation is de"ned as X
> "+x
y " x3X , y3> ,. If x
y is unde    "ned for any x3X or y3> , the corresponding interval   operation is unde"ned. For example, interval division is unde"ned if the denominator contains zero. Similarly, elementary functions can be extended to intervals. Let func denote some function (e.g., sin, cos, etc.). The interval value for func, denoted by Func, is de"ned by Func(X )"+func(x) " x3X ,. Again, Func(X ) is unde   "ned if func(x) is unde"ned for any x3X . An interval  extension of a general nonlinear function f : DL1LP1L, denoted by F : (1LP(1L, can be constructed by replacing all real operators and elementary functions in f with the corresponding interval operations. The interval extension of f has the important property that given X3(1L (X-D), F(X) contains the image set of X under f, that is, f(X)-F(X). The cost of evaluating the interval extension of a function is a small multiple of the cost of evaluating the function using real operations. This provides a rapid way of answering question (1) above: If 0,F(X) then there is no solution to f (x)"0 for any x3X. Unfortunately, since the interval extension of a function over-approximates the image set, 03F(X) does not provide any information about solutions within X. Question (2) above can be answered by a theorem due to Moore (1977). First, the Krawczyk operator is de"ned: *(X),y!>f (y)#[I!> F(X)](X!y),

(8)

(9)

where y3X-D, >31L"L is an arbitrary nonsingular, real matrix, and F(X)3(1L"L is the interval extension of the Jacobian matrix of f.

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

Theorem 1 (Moore, 1979). If *(X)Lint(X),

(10)

then there exists a unique solution xH3X to f (x)"0. Furthermore, Newton's method will converge to this unique solution from any starting point within X. In practice, y"m(X) and >"(m( F(X)))\ where m( ) ) denotes the midpoint of an interval vector, i.e.,






(XS #XJ )/2   m(X)" $ . (XS #XJ )/2 L L The intersection "nding algorithm can be summarized as follows. At the jth step during the continuation on a spurious k-ary homogeneous branch, a search domain is constructed around the continuation point, (xH, jH). Let DH"+(x,j) " jH!djH)j)jH#djH and 0( x (1, i"1,2, k, denote this initial domain, where G djH is based on the size of the current continuation step. Since the trivial solution, x''"xH, always satis"es the intersection equations (system (8)), this point is removed by partitioning DH into 2I subdomains by bisecting through xH. These subdomains are placed into a list for further analysis. The root-exculsion test is applied to each of the subdomains in this list and if condition 1 occurs (no root exists) the subdomain is deleted from the list, if condition 2 occurs (unique root exist) the intersection point has been identi"ed and the process is stopped and this subdomain is returned, or if condition 3 occurs (unknown), the subdomain is bisected through its longest edge and these two new subdomains are inserted into the list. Except near the intersection point, condition 1 occurs often and the subdomains are rapidly excluded. A detailed description of the intersection point identi"cation algorithm is presented in Tolsma (1999). Once the intersection point has been identi"ed, the missing components associated with the heterogeneous branch lost in the projection must be computed in order for the heterogeneous branch to be tracked to the heteroazeotrope at j"1. The return value from the intersection "nding algorithm is either the empty set indicating no root exists near the current continuation point or a bounded region containing a unique root to the intersection equations. Furthermore, this nonempty region, if found, satis"es the condition in Moore's theorem and thus Newton's method can be guaranteed to compute the intersection point. In the examples tested in this paper, the branch tracking of the homogenous and heterogeneous branches was very rapid, taking less than 2 s on each branch. The process took longer when the intersection "nding algorithm was applied during the continuation (between 10 s and approximately 1 min for the examples tested), however, the heteroazeotropes were also obtained

3839

through bifurcations on lower-dimensional heterogeneous branches and through the use of the third homotopy map mentioned above, making the intersection "nding problem redundant except for the binary branches. 2.1.4. Phase stability test Since only necessary conditions for azeotropy and heteroazeotropy are satis"ed at j"1, all solutions must be checked for stability. A phase stability test based on the Gibbs tangent plane analysis has been applied in this work (Baker, Pierce & Luks, 1982; Michelsen, 1982). The global minimum of the tangent plane distance function, D(x), is identi"ed by computing all stationary points of D(x) using interval arithmetic (Stadtherr, Schnepper & Brennecke, 1995).

3. Analysis of phase equilibrium structure The condition of azeotropy occurs when there is an extremum in the equilibrium temperature and pressure surface of a mixture. At such a point, the composition of the equilibrium vapor and liquid are equal, and thus, this point acts as a barrier to separations exploiting composition gradients (or, more strictly, chemical potential gradients). Azeotropic points are univariate; specifying one intensive property uniquely de"nes the intensive state of the system. In some situations, changing a system parameter (e.g., pressure) results in the appearance or disappearance of an azeotrope. For example, Fig. 1 contains a schematic of a ¹}xy diagram for a binary mixture where the azeotropic state does not exist at low pressure, but emerges from a pure component vertex as pressure is increased. The pressure at which the azeotrope appears or disappears is referred to as a bifurcation pressure. This is illustrated for the acetone}ethanol system and the tetrahydrofuran}water system in Sections 4.2.2 and 4.2.6, respectively. Heteroazeotropy occurs when the azeotropic composition on the vapor}liquid equilibrium surface intersects the liquid}liquid binodal. The locus of azeotropic and heteroazeotropic points as a function of pressure are shown in the schematic in Fig. 2. Pressure does not strongly e!ect the liquid}liquid region, however, decreasing pressure reduces the boiling temperatures of the components present in the mixture, thereby causing the vapor}liquid butter#y to move into the liquid}liquid region. As a result, the azeotrope becomes a heteroazeotrope. The pressure where the azeotropic point on the vapor}liquid butter#y just touches the liquid}liquid binodal is also referred to as a bifurcation pressure. In many cases, the azeotropic and heteroazeotropic conditions are highly sensitive to changes in the system parameters. Knowing how the conditions of the system change is important because in some cases avoiding

3840

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

Fig. 1. Binary homogeneous ¹xy diagram.

3.1. Examination of the phase equilibrium structure There are three types of points of interest on the bifurcation diagrams generated using the approach described above: solution points (where the j component of the homotopy branch crosses one), bifurcation points (intersections of k-ary and (k#1)-ary branches), and intersection points (intersection of heterogeneous and spurious homogeneous branches). The number of variables in the homogeneous and heterogeneous homotopy maps is one greater than the number of equations and thus, there is one degree of freedom. The basic idea described in this section is to append an additional equation to each of the systems describing the points of interest on the bifurcation diagram, thereby making the system fully determined. This allows an additional parameter (e.g., pressure) to be treated as a variable so that standard continuation methods can be applied to the new system of equations to determine how the point varies. Fig. 2. Binary heterogeneous ¹}xy diagram.

a liquid}liquid phase split is desirable whereas in other cases, it may be advantageous to exploit this condition. This paper describes an e$cient procedure for analyzing the changes in phase equilibrium structure with a parameter such as pressure. The approach is an extension of an algorithm for the computation of the azeotropes and heteroazeotropes present in a multicomponent mixture described in Tolsma and Barton (2000a). Thus, it is possible to compute simultaneously the azeotropes and heteroazeotropes present under one set of conditions (e.g., speci"ed temperature or pressure) as well as predicting the e!ect on the azeotropic states of the system under parameter variation. The approach developed in this paper can also be used to examine the e!ect of physical property parameter variation on the prediction of phase equilibrium behavior. This application is described in Tolsma and Barton (2000b).

3.1.1. Solution points A solution point is simply a point on the homotopy branch that satis"es the necessary conditions for homoazeotropy or heteroazeotropy (i.e., a point at which the homotopy branch crosses j"1). In the homogeneous case, by treating y, x, ¹, and P as variables, the necessary conditions for azeotropy become 2n#1 equations in terms of 2n#2 variables. Similarly, in the heterogeneous case, treating y, x, x', x'', s, ¹, and P as variables, the necessary conditions for heteroazeotropy are 4n#2 equations in terms of 4n#3 variables. In both cases, we are e!ectively appending the equation j"1 to the homotopy maps. Standard continuation methods (Rheinboldt & Burkardt, 1983; Watson, Billups & Morgan, 1987) can then be applied to the underdetermined set of equations to examine how the azeotrope and heteroazeotrope conditions vary. Since only the necessary conditions are satis"ed at each continuation point, a phase stability test must be used during the continuation to determine when the azeotrope or heteroazeotrope

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

disappears or switches between each other. The need for a phase stability test makes the tracking of solution branches computationally unattractive. The examination of the phase equilibrium structure using the approach described above is straightforward and, with the exception of the phase stability test, quite e$cient. However, an azeotrope or heteroazeotrope must be known a priori in order to start the continuation, that is, we must start o! with a branch that crosses j"1. The following two sections describe how the existence of incipient azeotropes and heteroazeotropes can be predicted by examining the bifurcation and intersection points associated with nonphysical branches, branches that do not cross j"1 or do so at points that do not satisfy the necessary conditions due to a variable bound violation. 3.1.2. Bifurcation points A bifurcation point is where a k-ary branch intersects a (k#1)-ary branch. Without loss of generality, the following conditions are satis"ed on a k-ary homogeneous branch (Tolsma & Barton, 2000a): x O0, j"1,2, k, H a "0, j"1,2, k, H x "0, j"k#1,2, n, H and

(11) (12) (13)

a O0, j"k#1,2, n (14) H where a "1!(jK #(1!j)PQ /P)"1!KM . It is also H H H H shown in Tolsma & Barton (2000a) that a necessary condition for a bifurcation onto a (k#1)-ary branch at a point mI is a (mI )"0 for some j"k#1,2, n. (15) H A bifurcation from a k-ary branch onto a (k!1)-ary branch occurs when condition (11) above is violated for some component. Again, without loss of generality, assume for some mI , a (mI )"0. At a bifurcation point I> corresponding to the intersection of a k-ary and a (k#1)-ary branch the following set of k#2 equations are satis"ed: "@(x, ¹,j)






x (1![jK (¹, P, x)#(1!j)PQ (¹)/P])    $

" x (1![jK (¹, P, x)#(1!j)PQ (¹)/P]) I I I I x !1 G G a (¹, P, x,j) I> "0.

(16)

Treating pressure as a variable, Eq. (16) is k#2 variables in terms of k#3 variables (only the k nonzero mole

3841

fraction elements are considered). Again, standard continuation methods can be applied to this underdetermined system to examine how the bifurcation point varies. We are interested in under what conditions the bifurcation point occurs at j"0 or 1. A bifurcation point at j"1 corresponds to the case where the (k#1)ary azeotrope emerges from a k-ary azeotrope. This is illustrated in an example in Section 4.2.2 where an acetone}ethanol homogeneous azeotrope emerges from the pure acetone vertex as pressure is increased above approximately 8.6 bar. A bifurcation at j"0 will most commonly occur between a pure component and a binary homogeneous branch and it is a consequence of two components, which form an azeotrope, having the same boiling temperature at a certain pressure. This point corresponds to the case where the bifurcation point associated with the physical binary azeotrope switches the pure component branch from which it emerges from. This phenomenon is illustrated in Section 4.2.1 for the chloroform}methanol system. In the heterogeneous case, the necessary condition for a bifurcation from a k-ary branch onto a (k#1)-ary branch at a point mI is d (mI )"[s(mI )!1]c' (mI )![s(mI )!KM ' (mI )]c''(mI )"0 H H H H for some j"k#1,2, n

(17)

where c' "jc' #(1!j), c''"jc''#(1!j), and KM ' " H H H H H jK' #(1!j)PQ /P (Tolsma & Barton, 2000a). Again, H H without loss of generality, we can assume that condition (17) is satis"ed for j"k#1 at mI . Similar to the homogeneous case, the following set of equations are satis"ed at the heterogeneous bifurcation point: *M@(x, x', x'', s, ¹, P, j) x !KM ' (x', ¹, P,j)x'    $ x !KM ' (x', ¹, P,j)x' I I I c' (x', ¹, P,j)x' !c'' (x'', ¹, P,j)x''     $ c' (x', ¹, P,j)x' !c''(x'', ¹, P,j)x'' I I I I " x !sx' !(1!s)x'' "0. (18)    $ x !sx' !(1!s)x'' I I I I x !1 G G I x' !1 G G d (x, x', x'', ¹, P, s,j) I> Allowing pressure to vary, system (18) is 3k#3 equations in terms of 3k#4 variables (again, only the k nonzero mole fraction elements are treated as variables). As above, standard continuation methods can be applied to examine how the bifurcation point varies. Since a heteroazeotrope is typically formed when the azeotropic

3842

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

composition on the vapor}liquid butter#y moves into the liquid}liquid region, we are more interested in how the intersection points vary with parameters such as pressure. This is described in the following section. 3.1.3. Intersection points An intersection point is an intersection of a spurious homogeneous branch and a projection of a heterogeneous branch into (x, ¹,j)-space. Alternatively, the intersection point can be viewed as an intersection between a heterogeneous branch and a spurious homogeneous surface in (x, x', x'', s, ¹, j)-space. As shown in Tolsma & Barton (2000a), the rank of the Jacobian of the heterogeneous homotopy map, Eq. (2) in this paper, is de"cient by at least one. The case where the rank of this matrix is de"cient by exactly one corresponds to the case of a `normala intersection point while the others are degenerate and require further analysis. The set of constraints satis"ed at a k-ary intersection point are: "G(x, x'', ¹,j)"




x (1![jK (¹, P, x)#(1!j)PQ (¹)/P])    $



x (1![jK (¹, P, x)#(1!j)PQ (¹)/P]) I I I j[c (¹, P, x)x !c'' (¹, P, x'')x'' ]#(1!j)[x !x'' ]       $

j[c (¹, P, x)x !c''(¹, P, x'')x'']#(1!j)[x !x''] I I I I I I I x !1 G G I x''!1 G G "0. (19)

System (19) was formed from system (2) by removing x', s, the n!k zero mole fraction vector elements, and the equations de"ning the overall liquid composition, resulting in one less degree of freedom. The singularity of FM M at the intersection point has been removed by eliminating x!sx'!(1!s)x''"0 from (2). Allowing pressure to vary, standard continuation methods can be used determine how the intersection points vary with pressure.

4. Numerical examples This section contains numerical examples illustrating the approach described in Tolsma & Barton (2000a) for the computation of the azeotropes and heteroazeotropes present in a multicomponent mixture as well as examples illustrating the extensions described in this paper for examining changes in phase equilibrium structure. In the "rst subsection, the azeotropes and heteroazeotropes present in "ve ternary mixtures are computed. The second subsection contains seven examples illustrating

the approaches discussed in this paper. All numerical calculations were performed using the NRTL equation to model the liquid-phase activity coe$cients and the extended Antoine correlation was used to compute the vapor pressure. Constants for both models were obtained from Aspen Technology Inc. (1995). In all cases, the vapor phase was treated as an ideal gas. The NRTL model and extended Antoine correlation parameters for the systems examined can be found in Appendices 1 and 2. 4.1. Computation of azeotropes and heteroazeotropes The following ternary heterogeneous systems were examined: (1) benzene, ethanol, and water, (2) ethyl acetate, ethanol, and water, (3) water, acetone, and chloroform, (4) toluene, ethanol, and water, and (5) benzene, isopropanol, and water. Table 1 contains a summary of the experimental and computed values for the azeotropes and heteroazeotropes. Experimental values were obtained from Horsley (1973). The thick lines in Figs. 3}7 correspond to the heterogeneous branches and the dashed regions denote where the liquid-phase fraction is outside its physical bounds of zero and unity. The thin lines in these "gures correspond to homogeneous branches. All bifurcation diagrams in this section were constructed at a pressure of 1 bar. 4.1.1. Benzene}ethanol}water system At 1 bar, this system exhibits two binary homogeneous azeotropes, one binary heteroazeotrope, and a ternary heteroazeotrope. The binary homogeneous water} ethanol branch bifurcates o! the lower boiling ethanol branch (the binary azeotrope is minimum boiling). Similarly, the minimum boiling binary homogeneous benzene}ethanol branch also bifurcates o! the lower boiling ethanol branch. As shown in Tolsma & Barton (2000a) this must occur. Fig. 3 contains a bifurcation diagram showing the branches bifurcating o! the pure benzene branch. A spurious homogeneous benzene}water branch bifurcates o! the pure benzene branch. Along this branch, an intersection point is identi"ed from which the actual binary benzene}water heteroazeotrope is obtained using the heterogeneous homotopy map. Also identi"ed on the spurious binary branch is a bifurcation point onto a spurious homogeneous ternary branch. Moving forward in j, a spurious homogeneous ternary azeotrope is obtained. Moving in the reverse direction from the bifurcation point, an intersection with a heterogeneous ternary branch is identi"ed from which the actual ternary heteroazeotrope is obtained. In this example, the ternary heteroazeotrope can also be obtained through a bifurcation on the binary heterogeneous branch. There are two additional bifurcation points on the spurious ternary branch not shown in this "gure. The "rst occurs at the turning

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

3843

Table 1 Experimental and computed values for azeotropes and heteroazeotropes. Heteroazeotropes are denoted by boldface Experimental composition

Experimental temp. (K)

Computed composition

Computed temp. (K)

Benzene, ethanol, water

(0.5607,0.4393,0) (0.7003,0,0.2997) (0,0.9037,0.0963) (0.5387,0.2282,0.2331)

341.25 342.40 351.32 338.01

(0.5572,0.4428,0) (0.7021,0,0.2979) (0,0.8945,0.1055) (0.5293,0.2799,0.1908)

340.52 342.12 350.99 336.84

Ethyl acetate, ethanol, water

(0.5380,0.4620,0) (0.6885,0,0.3115) (0,0.9037,0.0963) (0.5789,0.1126,0.3085)

344.96 343.53 351.32 343.38

(0.5587,0.4413,0) (0.6755,0,0.3245) (0,0.8945,0.1055) (0.5415,0.1757,0.2828)

344.56 344.31 350.99 343.23

Water, acetone, chloroform

(0.1604,0,0.8396) (0,0.3397,0.6603) (0.1626,0.4844,0.3530)

329.25 337.15 333.55

(0.1669,0,0.8331) (0,0.3341,0.6659) (0.1855,0.3891,0.4255)

329.40 337.38 332.79

Toluene, ethanol, water

(0.1905,0.8095,0) (0.4439,0,0.5561) (0,0.9037,0.0963) (0.2736,0.3971,0.3293)

349.85 357.25 351.32 347.55

(0.1919,0.8081,0) (0.4414,0,0.5586) (0,0.8945,0.1055) (0.2607,0.4660,0.2732)

349.69 357.30 350.99 346.67

Benzene, isopropanol, water

(0.6022,0.3978,0) (0.7003,0,0.2997) (0,0.6875,0.3125) (0.5650,0.1861,0.2489)

344.89 342.40 353.25 339.66

(0.5991,0.4009,0) (0.7021,0,0.2979) (0,0.6356,0.3644) (0.5047,0.2597,0.2355)

344.49 342.12 347.08 338.14

Components

branches. However, the theory developed in Tolsma and Barton (2000a) indicates where the bifurcations can be found, limiting the amount of branch tracking required. Finally, since the spurious ternary azeotrope lies outside the composition space (i.e., the set C"+x31L " L x "1 G G and x *0 ∀i"1,2, n,), the ternary heteroazeotrope G can also be obtained by a third, independent mechanism by using the spurious azeotrope homotopy map described in Tolsma and Barton (2000a).

Fig. 3. Benzene mole fraction versus the homotopy parameter for ternary system: benzene, ethanol, and water. Thin lines are homogeneous curves and thick lines are heterogeneous curves.

point of the spurious ternary branch where there is an intersection with a homogeneous benzene}ethanol branch. In addition, the point at which the spurious ternary branch crosses x "0 corresponds to an intersection with a homogeneous ethanol}water branch. Higher-dimensional branches can typically be obtained through bifurcations on several lower-dimensional

4.1.2. Ethyl acetate}ethanol}water system At 1 bar, this system exhibits two homogeneous binary azeotropes, one binary heteroazeotrope, and a homogeneous ternary azeotrope. As in the system above, the homogeneous ethanol}water branch bifurcates o! the pure ethanol branch. Fig. 4 contains a bifurcation diagram showing the branches bifurcating o! the pure ethyl acetate branch. The "rst bifurcation point on the pure ethyl acetate branch corresponds to the ethyl acetate}ethanol branch from which the homogeneous binary azeotrope is obtained. A bifurcation point is identi"ed on this binary branch which corresponds to a ternary branch from which the homogeneous ternary azeotrope is obtained. A second bifurcation point is identi"ed along the pure ethyl acetate branch which corresponds to a spurious homogeneous ethyl acetate}water branch. Along this spurious branch an intersection point is identi"ed from which the actual ethyl acetate}water heteroazeotrope is obtained.

3844

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

Fig. 4. Ethyl acetate mole fraction versus the homotopy parameter for ternary system: ethyl acetate, ethanol, and water. Thin lines are homogeneous curves and thick lines are heterogeneous curves.

Fig. 5. Chloroform mole fraction versus the homotopy parameter for ternary system: water, acetone, and chloroform. Thin lines are homogeneous curves and thick lines are heterogeneous curves.

Section 4.2.4 also contains an example of this ternary mixture. There it is shown how a ternary heteroazeotrope is readily predicted at a di!erent pressure by moving slightly past j"1 on the ternary homogeneous branch.

azeotrope, like the ternary benzene}ethanol}water heteroazetrope, is also obtained through the spurious azeotrope homotopy map.

4.1.3. Water}acetone}chloroform system At 1 bar, this system exhibits one homogeneous binary azeotrope, one binary heteroazeotrope, and a ternary heteroazeotrope. Fig. 5 contains a bifurcation diagram showing the branches bifurcating o! the pure chloroform branch from which all solutions are obtained. The "rst bifurcation point identi"ed along the pure chloroform branch corresponds to the spurious chloroform}water branch. Along this branch an intersection point is identi"ed from which the actual binary heteroazeotrope is obtained. In addition, a bifurcation point is identi"ed along the spurious binary branch which corresponds to a spurious ternary branch. Three intersection points are identi"ed along the spurious ternary branch, all of which correspond to the same heterogeneous ternary branch. Two heterogeneous ternary solutions are identi"ed on this branch, one physical and the other nonphysical (the liquid phase fraction lies outside the range zero and unity). Note that in this case, the ternary heterogeneous branch cannot be obtained through a bifurcation on the binary heterogeneous branch. However, this is consistent with the theory in Tolsma and Barton (2000a) in that the ternary heterogeneous branch is obtainable through an intersection with a ternary spurious branch (under reasonable assumptions, all heteroazeotropes will be obtained either through intersections with spurious homogeneous branches or through bifurcations on lower dimensional heterogeneous branches). Furthermore, the ternary hetero-

4.1.4. Toluene}ethanol}water system At 1 bar, this system exhibits two homogeneous binary azeotropes, one binary heteroazeotrope, and one ternary heteroazeotrope. The binary ethanol}water azeotrope is obtained just as it was in the previous systems containing these two components, via a bifurcation on the pure ethanol branch. In addition, the homogeneous toluene}ethanol branch bifurcates o! the pure ethanol branch. Fig. 6 contains a bifurcation diagram showing the branches bifurcating o! the pure water branch. The spurious homogeneous toluene}water branch bifurcates o! the pure water branch very close to j"0. Along this branch an intersection with the binary heterogeneous branch and another bifurcation point are identi"ed. The actual binary heteroazeotrope and a spurious homogeneous ternary azeotrope are computed from these two points, respectively. A second intersection point is identi"ed along the spurious ternary branch from which the corresponding ternary heteroazeotrope is obtained. As with the benzene}ethanol}water system, the ternary heteroazeotrope can also be obtained through a bifurcation on the binary heterogeneous branch and through the spurious azeotrope homotopy map. Furthermore, the point at which the spurious ternary branch crosses x "0 corresponds to an intersection with the 5 homogeneous toluene}ethanol branch which bifurcates o! the pure ethanol branch and leads to the toluene}ethanol azeotrope at j"1. The dramatic di!erence between this bifurcation diagram and bifurcation diagram of the physically similar mixture, benzene,

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

3845

Fig. 6. Water mole fraction versus the homotopy parameter for ternary system: toluene, ethanol, and water. Thin lines are homogeneous curves and thick lines are heterogeneous curves.

Fig. 7. Benzene mole fraction versus the homotopy parameter for ternary system: benzene, isopropanol, and water. Thin lines are homogeneous curves and thick lines are heterogeneous curves.

ethanol, and water, is due to the fact that at a pressure of 1 bar, ¹Q (¹Q (¹Q , where ¹Q denotes the boiling 5 2 temperature.

4.2. Changes in phase equilibrium structure

4.1.5. Benzene}isopropanol}water system At 1 bar, this system exhibits two binary azeotropes, a binary heteroazeotrope, and a ternary heteroazeotrope. Fig. 7 contains a bifurcation diagram showing the branches bifurcating o! the pure benzene branch. The "rst bifurcation point on the pure benzene branch corresponds to the spurious benzene}water branch, along which an intersection point is identi"ed that leads to the actual binary heteroazeotrope. The second bifurcation point on the pure benzene branch corresponds to a benzene}isopropanol branch, leading to an actual binary azeotrope. Along the binary heterogeneous branch, a bifurcation point is identi"ed which corresponds to a ternary heterogeneous branch, from which the actual ternary heteroazeotrope is obtained. Note, however, in this case the spurious ternary branch is not obtainable through bifurcations on lower dimensional homogeneous branches. The spurious ternary branch is identi"ed (drawn as a medium thickness line in Fig. 7) by monitoring the s-component on the ternary heterogeneous branch and this spurious ternary branch lies completely outside the physical composition space, C. As with the water}acetone}chloroform system above, this is consistent with the theory in Tolsma and Barton (2000a) in that the ternary heterogeneous branch is obtainable through a bifurcation from a lowerdimensional heterogeneous branch. Finally, as with the other examples in this paper, the ternary heteroazeotrope is also obtained with the spurious azeotrope homotopy map.

In this section, examples associated with the extensions of the heteroazeotrope "nding algorithm discussed in this paper are presented. There is limited experimental data on azeotropes and heterazeotropes at the pressures computed in this section. However, all solutions obtained satisfy the necessary and su$cient conditions for homoazeotropy and heteroazeotropy. 4.2.1. Chloroform}methanol system A mixture of chloroform and methanol will form a homogeneous azeotrope at a pressure of 1 bar. The homogeneous homotopy branch associated with this minimum boiling binary azeotrope bifurcates o! the pure chloroform branch at j"0.028. However, a chloroform}methanol bifurcation point also appears on the pure methanol branch at a value of j"!0.073. This branch does not lead to a binary azeotrope since it is unable to cross j"0 in order to reach j"1. As pressure is increased, the value of j at the bifurcation point decreases on the pure chloroform branch and increases on the pure methanol branch. At a pressure of 1.8 bar, the bifurcation points on both branches occur at j"0. At higher pressures, the bifurcation point on the chloroform branch occurs at j(0 and the bifurcation point on the methanol branch occurs at j'0. Thus, the branch associated with the binary azeotrope now bifurcates o! the pure methanol branch. Fig. 8 contains the bifurcation diagram for this system at three di!erent pressures. This behavior can be explained physically by looking at the boiling temperatures of pure chloroform and methanol as a function of pressure. As shown in Tolsma and Barton (2000a), the minimum boiling

3846

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

Fig. 8. Bifurcation diagram for the chloroform-methanol system at pressures of 1.5, 1.8 and 2.0 bar.

Fig. 9. ¹}xy diagram for the chloroform}methanol system at 1.5, 1.789, 2.0 bar.

chloroform-methanol azeotrope will be obtained through a bifurcation on the lower boiling component's branch. At pressures less than 1.8 bar, chloroform has a lower boiling point than methanol. At 1.8 bar, both chloroform and methanol have computed boiling points of 352.8 K. At pressures greater than 1.8 bar, methanol boils at a lower temperature than chloroform. Fig. 9 contains the T}xy diagram for the chloroform}methanol system illustrating this behavior. 4.2.2. Acetone}ethanol system At a pressure of 1 bar, a binary acetone}ethanol azeotrope is not predicted. However, if the pure acetone branch is extended beyond j"1, a bifurcation point corresponding to an acetone}ethanol branch is identi"ed at j"1.8. This value of j at the bifurcation point decreases with increasing pressure and equals unity at a pressure of 8.6 bar. At this point, the acetone}ethanol azeotrope emerges from the pure acetone vertex. Experimental data indicate this azeotrope appears around 6.7 bar (Gmehling, Menke, Krafezyk & Fischer, 1994a,b). Although the NRTL model and associated parameters do not quantitatively predict this bifurcation point to a su$cient degree of accuracy, the qualitative behavior is correct and quite useful for gaining insight into this system. 4.2.3. Water}2-butoxyethanol system A mixture of water and 2-butoxyethanol is predicted to form a homogeneous azeotrope at a pressure of 1 bar. An intersection with a heterogeneous branch is identi"ed along the binary homogeneous branch at j"1.03. Following this heterogeneous branch back to j"1, a solution is found that satis"es the necessary conditions for heteroazeotropy except for the fact that the liquid-phase fraction equals 1.01. Fig. 10 contains a diagram showing

Fig. 10. Pressure versus the value of j at the intersection point for the water}2-butoxyethanol system.

how the intersection point varies with pressure. The intersection point crosses j"1 at a pressure of 1.18 bar. At this pressure, the azeotropic point on the vapor}liquid butter#y just touches the liquid}liquid region. The water}2-butoxyethanol mixture exhibits both an upper and a lower critical solution temperature. However, the point at which the heteroazeotrope leaves through the top of the liquid}liquid region was not predicted with this model. A better description of the vapor phase nonideality may rectify this de"ciency. The analysis above indicates that near a pressure of 1 bar, the homogeneous water and 2-butoxyethanol azeotrope lies very close to the liquid}liquid region (and moves into this region at a slightly higher pressure). The experimental data indicate that at 1 bar, water and

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

3847

2-butoxyethanol form a heteroazeotrope which lies just inside the liquid}liquid region (Gmehling et al., 1994a,b). This heteroazeotrope becomes a homogeneous azeotrope below approximately 0.25 bar and above 1 bar (the bifurcation pressures where the heteroazeotropes become homogeneous azeotropes and vice versa were not found in the literature). This example illustrates the problems encountered when estimating physical property parameters particularly near complicated regions (e.g., near the liquid}liquid boundary); very small errors can lead to qualitatively di!erent behavior (e.g., heteroazeotropy versus homogeneous azeotropy). 4.2.4. Ethyl acetate}ethanol}water system As shown in Section 4.1.2, at a pressure of one bar, a mixture of ethyl acetate, ethanol, and water exhibits a homogeneous ethanol}water azeotrope, a heterogeneous ethyl acetate}water azeotrope, and a homogeneous ternary azeotrope. In the bifurcation diagram shown in Fig. 4, the continuation was stopped at j"1. Fig. 11 contains the same bifurcation diagram except with the continuation performed beyond j"1 on the ternary homogeneous branch. If an intersection search is performed beyond j"1, an intersection with a ternary heterogeneous branch is located at j"1.013. The corresponding heterogeneous ternary branch crosses j"1 with a value of s equal to 1.008. Fig. 12 contains a diagram showing how this intersection point varies with pressure. At a pressure of approximately 0.67 bar, the intersection point occurs at j"1 corresponding to the pressure where the azeotropic composition on the liquid}vapor butter#y just touches the liquid}liquid binodal (which exhibits an upper critical solution temperature for this system). The equilibrium temperature at this point is approximately 593C. Experimental data found in Gmehling et al. (1994a,b) indicate that a ternary heteroazeotrope is found at a temperature less than 403C. Again, the qualitative behavior of the predicted solutions provide useful insights into the actual physical behavior of this system. 4.2.5. Water}1,2-dichloroethane system A mixture of water and 1,2-dichloroethane forms a heteroazeotrope at 1 bar. The intersection point varies as a function of pressure as shown in Fig. 13. At a pressure of approximately 0.057 bar, the heteroazeotrope leaves the liquid}liquid region and the heteroazeotrope becomes a homogeneous azeotrope. Unfortunately, the lowest pressure data found in the literature was 0.17 bar, where a heteroazeotrope is reported (Gmehling et al., 1994a,b). 4.2.6. Tetrahydrofuran}water system The NRTL model parameters found in Aspen Technology Inc. (1995) predict the behavior of the tetrahyd-

Fig. 11. Ethyl acetate mole fraction versus the homotopy parameter at 1.0 bar for ternary system: ethyl acetate, ethanol, and water. Thin lines are homogeneous curves and thick lines are heterogeneous curves. The continuation is performed beyond j"1 in this "gure.

Fig. 12. Pressure versus the value of j at the intersection point for the ethyl acetate, ethanol, water system.

rofuran and water mixture quite well over a range of pressures compared to the experimental data summarized in Gmehling et al. (1994a,b). At a pressure of 1.01 bar, the predicted binary homogeneous azeotrope is x "0.8287 and ¹"63.43C. Several di!erent values 2&$ are reported at this pressure in Gmehling et al. (1994a,b) with 0.80)x )0.83 and 63.253C)¹)703C. 2&$ At a pressure of 3.4996 bar, the computed value of the homogeneous azeotrope is x "0.7320 and ¹" 2&$ 103.923C and the experimental value is x "0.6975 2&$ and ¹"101.753C. Finally, at a pressure of 0.2181 bar,

3848

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

Fig. 13. Pressure versus the value of j at the intersection point for the water}1,2-dichloroethane system.

the computed value of the homogeneous azeotrope is x "0.9194 and ¹"24.783C and the experimental 2&$ value is x "0.9410 and ¹"253C. Using the 2&$ approach described in this paper, it was found that the homogeneous azeotrope will disappear (i.e., leave 0)x )1) at a pressure of 0.0315 bar, however, 2&$ experimental data was not found that supports this prediction. An intersection point search was performed along the 1.01 bar homogeneous binary branch and an intersection with a heterogeneous branch was identi"ed at j"1.32. Fig. 14 shows how the position of the intersection point varies with pressure. This intersection point does not cross j"1 (where the homogeneous azeotrope enters the liquid}liquid binodal) at any pressure, however, at high pressures the j-position of the intersection point approaches unity. This may be due to the fact that a liquid}liquid binodal, exhibiting both an upper and a lower critical solution temperature, exists for the THF}water mixture (Sorensen & Arlt, 1979), but the homogeneous azeotropic point does not move into this region. A much richer, yet incorrect, phase behavior of the THF}water mixture was predicted using a di!erent set of NRTL model parameters (found in an earlier version of Properties Plus). Using these parameters, the model predicted that a mixture of THF and water does not form an azeotrope or heteroazeotrope at pressures less than approximately 0.64 bar, forms a homogeneous azeotrope at pressures between 0.64 and 2.3 bar, forms a heteroazeotrope between 2.3 and 18.5 bar, and forms a homogeneous azeotrope above 18.5 bar. Fig. 15 shows how the location of the intersection point on the spurious homogeneous binary branch varies with pressure. Fig. 16 shows the

Fig. 14. Pressure versus the value of j at the intersection point for the THF}water system.

Fig. 15. Pressure versus the value of j at the intersection point for the THF}water system.

same diagram with the region around j"1 (the region of interest) enlarged. This system is di!erent from the others in this paper in that there are two disconnected intersection point curves. The curve associated with the intersection points found at pressures less than approximately 2.3 bar lie on a curve that does not cross j"1 except at a non-physical solution at approximately 89 bar. This curve does not provide much useful information except for suggesting a second curve exists. However, it may be possible that the two intersection point curves are connected via a complex curve that bifurcates o! the low pressure curve at the turning point at j+1.07. This

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

Fig. 16. Pressure versus the value of j at the intersection point for the THF}water system. Enlarged region of interest at j"1.

connection in the complex domain was not explored, however, since the second curve is easily found by performing the intersection point search at a single pressure above 2.3 bar. The NRTL parameters used in this second set of calculations did a very good job at predicting the azeotropic conditions around 1 bar, but failed at other pressures. The techniques developed in this paper were used to examine rapidly the phase behavior predicted with a set of parameters and, upon comparison with experimental data, the analysis suggested other parameters should be found. 4.2.7. Water}acetone}chloroform system As shown in Section 4.1.3, at a pressure of 1 bar, a mixture of water, acetone, and chloroform forms a water}chloroform heteroazeotrope, a homogeneous acetone}chloroform azeotrope, and a ternary heteroazeotrope. Fig. 5 contains the bifurcation diagram for this system constructed at 1 bar. The single intersection point on the spurious water}chloroform homogeneous branch does not cross j"1 for any pressure. There are three intersection points on the ternary homogeneous branch. Fig. 17 shows how these intersection points vary with pressure. At 1 bar, the intersection points on the ternary branch at j"0.74 and 0.95 actually lie on the same intersection point curve. At pressures above approximately 0.65 bar, three intersection points appear on the ternary branch. At a pressure of 0.65 bar, two of the intersection points disappear and as pressure is decreased, the intersection point (originally at j"0.67 at 1 bar) crosses j"1 at approximately 0.32 bar. At this point, the heteroazeotrope becomes a homogeneous azeotrope. Unfortunately,

3849

Fig. 17. Pressure versus the value of j at the intersection point for the water, acetone, chlorofrom system.

no data were found at pressures other than 1 atm to support these predictions. 5. Conclusions This paper describes how the approach described in Tolsma & Barton (2000a) can be extended to compute e$ciently the changes in phase equilibrium structure with parameters such as pressure. The approach has the advantage of identifying incipient azeotropes and heteroazeotropes, that is, azeotropes and heteroazeotropes that do not exist under given conditions, but may exist at other conditions. Finally, the approach described in Tolsma & Barton (2000a) was used to successfully compute all azeotropes and heteroazeotropes present in several multicomponent mixtures. In addition, the examples illustrate three mechanisms by which the heteroazeotropes can be obtained: intersections with spurious homogeneous branches obtained through bifurcations on lower dimensional homogeneous branches, intersections with spurious homogeneous branches isolated from the other homogeneous branches but obtained through a spurious homogeneous homotopy map described in Tolsma & Barton (2000a), and bifurcations on lower dimensional heterogeneous branches. Acknowledgements This research was supported by the United States Department of Energy under grant number DE-FGO294ER14447. We are grateful to Aspen Technology Inc for granting permission to reproduce the parameters in the Appendices.

3850

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

Appendix A. Physical property constants for the NRTL model The NRTL model used in the calculations has the following form:






xq G xG x q G H GH q ! I I IH IH , ln c " H H HG HG # G GH xG x G x G H H HG I I IH H I I IH G "exp(!a q ), GH GH GH q "A #B /¹#E ln¹#F ¹, GH GH GH GH GH a "C #D (¹!273.15), GH GH GH q "0, GG G "1, GG where ¹ is in Kelvins. The binary parameter arrays A, B, C, D, E, and F were obtained from Aspen Technology, Inc.'s Properties Plus (1995). The values used for the calculations presented in this paper are shown in Tables 2}11.

Table 3 Physical property constants for the NRTL model from properties plus database. Ethyl acetate, ethanol, and water system Constants

Ethyl Acetate

A

Ethyl acetate Ethanol Water

B

Ethyl acetate * Ethanol 524.42 Water !1705.7

C

Ethyl acetate Ethanol Water

D

Ethyl acetate Ethanol Water

E

F

* !1.1512 9.4632

0.30 0.30 0.20

Ethanol !0.24310 * 3.4578 282.96 * !586.0809

Water !3.7198 !0.8009 * 1286.1 246.18 *

0.30 0.30 0.30

0.20 0.30 0.30

* * *

* * *

* * *

Ethyl acetate Ethanol Water

* * *

* * *

* * *

Ethyl acetate Ethanol Water

* * *

* * *

* * *

Table 2 Physical property constants for the NRTL model from properties plus database. Benzene, ethanol, and water system Constants

Benzene

A

Benzene Ethanol Water

* 0.5686 140.0874

B

Benzene Ethanol Water

* !54.8044 !5954.307

C

Benzene Ethanol Water

0.30 0.30 0.20

D

Benzene Ethanol Water

E

Benzene Ethanol Water

F

Benzene Ethanol Water

* * * * * !20.0254 * * *

Ethanol

Water

Table 4 Physical property constants for the NRTL model from properties plus database. Water, acetone, and chloroform system

!0.9155 * 3.4578

45.1905 !0.8009 *

882.0288 * !586.0809

591.3676 246.18 *

A

Water Acetone Chloroform

0.20 0.30 0.30

B

Water Acetone Chloroform

C

Water Acetone Chloroform

D

Water Acetone Chloroform

* * *

* * *

* * *

E

Water Acetone Chloroform

* * *

* * *

* * *

F

Water Acetone Chloroform

* * *

* * *

* * *

0.30 0.30 0.30 * * * * * * * * *

* * * !7.5629 * * * * *

Constants

Water

Acetone

* 6.3981 !7.3519 * !1808.991 3240.688 0.30 0.30 0.20

0.0544 * 0.5382

Chloroform 8.8436 0.9646 *

419.9716 !1140.115 * !590.026 !106.4216 * 0.30 0.30 0.30

0.20 0.30 0.30

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

3851

Table 5 Physical property constants for the NRTL model from properties plus database. Toluene, ethanol, and water system

Table 7 Physical property constants for the NRTL model from properties plus database. Chloroform and methanol system

Constants

Constants

Toluene

A

Toluene Ethanol Water

B

Toluene Ethanol Water

C

Toluene Ethanol Water

D

Toluene Ethanol Water

E

Toluene Ethanol Water

F

Toluene Ethanol Water

* 1.1459 627.0528

Ethanol

Water

!1.7221 * 3.4578

* 992.7367 !113.4658 * !27269.36 !586.0809 0.30 0.30 0.20 * * * * * !92.7182 * * *

0.30 0.30 0.30

!247.8792 !0.8009 * 14759.76 246.18 * 0.20 0.30 0.30

* * *

* * *

* * *

35.5820 * *

* * *

* * *

Table 6 Physical property constants for the NRTL model from properties plus database. Benzene, isopropanol, and water system Constants

A

Benzene Benzene Isopropanol Water

* !0.1279 140.087

Isopropanol !0.74840 * 6.82840

B

Benzene * Isopropanol 148.121 Water !5954.31

C

Benzene Isopropanol Water

D

Benzene Isopropanol Water

* * *

* * *

* * *

E

Benzene Isopropanol Water

* * !20.0254

* * *

!7.5629 * *

F

Benzene Isopropanol Water

* * *

* * *

* * *

0.30 0.30 0.20

713.312 * !1483.46 0.30 0.30 0.30

14759.76 426.398 * 0.20 0.30 0.30

Methanol

A

Chloroform Methanol

* *

B

Chloroform Methanol

* !71.903

C

Chloroform Methanol

D

Chloroform Methanol

* *

* *

E

Chloroform Methanol

* *

* *

F

Chloroform Methanol

* *

* *

0.30 0.30

* * 690.07 * 0.30 0.30

Table 8 Physical property constants for the NRTL model from properties plus database. Acetone and ethanol system Constants

Acetone

Ethanol

* !1.0787

!0.34712 *

A

Acetone Ethanol

B

Acetone Ethanol

* 479.05

C

Acetone Ethanol

0.30 0.30

0.30 0.30

D

Acetone Ethanol

* *

* *

E

Acetone Ethanol

* *

* *

F

Acetone Ethanol

* *

* *

Water 45.1905 !1.31150 *

Chloroform

206.6 *

Table 9 Physical property constants for the NRTL model from Properties Plus database. Water and 2-butoxyethanol system Constants

Water

2-Butoxyethanol

A

Water 2-Butoxyethanol

* !0.6956

6.0952 *

B

Water 2-Butoxyethanol

* 138.501

C

Water 2-Butoxyethanol

0.30 0.30

0.30 0.30

D

Water 2-Butoxyethanol

* *

* *

E

Water 2-Butoxyethanol

* *

* *

F

Water 2-Butoxyethanol

* *

* *

!940.369 *

3852

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

Table 10 Physical property constants for the NRTL model from Properties Plus database. Water and 1,2-dichloroethane system

Table 11 Physical property constants for the NRTL model from Properties Plus database. Tetrahydrofuran (THF) and water system

Constants

Constants

Water

1,2-Dichloroethane

THF

Water

A

Water 1,2-Dichloroethane

* *

* *

A

THF Water

1.21416 *

B

Water 1,2-Dichloroethane

* 367.201

1848.53 *

B

THF Water

C

Water 1,2-Dichloroethane

0.30 0.30

0.30 0.30

C

THF Water

D

Water 1,2-Dichloroethane

* *

* *

D

THF Water

* *

* *

E

Water 1,2-Dichloroethane

* *

* *

E

THF Water

* *

* *

F

Water 1,2-Dichloroethane

* *

* *

F

THF Water

* *

* *

* 4.76015 * !73.3402

157.781 *

0.30 0.3472553

0.472553 0.30

Table 12 Physical property constants for the extended Antoine correlation Constants

Water

Methanol

Ethanol

Isopropanol

C  C  C  C  C  C  C 

73.6490 !7258.20 0.0 0.0 !7.30370 4.16530E!6 2.00000

8.17680e#01 !6.87600e#03 0.0 0.0 !8.70780 7.19260e!06 2.00000

7.44750e#01 !7.16430e#03 0.0 0.0 !7.32700 3.13400e!06 2.00000

7.69640E#01 !7.62380E#03 0.0 0.0 !7.49240 5.94360E!18 6.00000

Table 13 Physical property constants for the extended Antoine correlation Constants

Chloroform

1,2-Dichloroethane

2-Butoxyethanol

Ethyl acetate

C  C  C  C  C  C  C 

1.46430e#02 !7.79230e#03 0.0 0.0 !2.06140e#01 2.45780e!02 1.00000

92.355 !6920.4 0.0 0.0 !10.651 9.1426e!06 2.00000

1.50800E#02 !1.17280E#04 0.0 0.0 !1.88830E#01 1.1294E!05 2.00000

6.68240e#01 !6.22760e#03 0.0 0.0 !6.41000 1.79140e!17 6.00000

Table 14 Physical property constants for the extended Antoine correlation Constants

Acetone

THF

Benzene

Toluene

C  C  C  C  C  C  C 

6.90060e#01 !5.59960e#03 0.0 0.0 !7.09850 6.22370e!06 2.00000

54.898 !5305.4 0.0 0.0 !4.7627 1.4291e!17 6.00000

8.39180E#01 !6.51770E#03 0.0 0.0 !9.34530 7.11820E!06 2.00000

80.877 !6902.4 0.0 0.0 !8.7761 5.8034e!06 2.00000

J.E. Tolsma, P.I. Barton / Chemical Engineering Science 55 (2000) 3835}3853

Appendix B. Physical property constants for the extended Antoine correlation The extended Antoine Correlation used to compute vapor pressures in the calculations has the following form: C G #C #C ln ¹C ¹! , ln pQ "C # (20) G G ¹#C G G G G G where pQ is in Pascals and ¹ is in Kelvins. The parameter G values were obtained from Aspen Technology, Inc.'s Properties Plus (1995). The values used for the calculations presented in this paper are shown in Tables 12}14.

References Aspen Technology Inc., 1995. Aspen Plus reference manual release 9. Cambridge, MA 02139, USA. Baker, L. E., Pierce, A. C., & Luks, K. D. (1982). Gibbs energy analysis of phase equilibrium. Society of Petroleum Engineers Journal, 22, 731}742. Fidkowski, Z. T., Malone, M. F., & Doherty, M. F. (1993). Computing azeotropes in multicomponent mixtures. Computers and Chemical Engineering, 17, 1141}1155. Gmehling, J., Menke, J., Krafczyk, J., & Fischer, K. (1994a). Azeotropic data, Part I. Weinheim: VCH. Gmehling, J., Menke, J., Krafczyk, J., & Fischer, K. (1994b). Azeotropic data, Part II. Weinheim: VCH.

3853

Horsley, L. H. (1973). Azeotropic data III, Advances in Chemistry Series, vol. 116. American Chemical Society, Washington, DC. Michelsen, M. L. (1982). The isothermal #ash problem: I. Stability. Fluid Phase Equilibria, 9, 1}19. Moore, R. E. (1977). A test for existence of solutions to nonlinear systems. SIAM Journal on Numerical Analysis, 14, 611}615. Moore, R. E. (1979). Method and applications of interval analysis. Philadelphia: SIAM. Park, T., & Barton, P. I. (1996). State event location in di!erential algebraic models. ACM Transactions on Modelling and Computer Simulation, 6, 137}165. Rheinboldt, W. C, & Burkardt, J. V. (1983). A locally parameterized continuation process. ACM Transactions on Mathematical Software, 9, 215}235. Sorensen, J., & Arlt, W. (1979). Liquid}liquid equilibrium data collection. Binary systems, Part 1.. In: Behrens, D., & Eckermann, R, Chemistry Data Series, vol. V Frankfurt: Dechema. Stadtherr, M. A., Schnepper, C. A., & Brennecke, J. F. (1995). Robust phase stability analysis using interval methods. A.I.Ch.E. Symposium Series 304, 91, 356}359. Tolsma, J. E. (1999). Simulation and analysis of heteroazeotropic systems. Ph.D Thesis, Massachusetts Institute of Technology, Cambridge, MA, June. Tolsma, J. E., & Barton, P. I. (2000a). Computation of heteroazeotropes, Part I: Theory. Chemical Engineering Science, in press. Tolsma, J. E., & Barton, P. I. (2000b). Sensitivity of predicted equilibrium behavior with respect to property model parameters, in preparation. Watson, L. T., Billups, S. C., & Morgan, A. P. (1987). Algorithm 652. Hompack: A suite of codes for globally convergent homotopy algorithms. ACM Transactions on Mathematical Software, 13, 281}310.