Computing All Homogeneous and Heterogeneous Azeotropes in

Ind. Eng. Chem. Res. 1999, 38, 4901-4912

4901

Computing All Homogeneous and Heterogeneous Azeotropes in Multicomponent Mixtures Stanislaw K. Wasylkiewicz,* Michael F. Doherty, and Michael F. Malone Department of Chemical Engineering, University of Massachusetts, Amherst, Massachusetts 01003

A homotopy method for computing azeotropes in homogeneous multicomponent mixtures, developed by Fidkowski et al.,1 has been generalized to include heterogeneous liquids. The method together with an arc length continuation and a rigorous stability analysis and calculation of liquid-liquid equilibrium in multicomponent mixtures, developed by Wasylkiewicz et al.,2 gives an efficient and robust scheme for finding all homogeneous as well as heterogeneous azeotropes predicted by a thermodynamic model. A new, more general homotopy function has been proposed which reduces to the form proposed by Fidkowski et al.1 for homogeneous liquid mixtures. Unlike the homogeneous case, turning points can often be found on solution branches leading to the heterogeneous azeotropes. A few examples for heterogeneous systems are shown and compared with results obtained using the homogeneous models. Introduction Because the azeotropes and their types determine the distillation boundaries, the existence of azeotropes among a considered component set is of critical importance to the design of separation equipment.3 For reliable design we need to know all azeotropes in the mixture. The problem of computing temperatures, compositions, and the types (stable node, unstable node, or saddle) of all azeotropes predicted by thermodynamic models for nonideal, multicomponent, heterogeneous mixtures is complicated by nonlinearity and complexity of realistic vapor-liquid-liquid equilibrium (VLLE) descriptions, the presence of multiple solutions, and the constraints on the compositions. Unlike the homogeneous case, there are additionally spurious solutions for heterogeneous mixtures, which can be identified and eliminated only if reliable and rigorous stability analysis and calculation of LLE are performed. In the most reliable and robust method for azeotrope calculations, Fidkowski et al.1 applied the homotopy method together with an arc length continuation to find temperatures and compositions of all the azeotropes predicted by the model. In their method all homotopy branches are connected; i.e., a branch ending at an azeotrope containing c components starts at a bifurcation point at the solution branch of order c - 1 and eventually all branches are directly or indirectly connected to one of the pure component branches. There are no isolated solutions, for which it would be difficult to provide initial conditions and an additional mapping would be necessary.4 It makes the method very rigorous and efficient. The method has been developed for homogeneous mixtures, and for heterogeneous liquids it can provide spurious solutions or some azeotropes would not be found (see the examples). To overcome these limitations, we have generalized the homotopy method by suitable modification of the homotopy func* Corresponding author. Present address: AEA Technology Engineering Software, Hyprotech, 800, 707 8th Avenue SW, Calgary, Alberta, Canada T2P 1H5, E-mail: stan. [email protected]. Phone: (403) 520-6000. Fax: (403) 520-6060.

tion and by imbedding the rigorous stability analysis2 in the calculations of homotopy bubble point temperature. The modified method gives an efficient and robust scheme for finding all homogeneous as well as heterogeneous azeotropes predicted by a thermodynamic model. The stability analysis imbedded in the azeotrope calculation method is not restricted to the method developed by Wasylkiewicz et al.2 It can be replaced by other algorithms, e.g., homotopy algorithm,5 global optimization technique,6-8 interval analysis,9 simulated annealing algorithm,10 or quasi-Newton successive substitution method11 (see the second example). Harding et al.12 applied a global optimization method to find all homogeneous azeotropes in multicomponent mixtures. Their approach is based on developing convex underestimators for each of the examined thermodynamic models. The models are restricted to activity coefficient types. The method is limited only to homogeneous mixtures and is unable to find heterogeneous azeotropes. Also, the computational time required for finding an azeotrope increases rapidly with the number of components in the mixture. Maier et al.13 applied an interval Newton method with a generalized bisection algorithm to enclose all homogeneous azeotropes in multicomponent mixtures or to verify when none exists. The method guarantees finding all homogeneous azeotropes and does not require initial starting points or the construction of model-specific convex underestimating functions. It can be applied in connection with any thermodynamic model. The theoretical guarantees on solution enclosures come, however, at a very large cost in additional computation time relative to faster, but not guaranteed, local methods. Chapman and Goodwin14 applied the LevenbergMarquardt algorithm to find homogeneous azeotropes and then checked their stability using Michelsen’s15 implementation of the tangent plane method. Then they used unstable solutions as starting points in new searches for heterogeneous azeotropes. Unfortunately, heterogeneous azeotropes do not necessarily correspond to unstable solutions and vice versa. We show later in the paper that if a search for azeotropes is performed without a stability test (liquid assumed to be homogeneous in the whole composition space), only three

10.1021/ie990214t CCC: $18.00 © 1999 American Chemical Society Published on Web 12/06/1999

4902

Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999

azeotropes for the ternary mixture of n-propanol, water, and toluene can be found, whereas the rigorous heterogeneous model predicts four azeotropes. There is no unstable solution, which corresponds to the heterogeneous ternary azeotrope. Eckert and Kubicek16 have done a generalization of the method developed by Fidkowski et al.1 for the calculation of homogeneous azeotropes to treat heterogeneous mixtures. Analogous to Fidkowski’s method, they use a continuous deformation of equilibrium relations by changing an artificial parameter (homotopy parameter) from 0 (ideal case) to 1 (nonideal case). They “deform” vapor-liquid equilibrium (VLE) equations as well as liquid-liquid equilibrium (LLE) equations and then solve the deformed bubble point problem.17 They have developed an adaptive algorithm for switching between VLE and VLLE calculations to overcome difficulty in properly choosing the initial approximation for the Newton-Raphson method. If multiple solutions (stable and unstable) coexist (see e.g. Figure 9 in2), a rigorous stability test (e.g., Gibbs tangent plane test) is necessary to distinguish between these solutions. In our approach we gradually deform only the VLE equations and imbed a rigorous LLE calculation algorithm without any deformation. This makes it possible to use our rigorous phase stability algorithm2 without any modifications. Because homotopy paths can cross boundaries between homogeneous and heterogeneous liquid regions, a rigorous stability analysis must be performed if a mixture is prone to split into two or more liquid phases. In the paper, we compare some results obtained without stability calculations to those obtained with the rigorous ones.

Figure 1. T-x,y diagram for a mixture of heptane and methanol at 0.1 atm. NRTL model.

Motivation If a system is prone to form two or more liquid phases in some areas of composition space, it is very important to perform stability tests to check whether the selected liquids are stable or split into multiple liquid phases. Even if we expect that a liquid should be stable, the stability analysis in any equilibrium calculations is necessary to find the correct stable solution. The following examples explain the problem in details. Let us consider a dew point calculation algorithm for the following two mixtures: (1) Binary mixture of methanol (1) and heptane (2). NRTL model. Pressure, 0.1 atm. Vapor composition, y1 ) 0.72. (2) Ternary mixture of water (1), acetic acid (2), and butyl acetate (3). UNIQUAC model. Pressure, 1.0 atm. Vapor composition, y1 ) 0.7285 and y2 ) 0.0081. At first, we do VLE calculations for the binary mixture with homogeneous liquid assumed throughout the algorithm (i.e., stability test not used). The solution is as follows: temperature, T ) 11.11 °C; liquid composition, x1 ) 0.8599. If we apply the stability test, the solution will be as follows: T ) 12.10 °C, x1 ) 0.1268 (see Figure 1). Analogously, for the ternary mixture, the first (incorrect) solution is as follows: T ) 90.86 °C, x1 ) 0.18865; x2 ) 0.03726. The second (correct) solution is as follows: T ) 91.85 °C; x1 ) 0.98083; x2 ) 0.01721 (see Figure 2). It is worth noticing that all the solutions are homogeneous. But only the second ones are stable. It is wellunderstood that it is very unlikely to find multiple-liquid solutions in the dew point calculations. So why should we perform a stability test if we expect only one liquid

Figure 2. Isobaric liquid-liquid equilibrium surface for the mixture of water, acetic acid, and n-butyl acetate at 1 atm. UNIQUAC model. / indicates vapor compositions.

phase in the equilibrium? The answer is we still MUST have the stability test to find the correct homogeneous solution. Let us look at the details of the dew point calculation algorithm. There are two loops in the dew point calculation algorithm. In the outer loop we search for the dew point temperature. In the inner loop, the unknown liquid compositions are calculated. In Figure 3, circles correspond to each step of the homogeneous algorithm. We start with the initial guess: T ) 10.0 °C; x1 ) y1 ) 0.72. At first, the inner loop of the algorithm converges to x1 ) 0.8631 at T ) 10.0 °C. Then, in the outer loop, the temperature is changed to 11.0, 13.0, 11.11 °C, etc., while in the inner loop the liquid composition is slightly adjusted. The homogeneous algorithm converges to T ) 11.11 °C and x1 ) 0.8599. It is a perfect solution if we do not do any stability tests. However, if we do the test, we can reveal that it is not a stable solution. In

Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 4903

algorithm eventually converges to the stable, correct solution: T ) 91.85 °C; x1 ) 0.98083; x2 ) 0.01721. It is worth mentioning that if we start calculations with the initial guess T ) 90.0 °C, x1 ) y1, x2 ) y2, the homogeneous algorithm alone will converge to the correct solution. As we have seen in the above examples, even if the correct solution in dew point calculations is homogeneous, the calculations without the stability test can give an incorrect solution, which is thermodynamically unstable. The difference between solutions can be very large. Even if the correct solution is actually homogeneous, we still MUST have the stability test. Therefore, if we are searching for azeotropes in a mixture prone to form multiple-liquid phases, we MUST apply a rigorous stability test to make our search rigorous and complete. Computing Azeotropes in Heterogeneous Mixtures Figure 3. T-x,y diagram for a mixture of heptane and methanol at 0.1 atm. NRTL model. A circle corresponds to one step of the homogeneous algorithm and a square to the heterogeneous one.

In a c-component mixture, heterogeneous as well as homogeneous azeotropes can be calculated by solving the equation

f(z) ≡ y - z ) 0

(1)

subject to the constraints on mole fractions (x, y) and liquid-phase fractions (φ) c

c

xki ) 1, ∑yi ) 1, ∑ i)1 i)1

p

zi )

p

∑ φkxki , k)1 ∑ φk ) 1 k)1

xkl g 0, yi g 0, φk g 0,

(2)

i ) 1, ..., c, k ) 1, ..., p (3)

and the vapor-liquid equilibrium constraints:

yi )

Figure 4. Dew point calculations for the mixture of water, acetic acid, and n-butyl acetate at 1 atm. UNIQUAC model. A circle corresponds to one step of the homogeneous algorithm and square to the heterogeneous one.

that situation we perform rigorous bubble point calculations (with stability test), which provides us with a stable liquid composition (x1 ) 0.1520) as a new initial guess (see a square in Figure 3). Starting from this point, the dew point calculation algorithm eventually converges to the stable, correct solution: T ) 12.10 °C; x1 ) 0.1268. The solution path of the dew point algorithm for the ternary mixture is shown in Figure 4. We start with the initial guess: T ) 100.0 °C; x1 ) y1 ) 0.7285; x2 ) y2 ) 0.0081. At first, the homogeneous algorithm converges to the unstable solution: T ) 90.86 °C; x1 ) 0.18865; x2 ) 0.03726. Then, we perform rigorous bubble point calculations (with stability test), which provide us with a stable liquid composition (x1 ) 0.98076; x2 ) 0.01723) as a new initial guess (see a square in Figure 4). Starting from this point, the dew point calculation

γki (T, xk) Psi (T) k xi , φi(T,y)P

i ) 1, ..., c, k ) 1, ..., p (4)

z is a vector of the overall compositions of a liquid, which splits into p liquid phases, xk is a vector of the kth liquidphase compositions, and y is a vector of compositions of the vapor phase at equilibrium with each of the liquid phases. Fugacity (φ) and activity coefficients (γ) are nonlinear functions of temperature and compositions of the vapor or liquid phases. Pure component saturation pressures (Ps) are also strongly nonlinear functions of temperature. To find a temperature and composition of an azeotrope, one has to solve the above set of equations together with appropriate correlations for physical properties (fugacity, activity, and saturation pressure). This is a strongly nonlinear, constrained problem with multiple solutions. Not only azeotropes but also all pure components are solutions to this problem. There is no conventional root-finder, which can be used robustly to find all of the solutions for a multicomponent mixture. An efficient and robust scheme for finding all homogeneous (p ) 1) azeotropes predicted by a thermodynamic model has been proposed by Fidkowski et al.1 They introduced an artificial equilibrium relation

[

]

γi Psi x, y˜ i ) (1 - t) + t φi P i

i ) 1, ..., c

(5)

4904


way. It has also been shown1 that all solution branches are connected to the branches for which a starting point can be identified a priori (pure components). This makes the method very robust and guarantees finding all the homogeneous azeotropes predicted by a model. In this paper, the homotopy method for computing azeotropes in homogeneous mixtures has been generalized to include heterogeneous liquids. If the liquid is...(1) homogeneous, then p ) 1 and x in eq 4 is equal to z; (2) heterogeneous, then p > 1 and x is the composition of one liquid phase chosen from liquids into which the liquid of the overall composition z is splitting. The homotopy continuation method for heterogeneous liquids can be summarized as follows: (1) Homotopy starts at the ideal equilibrium described by Raoult’s law:

yideal ) i

Psi z, P i

i ) 1, ..., c

(6)

and ends at the final nonideal vapor-liquid or vaporliquid-liquid equilibrium:

yreal i

γiPsi ) x, φ iP i

i ) 1, ..., c

(7)

(2) Solutions are tracked as the equilibrium is deformed smoothly from ideal to real by changing the homotopy parameter t from 0 to 1 in the homotopy function:

h(z, t) ) y˜ - z

(8)

y˜ ) (1 - t)yideal + tyreal

(9)

where

Figure 5. (a) y˜ -x and (b) T-y,x diagrams at various values of the homotopy parameter t for a mixture of chloroform and methanol at 1 atm. Wilson model.

which gives Raoult’s law solution when the homotopy parameter t ) 0 and the real nonideal equilibrium relation when t ) 1. Examples of equilibrium lines given by the solutions of eq 5 for various homotopy parameters are shown in Figure 5. By changing homotopy parameter from 0 to 1, the equilibrium surface can be gradually deformed, beginning from the hypothetical ideal equilibrium described by Raoult’s law (no azeotropes; only pure component compositions are solutions of the set of equations (1)-(4)). An emerging azeotrope is then gradually “moved” from the point corresponding to the pure chloroform and ended at the final, nonideal equilibrium surface with the real azeotropic point. To start, there are as many branches of solutions as there are pure components and at least one of the branches will show a bifurcation when an azeotrope is present in the mixture. A homotopy continuation method has been applied1 to do this gradual deformation in a systematic

Starting points for the homotopy continuation method are known a priori. These are pure components, which are trivial solutions to the set of equations (1)-(4). Actually, they are the same as those in the case of the homogeneous model because pure component liquids cannot split into two or more liquid phases. During branch tracking the stability test has to be performed, and if liquid-phase splitting has been detected, not the overall liquid composition but the composition of a chosen liquid phase in equilibrium must be used in eq 7. To ensure proper switching between homogeneous and heterogeneous regions, the rigorous stability analysis and LLE calculation2 are imbedded in the calculations of homotopy bubble point temperature. Theoretically, there should be no difference in which one of the liquid phases is used to calculate yreal i . All of them are in liquid-liquid equilibrium and should give the same results. During the continuation calculations, however, these liquids not only change their composition but also appear and disappear. To avoid discontinuities in calculations and low accuracy for disappearing phases, we select the liquid phase, which is in the greatest amount in the mixture of liquids in the liquid-liquid equilibrium. Examples of equilibrium lines given by solutions of eq 9 for various homotopy parameters are shown in Figure 6. Changing the homotopy parameter from 0 to 1, we gradually deform the equilibrium lines beginning from the hypothetical ideal equilibrium described by Raoult’s law. There is no difference between the het-


slightly different and the temperature is almost 20 °C lower. The correct heterogeneous model applies the stability test to cut off the unstable phases represented by the curves below the liquid-liquid tie line in Figure 7c. Equilibrium lines for both the homogeneous and heterogeneous models calculated at t ) 0.5 are shown in Figure 7b. The parts of the liquid and vapor lines, which correspond to the homogeneous liquid region, are common for both models. But in the heterogeneous region the solutions of eq 9 are different because they combine a part from the homogeneous model (eq 6) with the heterogeneous one (eq 7). In the two-liquid region, the pseudo-liquid lines are not parallel to the composition axes and even are not straight lines. But they become more and more parallel and straight as the homotopy parameter approaches 1. The pseudo-vapor lines given by the solution of eq 9 and corresponding to the heterogeneous liquid region (Figure 6b) become shorter and shorter as the homotopy parameter increases and eventually disappear at t ) 1. Tracking Homotopy Branches Efficient tracking of all homotopy branches is crucial for the success of the azeotrope calculation method. One can attempt to track the branches of eq 8 originating from pure components, where they are located at t ) 0, by increasing gradually the value of the homotopy parameter and using the most recently calculated z as an initial guess for the next point on the branch. This type of parameter continuation fails, however, at turning points (or other singularities) where the homotopy parameter must start to decrease to follow the branch. Turning points are common occurrences in heterogeneous azeotrope calculations (see the examples). To resolve these difficulties, an arc length continuation has been used, where an additional equation defining the arc length g(z, t, s) has been added to the system of equations (8) to form an augmented system:1

G(z,t,s) )

Figure 6. (a) y˜ -x and (b) T-y,x diagrams at various values of the homotopy parameter t for a mixture of water and toluene at 1 atm. Heterogeneous UNIQUAC model.

erogeneous and homogeneous homotopy functions for t ) 0. Then, we gradually “add nonideality” to the system by increasing the value of the homotopy parameter. Because the water-toluene miscibility is so small, the deflection on the equilibrium lines can be seen, even for as small a value of the homotopy parameter as 0.05. Notice that only for t ) 1 the equilibrium line in Figure 6a is the straight line in the two-liquid region and it become more and more curved as the homotopy parameter decreases. The azeotrope emerges from the point corresponding to the pure water, gradually “moves” toward the toluene point as the homotopy parameter increases, and finally ends at the real heterogeneous azeotropic point. Examples of equilibrium lines given by solutions of the homogeneous model (eq 5) for various homotopy parameters are shown in Figure 7. The incorrect homogeneous model also predicts one azeotrope for the mixture of water and toluene, but the composition is

[

]

h(z,t) )0 g(z,t,s)

(10)

Equation 10 contains c unknowns (z1, z2, ..., zc-1, t) and one parameter s (arc length). The current branch can be tracked by gradually increasing parameter s and solving eq 10. Now, homotopy parameter t is one of the unknowns and can also decrease (e.g., near turning points) while the arc length always increases. The arc length relation4 c-1

∑ i)1

() () dzi

2

+

ds

dt

2

)1

ds

(11)

is a nonlinear function of the unknowns. In a more practical approach, for small step sizes, it is more useful to use, instead of (11), the pseudo arc length equation as defined by Keller:18 c-1

g(z, t, s) ) Θ

[zi - zi(sk)] ∑ i)1

dzi(sk)

+

ds

(1 - Θ)[t - t(sk)]

dt(sk) ds

- (s - sk) (12)

The tuning factor Θ ) 1 - (1/c) has been selected to

4906


Figure 7. (a) y˜ -x and (b and c) T-y,x diagrams at various values of the homotopy parameter t for a mixture of water and toluene at 1 atm. UNIQUAC model.

place more emphasis on composition variables than on the homotopy parameter. Note that eq 12 is a linear equation according to z and t because z(sk) and t(sk) are known from the previously calculated point (for an arc length sk) and the vector tangent to the homotopy path at point k [dzT(sk)/ds, dt(sk)/ds] can be approximated by the secant vector connecting points k - 1 and k. This vector is also used in a predictor step during branch tracking. In a corrector step, we solve eq 10 by a Newton-Raphson method including constraints (2)-(4) and an optimization of the Newton step length.1 The arc step size (∆s ) s - sk) is set initially to 0.01 and then reduced, when the corrector step has difficulties with convergence (e.g., for a strongly nonlinear branch), or increased for close-to-linear branches, far from any bifurcation point. At each point we check for bifurcation. If the determinant of the Jacobian matrix of eq 10 changes sign between two consecutive points on the branch, we calculate the bifurcation point (where Ja-

cobian is singular) precisely using a secant method. At the bifurcation point, two continuation branches intersect: the currently tracked branch (leading to an n-component azeotrope) and a new one of the higher order (n+1) (leading to an (n+1)-component azeotrope). The bifurcation point is stored together with initial direction of the new branch [dzT0 /ds, dt0/ds]. This direction is computed from the following approximation of the algebraic bifurcation equation,1,18

[ ][ ] [ ] [ ][ ] ∂hi ∂zj dzTb ds

∂hi ∂t dtb ds

dz0 0 ds dt0 ) 0 ds

(13)

where the subscript b indicates derivatives computed at the bifurcation point for the old branch and the subscript 0 for the new branch. Initial directions for


Figure 8. Branches of solutions (∇) for binary and ternary azeotropes on a composition triangle and the isobaric liquid-liquid tie lines with corresponding vapor compositions (/) for a mixture of n-propanol, water, and toluene at 1 atm. NRTL model.

binary azeotropes are calculated by analyzing boiling temperatures of pure components and azeotropes in binary systems (see ref 1 for details). For ternary and higher order azeotropes, the starting points and initial directions are determined at the corresponding bifurcation points at the lower order homotopy branches. The azeotrope search procedure is finished when all initial conditions have been explored for branch tracking up to t ) 1. At the beginning we know initial conditions only for binary branches. Then, during branch tracking, each new bifurcation point provides additional initial conditions, which is stored for later exploration. Consistency Test The eigenvalues of the Jacobian matrix [-∂hi/∂zj] at a solution t ) 1 give the stability (type) of an azeotrope or a pure component.1 When the azeotropes’ types are known, the topological consistency of the residue curve map (RCM) can be checked via the Zharov-Serafimov topological constraint:19 c

c-1 + 2k(N+ +1 ∑ k + Sk - Nk - Sk ) ) (-1) k)1

(14)

+ where N+ k and Sk are the numbers of k-component nodes and saddles, respectively, with index +1, while Nk and Sk are the corresponding singular points with index -1. If the number of negative eigenvalues of the Jacobian matrix is even, then the index is equal to +1. If the number is odd, then the index is -1.

Example 1: Ternary Mixture Branches of solutions for three binary and one ternary azeotropes and the isobaric liquid-liquid miscibility gap are shown on a composition triangle in Figure 8 for a mixture of n-propanol, water, and toluene at a pressure of 1 atm. The triangle symbols represent steps of the branch-tracking algorithm and point into the calculation direction. A total CPU time used for this example was

Figure 9. Branches of solutions for binary and ternary azeotropes on a composition triangle for a mixture of n-propanol (P), water (W), and toluene (T) at 1 atm. (a) Homogeneous NRTL model. (b) Heterogeneous NRTL model.

8 s on an Alpha Station 250. The calculations were performed using the NRTL model and the rigorous stability test (see also Figure 9b). For comparison, in Figure 9a, the branches of continuation solutions are shown which have been obtained under the assumption that liquid phases are homogeneous.1 Of course, the assumption is not correct for this system and has been made only for comparison purposes. The homogeneous model does not predict the existence of a ternary azeotrope. For the homogeneous model, there are two bifurcation points on binary branches WP (leading to the water-n-propanol azeotrope) and WT (leading to the water-toluene azeotrope) at which the ternary PWT branch is crossed. The bifurcation points are connected by the PWT branch and no ternary azeotrope has been detected at this branch. Notice that there is no bifurcation point at the WT branch (Figure 9b) for the rigorous heterogeneous model. The bifurcation point at the branch WP is the same as that for the homogeneous

4908


Figure 10. (a) Composition and (b) temperature bifurcation diagrams for a mixture of n-propanol (P), water (W), and toluene (T) at 1 atm. Heterogeneous NRTL model. One ternary azeotrope (PWT) has been found.

model and gives rise to the ternary branch PWT, which after several steps in the homogeneous region (see Figure 8), turns toward and ends in the heterogeneous ternary azeotrope. The turning point at the border between the homogeneous and heterogeneous regions is very noticeable in Figure 10 for the composition (a) and temperature (b) bifurcation diagrams. For homogeneous mixtures, branches with a decreasing homotopy parameter have been ruled out from further consideration1 because as a rule they merge with other branches attainable by forward branch-tracking calculations. No turning points have been detected for homogeneous mixtures. But for heterogeneous mixtures this is not the case. The distillation region diagrams (DRD) are shown in Figure 11 for the homogeneous (a) and heterogeneous (b) models. The structures of the diagrams, the number of azeotropes, and the stability of the water-toluene (WT) azeotrope are different. For the homogeneous

Figure 11. The distillation region diagrams for a mixture of n-propanol (P), water (W), and toluene (T) at 1 atm. (a) Homogeneous NRTL modelsno ternary azeotrope. (b) Heterogeneous NRTL modelsone ternary minimum boiling azeotrope.

model, the WT azeotrope is a minimum boiling one (unstable node) in Figure 11a whereas it becomes a saddle (intermediate boiling) azeotrope for the heterogeneous model (Figure 11b). Notice that both DRDs are topologically consistent because the Zharov-Serafimov topological constraint19 is fulfilled in both cases (see Table 1). The DRD structure gives rise to three separation regions for both models but they are very different. In Figure 11a, two regions are determined by one saddle, one unstable node, and one stable node, whereas in Figure 11b, all regions are defined by two saddles, one unstable mode, and one stable node. A separation system designed on the grounds of the homogeneous DRD (Figure 11a) could be very different from the system design using the correct heterogeneous DRD

Ind. Eng. Chem. Res., Vol. 38, No. 12, 1999 4909 Table 1. Results of Azeotrope Calculations (Pressure, 1 atm; NRTL Model) components and azeotropes

azeotrope composition (mole fraction)

boiling point (°C)

type

Homogeneous Model n-propanol (P) water (W) toluene (T) WT WP PT

0.0000 0.4096 0.6671

0.5550 0.5904 0.0000

97.15 100.00 110.61 72.05a 90.19 95.29

0.4450 0.0000 0.3329

index

SN SN SN UN SA SA

+1 +1 +1 +1 -1 -1

SN SN SN SA SA SA UN

+1 +1 +1 -1 -1 -1 +1

2 3-1 + + + N+ +1 1 ) 3, N2 ) 1, S2 ) 2, 2N1 + 2 (N2 - S2 ) ) 2 ) (-1)

Heterogeneous Model n-propanol (P) water (W) toluene (T) WT WP PT PWT

0.0000 0.4096 0.6671 0.1807

97.15 100.0 110.61 84.52 90.19 95.29 81.99

0.4414b 0.0000 0.3329 0.3408

0.5586 0.5904 0.0000 0.4785

2 3 + 3-1 + + N+ +1 1 ) 3, S2 ) 3, N3 ) 1, 2N1 - 2 S2 + 2 N3 ) 2 ) (-1) a

Underlines help to see differences between models. b Heterogeneous azeotrope composition printed in bold.

Table 2. Results of Azeotrope Calculations (Pressure, 1 atm; NRTL Model) components and azeotropes ethanol (E) water (W) cyclohexane (C) benzene (B) EW EC EB WC WB CB EWC EWB ECB WCB EWCB

azeotrope composition (mole fraction)

0.9187 0.4758 0.4435 0.0000 0.0000 0.0000 0.3098 0.2722 0.4326 0.0000 0.2672

0.0813 0.0000 0.0000 0.2977 0.2860 0.0000 0.1702 0.1891 0.0000 0.2703 0.1771

0.0000 0.5242 0.0000 0.7023 0.0000 0.4314 0.5200 0.0000 0.3493 0.3232 0.2816

boiling point (°C)

0.0000 0.0000 0.5565 0.0000a 0.7140 0.5686 0.0000 0.5387 0.2181 0.4065 0.2741

78.17 100.00 80.68 78.75 78.09 66.51 67.30 69.38 68.33 76.91 63.33 63.29 66.07 67.10 62.28

type SN SN SN SN SA SA SA SA SA SA SA SA SA SA UN

index -1 -1 -1 -1 +1 +1 +1 +1 +1 +1 -1 -1 -1 -1 +1

+ + N1 ) 4, S2 ) 6, S3 ) 4, N4 ) 1

2 + 3 4 + 4-1 -2N+1 1 + 2 S2 - 2 S3 + 2 N4 ) 0 ) (-1) a

Heterogeneous azeotrope composition printed in bold.

(Figure 11b). Of course, only the last system would work for the real heterogeneous liquid separation. Example 2: Quaternary Mixture There are 11 azeotropes in the quaternary mixture of ethanol (E), water (W), cyclohexane (C), and benzene (B) at 1 atm: 1 quaternary, 4 ternary, and 6 binary azeotropes (see Table 2). Binary WC and WB, ternary EWC, EWB, and WCB, and quaternary EWCB azeotropes are heterogeneous (printed in bold). The other five azeotropes are homogeneous. Boiling point temperatures, types, and indices of all azeotropes and pure components are shown in Table 2. The results of the calculations are topologically consistent because the Zharov-Serafimov topological constraint19 is fulfilled (see calculations on the bottom of Table 2). The mixture has a wide miscibility gap. Four ternary faces of a composition tetrahedron are shown in Figure

12 together with all binary and ternary azeotropes, distillation boundaries, liquid-liquid equilibrium regions, and vapor lines. Branches of solutions for all 11 azeotropes are shown on the composition tetrahedron in Figure 13a and on the temperature bifurcation diagram in Figure 13b. The calculations were performed using the NRTL model and the rigorous stability test. For this example, the homotopy function defined by eqs 6-9 has been generalized to include any thermodynamic model, not just the activity coefficients one. The vapor mole fractions have been calculated using equilibrium ratios (Ki),

yideal ) Kideal zi, i i

i ) 1, ..., c

(15)

) Kreal yreal i i xi,

i ) 1, ..., c

(16)

where the Kreal values are expressed as a function of

4910


Figure 12. Ternary faces of a composition tetrahedron for a mixture of ethanol (E), water (W), cyclohexane (C), and benzene (B) at 1 atm. Note azeotropes, distillation boundaries, liquid-liquid equilibrium regions, and vapor lines.

fugacity coefficients in vapor (φV) and liquid (φL) phases:

Kreal ) i

φVi (x, P, T) , φLi (y, P, T)

) Kideal i

Psi , P

i ) 1, ..., c

i ) 1, ..., c

(17)

(18)

The stability analysis11 imbedded in the azeotrope calculation method is based on Michelsen’s15 implementation of the tangent plane criterion and applies equilibrium ratios in their problem formulation and quasiNewtom successive substitution method to solve the equations. There are no bifurcation points on the binary branches EW in Figure 13 (leading to the ethanol-water azeotrope), CW (leading to the cyclohexane-water azeotrope), and BC (leading to the benzene-cyclohexane azeotrope). A bifurcation point on the binary branch EC gives rise to a new ternary branch, which ends at the

ECW azeotrope. At a bifurcation point on the binary branch BW, a new ternary branch starts and after several steps ends in the BWC azeotrope. The binary branch EB bifurcates twice to two ternary branches EBC and EBW. The last one bifurcates again in the point where the quaternary EBWC branch starts. Two turning points on the branches can be seen in Figure 13b, one on the ECW branch and one on the EBW branch. The turning points can be associated with crossing the borders between homogeneous and heterogeneous regions. The NRTL model has been used to describe liquidphase behavior. Binary interaction parameters for W-B and W-C have been taken from ref 20 for the E-C pair from ref 21 and for other pairs from ref 22. The parameters are summarized in Table 3. A total CPU time used for this example was 40 s on a PC (Pentium II, 400 MHz) running Windows NT. Maier et al.13 applied the interval method to find all homogeneous azeotropes in the ternary system benzeneethanol-water (a subsystem of our quaternary system)


did not find any ternary EBW azeotrope. It is correct because they used the homogeneous model, which does not predict any ternary azeotrope (similar to that in Figure 11a for the n-propanol-toluene-water mixture). On the contrary, the heterogeneous model describes the system correctly and predicts heterogeneous, ternary azeotrope EBW (see Table 2), similar to that of the ternary azeotrope in Figure 11b for the n-propanoltoluene-water system. Conclusions Knowledge of temperatures and compositions of all azeotropes in a multicomponent mixture is crucial for the design of distillation systems. A homotopy method for computing azeotropes in homogeneous multicomponent mixtures has been generalized for heterogeneous liquids. The method together with an arc length continuation and a rigorous stability analysis in multicomponent mixtures gives an efficient and robust scheme for finding all homogeneous as well as heterogeneous azeotropes predicted by a thermodynamic model. The method should be applied to all mixtures, which have multiple-liquid regions, because the homogeneous method is prone to miss some heterogeneous azeotropes or predict wrong type of the azeotropes. Consequently, the structure of the distillation region diagram can be determined incorrectly, giving rise to incorrect design of distillation systems. The topological constraint (eq 6) fulfillment does not guarantee that the DRD has been determined correctly. A rigorous stability analysis must be performed if a mixture is prone to split into two or more liquid phases in some regions of concentration space. Acknowledgment We are grateful for financial support from E. I. duPont Co., Eastman Kodak Co., ICI, and Union Carbide Corporation. Figure 13. (a) Branches of solutions for binary, ternary, and quaternary azeotropes on a composition tetrahedron and (b) temperature bifurcation diagram for a mixture of ethanol (E), water (W), cyclohexane (C), and benzene (B) at 1 atm. Heterogeneous NRTL model. Table 3. Binary Interaction Parameters for a Mixture of Ethanol (E), Water (W), Cyclohexane (C), and Benzene (B) (NRTL Model) aEW ) 1332.312a aWE ) -109.6339 aCE ) 790.5410 aBE ) 334.1524

aEC ) 993.6688 aWC ) 4069.410 aCW ) 2217.530 aBW ) 3719.395

aEB ) 991.5067 aWB ) 4843.376 aCB ) -33.39469 aBC ) 286.9768

REW ) 0.3031 RWC ) 0.2000

REC ) 0.4626 RWB ) 0.2000

REB ) 0.2911 RCB ) 0.3025

a

Binary interaction parameters aij are in cal/mol.12

at 1 atm using the UNIQUAC model. Their results are very close to our calculations for the binary azeotropes EB and EW, despite a different thermodynamic model. There are, however, significant differences for the binary BW azeotrope both in the temperature and the composition. Because Maier et al.13 used a homogeneous model, they found correct solutions for the homogeneous azeotropes EB and EW, but not for the heterogeneous BW azeotrope. They found a spurious homogeneous solution (similar to the one shown in Figure 7c for the toluenewater mixture), which has a much lower boiling temperature than the real heterogeneous azeotrope. They

Nomenclature c ) number of components G ) augmented homotopy function g ) arc length equation h ) homotopy function Ki ) equilibrium ratio of component i N+ k ) number of k-component nodes with index +1 Nk ) number of k-component nodes with index -1 P ) pressure Ps ) saturation pressure p ) number of liquid phases in equilibrium s ) arc length S+ k ) number of k-component saddles with index +1 Sk ) number of k-component saddles with index -1 SA ) saddle SN ) stable node t ) homotopy parameter T ) temperature UN ) unstable node xi ) mole fraction of component i in the liquid phase x ) vector of mole fractions in the liquid phase yi ) mole fraction of component i in the vapor phase y ) vector of mole fractions in the vapor phase y˜ i ) pseudo-mole fraction of component i in the vapor phase y˜ ) vector of pseudo-mole fractions in the vapor phase zi ) overall liquid mole fraction of component i z ) vector of overall liquid mole fractions

4912


γ ) activity coefficient φ ) fugacity coefficient φ ) liquid-phase fraction Θ ) tuning factor in pseudo arc length equation

Literature Cited (1) Fidkowski, Z. T.; Malone, M. F.; Doherty, M. F. Computing Azeotropes in Multicomponent Mixtures. Comput. Chem. Eng. 1993, 17, 1141-1155. (2) Wasylkiewicz, S. K.; Sridhar, L. N.; Malone, M. F.; Doherty, M. F. Gibbs Stability Analysis and Calculation of Liquid-Liquid Equilibrium in Multicomponent Mixtures. Ind. Eng. Chem. Res. 1996, 35, 1395-1408. (3) Wasylkiewicz, S. K.; Kobylka, L. C.; Satyro, M. A. Design of Complex Azeotropic Distillation Systems Using Geometric Methods. In Proceedings of the Chemical Engineering Chemputers Northeast Conference and Exhibition, Philadelphia, PA, 1996; Chemical Engineering, A Division of the McGraw-Hill Companies: New York, 1996. (4) Seader, J. D.; Kuno, M.; Lin, W. J.; Johnson, S. A.; Unsworth, K.; Wiskin, J. W. Mapped Continuation Methods for Computing all Solutions to General System of Nonlinear Equations. Comput. Chem. Eng. 1990, 14, 71-85. (5) Sun, A. C.; Seider, W. D. Homotopy-Continuation Method for Stability Analysis in the Global Minimization of the Gibbs Free Energy. Fluid Phase Equilibr. 1995, 103, 213-249. (6) McDonald, C. M.; Floudas, C. A. Decomposition Based and Branch and Bound Global Optimization Approaches for the Phase Equilibrium Problem. J. Global Optim. 1994, 5, 205-251. (7) McDonald, C. M.; Floudas, C. A. Global Optimization for the Phase Stability Problem. AIChE J. 1995, 41, 1798-1814. (8) McDonald, C. M.; Floudas, C. A. Global Optimization and Analysis for the Gibbs Free Energy Function using the UNIFAC, Wilson and ASOG Equations. Ind. Eng. Chem. Res. 1995, 34, 1674-1687. (9) Hua, J. Z.; Brennecke, J. F.; Stadtherr, M. A. Reliable Prediction of Phase Stability using an Interval Newton Methodol. Fluid Phase Equilibr. 1996, 116, 52-59. (10) Zhu, Y.; Xu, Z. A Reliable Prediction of the Global Phase Stability for Liquid-Liquid Equilibrium through the Simulated

Annealing Algorithm: Application to NRTL and UNIQUAC Equations. Fluid Phase Equilibr. 1999, 154, 55-69. (11) Nghiem, L. X.; Li, Y. K.; Heidemann, R. A. Application of the Tangent Plane Criterion to Saturation Pressure and Temperature Computations. Fluid Phase Equilibr. 1985, 21, 39-60. (12) Harding, S. T.; Maranas, C. D.; McDonald, C. M.; Floudas, C. A. Locating All Homogeneous Azeotropes in Multicomponent Mixtures. Ind. Eng. Chem. Res. 1997, 36, 160-178. (13) Maier, R. W.; Brennecke, J. F.; Stadtherr, M. A. Reliable Computation of Homogeneous Azeotropes. AIChE J. 1998, 44, 1745-1755. (14) Chapman, R. G.; Goodwin, S. P. A General Algorithm for the Calculation of Azeotropes in Fluid Mixtures. Fluid Phase Equilibr. 1993, 85, 55-69. (15) Michelsen, M. L. The Isothermal Flash Problem: Part I. Stability. Fluid Phase Equilibr. 1981, 9, 1-19. (16) Eckert, E.; Kubicek, M. Computing Heterogeneous Azeotropes in Multicomponent Mixtures. Comput. Chem. Eng. 1997, 21, 347-350. (17) Eckert, E.; Kubicek, M. Bubble and Dew Point Calculations for Multiple Liquid Cases. Comput. Chem. Eng. 1995, 19, 493494. (18) Keller, H. B. Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems. In Application of Bifurcation Theory; Rabinowitz, P. H., Ed.; Academic Press: New York, 1977. (19) Zharov, W. T.; Serafimov, L. A. Physicochemical Fundamentals of Distillation and Rectification; Khimiya: Leningrad, 1975. (20) Sorensen, J. M.; Arlt, W. Liquid-Liquid Equilibrium Data Collection; DECHEMA: Frankfurt am Main, Germany, 1979. (21) Connemann, M.; Gaube, J.; Karrer, L.; Pfennig, A.; Reuter, U. Measurement and Representation of Ternary Vapour-LiquidLiquid Equilibria. Fluid Phase Equilibr. 1990, 60, 99-118. (22) Gmehling, J.; Onken, U.; Arlt, W. Vapor-Liquid Equilibrium Data Collection; DECHEMA: Frankfurt am Main, Germany, 1981.

Received for review March 22, 1999 Revised manuscript received September 20, 1999 Accepted September 27, 1999 IE990214T

Computing All Homogeneous and Heterogeneous Azeotropes in

Recommend Documents