- Email: [email protected]

Contents lists available at ScienceDirect

Fluid Phase Equilibria journal homepage: www.elsevier.com/locate/fluid

Fluid phase stability and equilibrium calculations in binary mixtures Part II: Application to single-point calculations and the construction of phase diagrams A. Giovanoglou, C.S. Adjiman, G. Jackson, A. Galindo ∗ Department of Chemical Engineering, Centre for Process Systems Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, United Kingdom

a r t i c l e

i n f o

Article history: Received 29 April 2008 Received in revised form 13 August 2008 Accepted 29 August 2008 Available online 18 September 2008

a b s t r a c t The methodology presented in Part I of this work is applied to a large number of pressure–temperature ﬂash calculations, and to the automated construction of constant temperature pressure–composition phase diagrams, and constant pressure temperature–composition phase diagrams for binary mixtures modeled with an augmented van der Waals equation of state. An automated prototype implementation of the algorithm is developed for this purpose. We follow the classiﬁcation of Scott and van Konynenburg [R.L. Scott, P.H. van Konynenburg, Discuss. Faraday Soc. 49 (1970) 87] and present phase diagrams corresponding to non-azeotropic mixtures of the ﬁve main types of ﬂuid phase behavior (I–V), studying in detail representative diagrams at constant pressure and constant temperature. Special attention is given to the solution of numerically problematic equilibrium regions, such as those close to three-phase equilibria where metastable and unstable critical points can also be found. Of the order of 104 ﬂash calculations at varying temperatures and pressures, and for different intermolecular parameters of the components in the mixture, have been carried out. The algorithm provides the correct stable equilibrium state for all of the points considered. Despite the fact that our implementation is not optimised for performance, we ﬁnd that the algorithm identiﬁes the stable solution in difﬁcult regions of the phase space without any penalty in terms of computational time, when compared to simpler regions. © 2008 Elsevier B.V. All rights reserved.

1. Introduction The ﬂuid phase equilibrium problem is often considered as a ﬂash separation at given pressure P, temperature T, and total composition z. However, the construction of phase diagrams over ranges of temperature, or pressure, and composition, and in particular the study of global phase behavior in the pressure–temperature P–T space is also of important practical and theoretical interest. For example, in process design the P–T projection of the phase diagram can be used to analyze optimal conditions and solvents for separation [1,2]. From a theoretical stand-point, global phase diagrams, where the generic types of ﬂuid phase behavior described by an equation of state are represented in terms of the parameter space of the intermolecular potential model, provide an immediate view of the capabilities of a chosen model and equation of state. The study of ﬂuid phase behavior in mixtures from equation of state approaches and its representation in phase diagrams dates back to the Dutch school under van der Waals and his co-workers

∗ Corresponding author. Tel.: +44 207 594 5606. E-mail address: [email protected] (A. Galindo). 0378-3812/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.ﬂuid.2008.08.018

[3]. The ﬁrst systematic classiﬁcation was presented by Scott and van Konynenburg [4] who studied binary mixtures with the van der Waals equation of state and simple mixing rules to map a global phase diagram. Using pressure–temperature projections of the full pressure–temperature–composition surface and studying the topology (the number, type, and continuity) of the critical lines, they were able to identify ﬁve general types of phase behavior (or types of mixtures) for different regions of the intermolecular parameter space. A key result was the ﬁnding that the van der Waals equation of state for mixtures yields ﬁve out of the six types of phase behavior observed experimentally (types I–V). Shortly after this work, Furman and Grifﬁths [5] identiﬁed the “shield region” of the phase diagram, which helped in the understanding of the boundaries between the different types of behavior of the global phase diagram for van der Waals mixtures. Scott and van Konynenburg acknowledged the existence of a sixth type of phase behavior (type VI) which is related to the existence of a closed-loop liquid–liquid equilibrium region and is invariably associated with mixtures that exhibit hydrogen bonding. It is generally accepted that some type of directional interaction [3,6–9], or a temperature-dependent attractive parameter, need to be incorporated in the theoretical model in order to obtain closed-loop behavior. It is also important to note in

96

A. Giovanoglou et al. / Fluid Phase Equilibria 275 (2009) 95–104

this context that type VI, and a new type VII ﬂuid phase behavior also associated with closed-loop liquid–liquid immiscibility, have now been obtained with equations of state of the van der Waals form [10,11]; these new types of behavior are found for a very small region of the parameter space in the global phase diagram. Though useful in understanding the transformations between the various types of phase behavior, the temperature states at which the loops appear (particularly the lower boundary of immiscibility) would be well within the solid region of the phase diagram, and one would expect such ﬂuid phases to be metastable with respect to crystallization. Following the work of Scott and van Konynenburg [4,12], the global phase diagrams of binary ﬂuid mixtures have been determined with a number of other equations of state [13–17] and new types of phase behavior have been suggested. As far as we are aware only types I–VI have been conﬁrmed experimentally to date. A new classiﬁcation has been presented to include the new types of phase behavior [18]. In particular the global phase behavior of binary mixtures with the augmented van der Waals (AVDW) equation of state (EoS) [19–22], which combines the Carnahan–Starling expression for hard spheres [23] with an attractive term at the van der Waals mean-ﬁeld level in the one-ﬂuid approximation, has been studied in some detail [10,11,24,25]. Yelash and Kraska [10] and Scott [11] studied binary mixtures of molecules of equal size modeled with the AVDW EoS, and calculated the corresponding global phase diagram including type VI phase behavior in a small region of the parameter space. Kolafa et al. [24,25] studied mixtures with varying diameter ratios of the components using the same theory. The work of Wang and Sadus [26] provides a link between the generic types of phase behavior obtained with the AVDW EoS and experimental systems. Using only pure component properties and standard mixing rules to map the properties of mixtures on to those of an effective pure component system, they provide a comparison of the experimental and calculated phase behavior, ﬁnding that the global phase diagram of the AVDW EoS can be used to predict the qualitative phase behavior of experimental systems. This study is of interest as it conﬁrms the fact that the AVDW EoS, although quite simple, is of practical use in the study of the phase behavior of real mixtures. The algorithm presented in Part I of this work [27] can be used for a robust and efﬁcient investigation of binary mixture ﬂuid phase behaviour through P–T ﬂash calculations, and the construction of constant-temperature or constant-pressure slices of the phase surface. An automated GAMS [28] code has been developed, and we use it to assess the performance of the phase equilibrium approach in these types of calculations with the AVDW EoS. We also examine the existence of metastable and unstable equilibrium states which appear close to the conditions of stable three-phase equilibria. Imre and Kraska [29] have also presented an interesting analysis of the limits of stability in binary mixtures with the AVDW expressions. The projections of the phase diagram in the P–T space can also be obtained directly by calculating the critical lines and end-points of the mixture. A recent algorithm for such a calculation has been presented by Cismondi and Michelsen [30]. Though the determination of critical lines from an analysis of the full pressure–temperature–composition surface is less efﬁcient than these types of direct algorithms, an advantage of the detailed approach is that one also obtains the full ﬂuid phase behavior. The methodology followed in our current work is described in the next section. We review the ﬁve main types (I–V) of phase behavior in the classiﬁcation of Scott and van Konynenburg [4], building the P–T projections and studying in detail constant-pressure and constant-temperature phase diagrams for the AVDW EoS. We consider temperatures above 0.3 of the critical point of the heavy component, which is expected to be the limit of stability of ﬂuid phases; types VI and above appear for the AVDW equation only

if temperatures below this cut-off are considered. Our implementation follows the proposed three-step algorithm for the case of non-azeotropic mixtures presented in Part I of this work [27]. The key thermodynamic properties and derivatives for the AVDW EoS are implemented analytically. The GAMS/CONOPT 2 solver is used for the solution of each individual step. The GAMS/CONOPT 2 solver employs a Generalized Reduced Gradient (GRG) local search algorithm and is suitable for the solution of highly non-linear optimization problems as well as for the solution of non-linear square systems of equations. 2. Methodology 2.1. P–T ﬂash calculations In a single-point phase equilibrium calculation, the input speciﬁcations are the scaled parameters of the AVDW EoS that describe the binary mixture, namely the ratios of the van der Waals attractive constants ˛1,1 /˛2,2 and ˛1,2 /˛2,2 (in the following we assume that the two species are of the same diameter), and the conditions of ∗ , the reduced pressure P ∗ , interest: the reduced temperature Tspec spec and the total composition of the mixture z1 . The reduced quantities are deﬁned in Part I [27]. In a P–T ﬂash calculation the algorithmic steps are executed sequentially only once and the default output information includes the number and types of phases in equilibrium, the molar fraction of each phase k in the mixture, wk , and the composition, xik (i = 1, 2), and density, k , of each phase. Since the algorithm is based on an exhaustive analysis of various thermodynamic spaces some additional information is also available, including all of the limits of mechanical and material stability of the mixture at the speciﬁed temperature and pressure, the limits of mechanical stability of any subcritical pure component at the speciﬁed temperature, and the vapor pressure of one or both components if the speciﬁed pressure belongs to an appropriate pressure interval. 2.2. Isothermal pressure–composition phase diagrams The structure of the algorithm presented in Part I easily allows for the automatic construction of phase diagrams, because during a P–T ﬂash calculation, the entire density–composition pattern is analyzed in terms of material and mechanical stability (step II of the algorithm [27]) and bounds are placed on all possible stable and metastable phases. Furthermore, in step IIIa, all possible twophase equilibrium calculations are performed between the existing phases, and, more importantly, all of the stable two- and threephase regions are identiﬁed before the total composition of the mixture is ﬁxed. This means that all of the necessary information for the construction of a phase diagram is available. The input speciﬁcations for the construction of a pressure– composition phase diagram are the scaled parameters of the AVDW EoS that describe the binary mixture, namely ˛1,1 /˛2,2 and ˛1,2 /˛2,2 , as before, and the ﬁxed reduced temperature of inter∗ . Optional input speciﬁcations include the lower (P ∗ ) and est Tspec low ∗ ) bounds that deﬁne the pressure range over which the upper (Pup phase diagram is constructed and the number of point calculations (l) to be performed over the speciﬁed pressure range. Here l is chosen to be 500. Considering that step I of the algorithm depends solely on the ∗ . For this speciﬁed temperature, it is executed only once at Tspec reason, it is considered as a preprocessing step. If bounds on the pressure are not speciﬁed, they are generated automatically in a way that ensures that the entire vapor–liquid equilibrium region(s) and a substantial part of the liquid–liquid equilibrium region are reported in the output information.

A. Giovanoglou et al. / Fluid Phase Equilibria 275 (2009) 95–104

In Part I, we showed that the pressure–density–composition surface can be classiﬁed in terms of 4 different types of topology; this should not be confused with the Scott and van Konynenburg classiﬁcation of phase behavior discussed earlier. For a type 1, 2 or 3 surface (Figs. 2–4, respectively, in [27]), the lower bound in pressure is set to be equal to the pressure corresponding to the high+ density limit of mechanical stability of component 2 (Pmech (0)). If the pressure–density–composition surface is of type 1 or 2, the upper bound on pressure is set to be twice the pressure corresponding to the low-density limit of mechanical stability of component 1 − (Pmech (1)), while if the pressure–density–composition surface is of type 3, then it is set to be twice the pressure at the pseudo-critical point (Pp−c ). Finally, if the pressure–density–composition surface is of type 4 (Fig. 5 in [27]), phase separation is observed only at elevated pressures and always above the critical pressure of the heavy component. In the context of the algorithm and given the scaling applied to the AVDW EoS, component 2 is always the heavy com∗ = 0.0044. In ponent and its critical pressure in reduced units is Pc,2 this case the lower bound of the pressure range is set to the critical pressure and the upper bound is set to be three times that of the lower bound. To construct the phase diagram, the algorithm is made to proceed through l iterative executions of steps II and IIIa. The default output information for the construction of a pressure–composition phase diagram consists of a set of data ﬁles that contain the pressure and the corresponding equilibrium compositions and densities for each of the stable two- or three-phase regions. The structure of the algorithm allows all the limits of mechanical and material stability to be made available as well. The procedure followed for the construction of a pressure–composition phase diagram is summarized in Table 1. 2.3. Isobaric temperature–composition phase diagrams In the construction of temperature–composition phase diagrams at constant pressure, the iterative procedure includes all three steps of the algorithm presented in Part I at each temperature. The input speciﬁcations are the scaled parameters of the AVDW EoS that describe the binary mixture, namely ˛1,1 /˛2,2 and ˛1,2 /˛2,2 , ∗ , the lower as before, and the ﬁxed reduced pressure of interest Pspec ∗ ) and upper (T ∗ ) bounds that deﬁne the temperature range (Tlow up over which the phase diagram will be constructed, and ﬁnally the number of point calculations (l) to be performed over the temperature range. An automated generation of the temperature range cannot be achieved in this case. From a physical point of view, the Table 1 Construction of pressure–composition phase diagrams. 1. Read input speciﬁcations: (a) parameters for the AVDW equation of state: ˛1,1 /˛2,2 and ˛1,2 /˛2,2 ∗ (b) reduced temperature: Tspec ∗ , l) (c) optional speciﬁcations (P ∗ , Pup low

∗ 2. Preprocessing step: execute step I only once at Tspec ∗ if not speciﬁed: 3. Generate P ∗ and Pup low If type 1 or 2 P ∗ − − z1 surface, then P ∗

low

= P+

mech

(0) and

∗ = 2P − Pup (1) mech ∗ = 2P If type 3 P ∗ − − z1 surface, then P ∗ = P + (0) and Pup p−c mech low ∗ = 3P ∗ If type 4 P ∗ − − z1 surface, then P ∗ =0.0044 and Pup low

low

97

Table 2 Construction of temperature–composition phase diagrams. 1. Read input speciﬁcations: (a) parameters for the AVDW EoS: ˛1,1 /˛2,2 and ˛1,2 /˛2,2 ∗ (b) reduced pressure: Pspec ∗ (default value: (c) temperature range: T ∗ (default value: 0.02829), Tup low 0.0943) (d) number of point calculations: l (default value: 500) 2. For i = 1 to l ∗ − T ∗ )/l T ∗ = T ∗ + i(Tup low low Solve steps I, II and III 3. Report output information: number and type of phases in equilibrium at all temperatures mole fraction of each component, xk , and density of each phase, k , at all i pressures

temperature range over which most ﬂuid phase phenomena are observed is between 30 % and 100 % of the critical temperature of the heavy component (component 2 for the algorithm). Based on the scaling applied to the AVDW EoS, the critical temperature of ∗ = 0.0943. The temperature component 2 in reduced units is Tc,2 ∗ , T ∗ ] = [0.3 × 0.0943, 0.0943] is provided as default, range [Tlow up although in many cases a large number of unnecessary calculations in single phase regions is likely to be performed. Furthermore, if the speciﬁed pressure is greater than the critical pressure of component 2, parts of the phase space which are of interest may not be explored; such cases involve ﬂuid–ﬂuid separations at temperatures higher than the critical temperature of component 2 exhibited by type III mixtures according to the classiﬁcation of Scott and van Konynenburg [4,12]. The default number of point calculations is set to 500 as before. The procedure followed for the construction of a temperature–composition phase diagram is summarized in Table 2. The default output information consists of a set of data ﬁles that contain the temperature and the corresponding equilibrium compositions and densities for each of the stable two-phase regions. 3. Results and discussion 3.1. P–T ﬂash calculations Single point calculations are the central function of the algorithm and form the building block for any other application, such as the construction of phase diagrams. The algorithm has been tested extensively over the entire range of parameter values that correspond to non-azeotropic mixtures for the AVDW EoS (˛1,1 = [0.45, 0.95] and ˛1,2 = [0.452, 0.992]) over a wide range of temperature and pressure conditions. The 69 parameter combinations for which the algorithm has been tested are provided as supplementary material. By varying the temperature and the pressure speciﬁcations, the total number of single point calculations performed is of the order of tens of thousands. The algorithm proved reliable in all cases, including liquid–liquid equilibrium, vapor–liquid equilibrium, three-phase (vapor–liquid–liquid) equilibrium, and near-critical calculations. Reliability was the main concern with the prototype implementation developed in GAMS. Nevertheless performance in terms of CPU time also proved to be reasonable. The average CPU time on an Intel Pentium III 600 MHz processor was 1.5 s and never exceeded 2 s. This implies that, regardless of the difﬁculty of the problem, the CPU time is relatively constant.

4. For i = 1 to l ∗ − P ∗ )/l P ∗ = P ∗ + i(Pup low low Solve steps II and IIIa

3.2. Construction of phase diagrams

5. Report output information: (a) number and type of phases in equilibrium at each pressure (b) mole fraction of each component, xk , and density of each phase, k , at i each pressure

Each of the ﬁve types of phase behavior in the classiﬁcation of Scott and van Konynenburg [4,12] that can be observed with the AVDW equation is brieﬂy reviewed for the case of non-azeotropic mixtures.

98

A. Giovanoglou et al. / Fluid Phase Equilibria 275 (2009) 95–104

Fig. 1. Pressure–temperature projection of the phase diagram for a type I mixture (˛1,1 /˛2,2 = 0.5 and ˛1,2 /˛2,2 = 0.8). The continuous curves denote the vapor pressure curves of components 1 and 2. The dashed curve denotes the vapor–liquid critical line, and the black circles the pure component critical points. The dotted lines indicate the temperature and pressure projections at which phase diagrams were constructed.

3.2.1. Type I phase behavior Binary mixtures of components of similar size and chemistry, such as methane + ethane, exhibit pressure–temperature phase diagrams corresponding to type I phase behavior in the classiﬁcation of Scott and van Konynenburg [4,12]. Here, a type I model mixture with parameters ˛1,1 /˛2,2 = 0.5 and ˛1,2 /˛2,2 = 0.8 is examined as a representative example. The pressure–temperature diagram for this mixture is shown in Fig. 1. The main characteristic of type I phase behavior is the existence of a single vapor–liquid critical line that connects the pure component critical points. Three isothermal pressure–composition phase diagrams at T1 = 0.037720, T2 = 0.061295 and T3 = 0.084870 indicated in the ﬁgure were obtained;

a total of 2100 successful point calculations were performed following the procedure outlined in Table 1. The implementation of the algorithm was always successful in locating the two-phase region through the use of the pure component vapor pressure test (see Step IIb in Part I of this work [27]), and in solving the corresponding equilibrium problem (even in the vicinity of the critical point). With many algorithms, calculations close to a vapor–liquid critical point, and in particular at pressures higher than the pseudocritical point, are likely to converge to the trivial solution, not only because the two phases have similar densities and compositions, but more importantly because both phases belong to the same density branch (pattern g, e.g., Fig. 9 in [27]). Such a situation is problematic even when cubic equations of state are used; because a single density root characterizes both the vapor and the liquid phases, there is no analytical way to distinguish them. In Fig. 2(a) the locus of the limits of material and mechanical stability at T2 are shown (similar plots are also presented by Rowlinson and Swinton [3]). Such a ﬁgure can be automatically generated since the limits of mechanical and material stability are calculated in step II of the algorithm developed in Part I [27]. At pressures between the pseudo-critical point (point E) and the critical point (point Pc ), material instability is the only driving force for phase separation. The corresponding limits of stability lie on the single density branch of pattern g and are used in the algorithm to bound the two phases during the solution of the restricted necessary conditions for equilibrium in step IIIa. In this way, convergence to the trivial solution is prevented. At pressures lower than the pseudo-critical point, the presence of mechanical instability is associated with the existence of multiple density branches in the density composition pattern. In Fig. 2(a) and (b), the lower branch of the material stability curve is used as a bound for the vapor phase, while the upper branch of the curve is used as a bound for the liquid phase. Furthermore, as shown in the density–composition plot of Fig. 2(b), the limits of material stability also provide tighter bounds on the equilibrium phases than the limits of mechanical stability at pressures below the pseudo-critical point. In Fig. 3 isobaric density–composition patterns obtained for six different pressure slices at temperature T2 of Fig. 1, and the

Fig. 2. (a) Limits of material and mechanical stability in the pressure–composition phase diagram at T2 of Fig. 1, and, (b) the projection of the same phase diagram in the density–composition space. The thick continuous curves denote equilibrium boundaries, the continuous curves of intermediate thickness the limits of material stability, and thin continuous curves the limits of mechanical stability. The dashed lines denote tie-lines that connect equilibrium points at pressures P1 , P2 and P3 , and points E and Pc correspond to the pseudo-critical and the critical point of the mixture, respectively.

A. Giovanoglou et al. / Fluid Phase Equilibria 275 (2009) 95–104

99

Fig. 3. Three thermodynamic spaces at six different pressures based on the phase diagram presented in Fig. 1. The pressure increases from top to bottom. The ﬁgures in the ﬁrst column correspond to density–composition patterns, the ﬁgures in the second column to Gibbs free energy–composition (G − z1 ) spaces, and the ﬁgures in the third column to the variation of chemical potential of component 2 vs that of component 1 (2 − 1 ). Mechanically and materially unstable regions are denoted with thin continuous curves, and black circles are used to highlight the vapor–liquid equilibrium points. The thick continuous curves denote regions that are materially stable or metastable.

corresponding dependence of the Gibbs free energy on composition (G − z1 ), and of the chemical potential of component 2 on the chemical potential of component 1 (2 − 1 ) are presented. In the ﬁrst row of the ﬁgure, the pressure is less than the pressure at the − (0), low-density limit of mechanical stability of component 2, Pmech which gives rise to a density–composition pattern of type b. In this case the liquid branch and the mechanically unstable branch are connected, while the vapor branch is detached. The corresponding G − z1 and 2 − 1 spaces consist of multiple branches. The liquid branch is connected to the mechanically unstable one at a point where the ﬁrst derivative is discontinuous, while the vapor branch is detached from the other two. The point where the Gibbs free energy is non-differentiable is the limit of mechanical stability. In the isobaric and isothermal G − z1 space, the equilibrium compositions share a common tangent, while in the 2 − 1 representation the vapor and the liquid branches cross at the equilibrium point in order to satisfy the equality of chemical potentials. In the second − row of Fig. 3, the pressure is between Pmech (0) and the pseudocritical pressure Pp−c . The density–composition pattern is of type d and the three branches are connected. In this case there are two limits of mechanical stability at which differentiability is not satisﬁed in the G − z1 and 2 − 1 spaces. In addition, the size of the unstable region has decreased. The third row of the ﬁgure corresponds to the pseudo-critical pressure. Here the density composition pattern changes from type d to type g. At this pressure, there is a oneto-one mapping between density and composition for all of the G − z1 space; here equilibrium is not due to the existence of multiple branches but rather to variations in the curvature of G. On the other hand, the 2 − 1 curve still exhibits a loop and two nondifferentiable points, the limits of material stability. The 2 − 1 curve crosses back onto itself at the equilibrium point. As the pressure is increased above the pseudo-critical point (fourth row in the ﬁgure) the density–composition pattern is still of type g, but the

region of material instability and the loop in the 2 − 1 space become smaller until they disappear at the critical pressure (row ﬁve in the ﬁgure). Above the critical pressure (last row) there is no phase separation and both material and mechanical stability are preserved for all compositions. 3.2.2. Type II phase behavior Binary mixtures of components of similar size but with weak unlike interactions, such as mixtures of alkanes and perﬂuoroalkanes of the same carbon number, often exhibit type II phase behavior [3]. Here, a type II model mixture with parameters ˛1,1 /˛2,2 = 0.5 and ˛1,2 /˛2,2 = 0.69 is studied. The pressure–temperature phase diagram for this mixture is shown in Fig. 4. Type II phase behavior is characterized by a demixing in the liquid region at sufﬁciently low temperatures. A three-phase line close to the vapor–pressure curve of component 1 that terminates at an upper critical end point (UCEP) is observed. In addition to the vapor–liquid critical line, a second critical line, related to liquid–liquid equilibrium, starts from the UCEP and extends to higher pressures. This critical line is usually referred to as the upper critical solution temperature (UCST) curve, since at higher temperatures liquid–liquid separation no longer exists. In the P–T projection shown in Fig. 4 the UCEP occurs at a relatively high pressure causing the UCST to intersect with the vapor–liquid critical line Three isothermal pressure–composition phase diagrams corresponding to temperatures T1 = 0.042435, T2 = 0.051865 and T3 = 0.070725 of Fig. 4 were obtained. These three temperatures were selected since the corresponding pressure–composition phase diagrams are representative of the behavior that type II mixtures can exhibit. To construct the phase diagrams a total of 2400 point calculations were performed successfully. Liquid–liquid equilibrium is detected by the existence of limits of material stability on the liquid density branch and on the single branch of pattern g. These

100

A. Giovanoglou et al. / Fluid Phase Equilibria 275 (2009) 95–104

Fig. 4. Pressure–temperature projection of the phase diagram for a type II mixture (˛1,1 /˛2,2 = 0.5 and ˛1,2 /˛2,2 = 0.69). The continuous curves denote the vapor pressure curves of components 1 and 2 and the three-phase line that terminates at an UCEP denoted by a black triangle. The dashed curves denote critical lines and the black circles the pure component critical points. The dotted lines indicate temperatures and pressures at which phase diagrams have been constructed to test the algorithm.

limits are then used to prevent convergence of the subsequent liquid–liquid equilibrium calculation to the trivial solution or to an unstable solution. Furthermore, the pure component vapor pressure test is very useful to detect the appearance of the vapor phase and to avoid unnecessary calculations in its absence, especially in the very small two-phase L2V regions seen at T1 and T2 . At constant temperature, over a range of pressures close to the three-phase pressure, the equilibrium problem has a large number of solutions. The existence of these solutions can be understood by studying the locus of the limits of mechanical and material stability on the phase diagram. This information is readily available from the algorithm and an indicative example is depicted in Fig. 5. It should be noted here that the metastable and unstable equilibrium points shown in the ﬁgure cannot be obtained from the algorithm, as these are discarded. They were obtained with the gPROMS software [31] by solving the necessary conditions for phase equilibrium from a suitable starting point (see [27]). The branch of the material stability curve associated with the liquid phases exhibits a maximum (point C1 ) and a minimum (point C2 ) in pressure and ﬁnally meets the other branch of the curve at the vapor–liquid critical point (C3 ). Points C1 and C2 are also critical points: C1 corresponds to a UCST, while C2 is an unstable critical point located inside the stable L1V two-phase region. The origin of point C2 and the reason why it is unstable are more clearly understood by examining the behavior of the (∂1 /∂x1 )T ∗ ,P ∗ − x1 curve along the liquid density branch at pressures close to that of point C2 (see Fig. 6). At a pressure P1 well below that of point C2 the curve exhibits a minimum and a maximum, both of which are negative. This indicates that there is no liquid–liquid separation. The single root of this curve corresponds to the limit of material stability on the liquid branch at that pressure, and is associated with the stable equilibrium between liquid phase 1 and the vapor phase. As the pressure increases, the curve is shifted upwards. At pressure P3 the maximum point is exactly on the zero axis and corresponds to point C2 . At pressures higher than P3 , the maximum point is strictly positive and the curve has three roots, which correspond to the limits of material stability. This indicates phase separation on the liquid density branch, i.e., liquid–liquid demixing. The shape of the curve at pressure P3 suggests that at

Fig. 5. The locus of the limits of material and mechanical stability in the pressure–composition phase diagram at T2 of Fig. 4. The thick continuous curves correspond to stable equilibrium boundaries, dash-dotted curves to metastable equilibrium states, and dotted curves to unstable equilibrium states. The thin continuous curves denote the limits of mechanical stability and the dashed curves the limits of material stability. The pseudo-critical point is denoted by E, and C1 , C2 and C3 denote critical points.

point C2 the following conditions must hold:

∂1 ∂x1

=0,

∂2 G ∂x12

=0, T ∗ ,P ∗

∂x12

T ∗ ,P ∗

or equivalently:

∂2 1

∂3 G ∂x13

=0

and

∂x13

T ∗ ,P ∗

=0 T ∗ ,P ∗

and

∂3 1

∂4 G ∂x14

< 0, T ∗ ,P ∗

(1)

< 0.

(2)

T ∗ ,P ∗

The fact that the ﬁrst and the second derivatives of the chemical potential (or the second and third derivatives of the Gibbs free energy) are zero implies that C2 is a critical point, but the third derivative being negative implies that it is an unstable critical point [32]. An unstable liquid–liquid equilibrium region, as shown in Fig. 5, appears at point C2 and extends to slightly higher pressures.

Fig. 6. Form of the (∂1 /∂x1 )T ∗ ,P ∗ –x1 curves along the liquid density branch at four pressures close to point C2 of Fig. 5: P1 = 0.0016, P2 = 0.0018, P3 = 0.0020428 and P4 = 0.0022.

A. Giovanoglou et al. / Fluid Phase Equilibria 275 (2009) 95–104

101

Fig. 7. (a) The density–composition pattern, and (b) the 2 − 1 space corresponding to a reduced pressure P ∗ = 0.0022 in Fig. 5. In (b) the multiplicity in the solution space at these conditions is shown. The thick continuous curves denote stable and metastable states, the dashed curves denote unstable states, and the black circles denote equilibrium points. S denotes the stable solution, M the metastable solution, and U1 and U2 correspond to unstable solutions.

Both phases in this region of phase space are unstable since (∂1 /∂x1 )T ∗ ,P ∗ is negative along the equilibrium curves. There is an additional unstable two-phase region connecting the metastable and unstable liquid–liquid equilibrium regions originating from points C1 and C2 , respectively. In this case the phase which occurs at low compositions is metastable, (∂1 /∂x1 )T ∗ ,P ∗ > 0, while the one at high compositions is unstable, (∂1 /∂x1 )T ∗ ,P ∗ < 0. The existence of these unstable and metastable regions close to the three-phase line increases the difﬁculty of solving the phase equilibrium problem in this region with standard techniques. This is illustrated in Fig. 7, where the density–composition pattern and the 2 − 1 curve for a pressure slightly below the three-phase line of the phase diagram in Fig. 5 are shown. The density composition pattern is of type d. At this pressure, stable equilibrium occurs between liquid phase 1 and the vapor phase. In the 2 − 1 space, the solution corresponds to point S, which is located at the intersection of the liquid and vapor branches. The ﬁgure also reveals that at these conditions the problem has three additional solutions, in addition to the trivial solution. Two of them, U1 and U2, are unstable while the other, M, is metastable. These solutions are located on the liquid branch and occur due to the existence of multiple regions of instability. Solution U1 is located at the intersection of two unstable branches, solution U2 is located at the intersection of a metastable and an unstable branch, and ﬁnally solution M is located at the intersection of two metastable branches. The proposed algorithm converges to the correct solution S without any problem. Since all of the unstable regions of the density composition pattern are discarded in step II, both unstable solutions U1 and U2 would never occur. Only solutions M and S are obtained in step IIIa; of these, S is found to be stable. 3.2.3. Type III phase behavior The extent of the liquid–liquid region observed in mixtures of type II increases as the size difference of the two components, or their chemistry, further differs. The vapor–liquid critical line merges with the loci of UCSTs and the result is a critical line that originates from the critical point of the less volatile component and extends up to high pressures. Water + n-alkanes are typical binary mixtures exhibiting type III phase behavior [3]. A second critical line starts from the critical point of the more volatile component and extends up to a UCEP. Here, a type III model mixture with parameters ˛1,1 /˛2,2 = 0.5 and ˛1,2 /˛2,2 = 0.62 is examined. The pressure–temperature projection of the phase diagram for this mixture is shown in Fig. 8. Two isothermal pressure–composition phase diagrams at T1 = 0.048093 and T2 = 0.066010 (indicated on

the ﬁgure) where calculated for this mixture. To construct the phase diagrams, 1200 point calculations were required in total, all of which were carried out successfully. In Fig. 9 the locus of the limits of material and mechanical stability as well as the metastable and unstable equilibrium regions corresponding to the phase diagram at T1 of Fig. 8 are shown. The limits of stability are readily available from the algorithm, while the metastable and unstable equilibrium points were obtained using gPROMS [31] as explained previously. Unlike the phase diagram of Fig. 5, there are only two critical points in this case. This is because the liquid–liquid region is not bounded by a critical point at high pressure. Critical point C2 is unstable and is similar to that occurring in the type II phase diagram of Fig. 5. Point C3 is a stable vapor–liquid critical point. Another important characteristic of the phase diagram shown in Fig. 9 is that, in contrast to Fig. 5, the L1L2 and L1V phase regions are connected via metastable and unstable equilibrium lines. The shape of these lines is more clearly demonstrated in the insets shown in Fig. 10 (a and b). At low compositions, Fig. 10(a), the L1L2 and L1V stable equilibrium curves are connected due to the existence of two peaks. Close to pure component 1, Fig. 10(b),

Fig. 8. Pressure–temperature projection of the phase diagram for a type III mixture (˛1,1 /˛2,2 = 0.5 and ˛1,2 /˛2,2 = 0.62). The continuous curves denote the vapor pressure curves of components 1 and 2 and the three-phase line which terminates at an UCEP denoted by a black triangle. The dashed curves represent critical lines and the black circles denote pure component critical points. Finally, the dotted lines denote the temperature projections at which phase diagrams have been constructed.

102

A. Giovanoglou et al. / Fluid Phase Equilibria 275 (2009) 95–104

Fig. 9. The location of the limits of material and mechanical stability in the pressure–composition phase diagram at T1 = 0.048093 of Fig. 8. The thick continuous curves correspond to stable equilibrium states, the dash-dotted curves to metastable equilibrium states, and the dotted curves to unstable equilibrium states. The thin continuous curves represent the limits of mechanical stability and dashed curves the limits of material stability. C2 and C3 denote critical points.

Fig. 11. Pressure–temperature projection of the phase diagram for a type IV mixture (˛1,1 /˛2,2 = 0.43 and ˛1,2 /˛2,2 = 0.707852). The continuous curves correspond to the vapor pressure curves of components 1 and 2 and the three-phase lines which terminate at UCEP and LCEP denoted by black triangles. The dashed curves are critical lines and the black circles denote pure component critical points. The dotted line denotes the pressure at which the phase diagram in Fig. 12 was constructed.

the stable equilibrium lines are connected via an “S” shaped curve. Furthermore, a loop is formed between critical points C2 and C3 . The largest part of that loop is unstable or metastable, while only the L2V phase region, above the three-phase line, is stable. As the temperature is increased the entire loop shifts inside the L1V region while the two critical points come closer to each other. The loop ﬁnally disappears when the two critical points coincide.

found close to the critical point of component 1 and is bounded by an LCEP and a second UCEP. This results in three distinct stable critical lines as shown in the ﬁgure. The ﬁrst one emerges from the UCEP at low temperatures and extends up to very high pressures. The second critical line starts from the critical point of component 2 and terminates at the LCEP. A third critical line is located between the critical point of component 1 and the second UCEP. This topology of the critical lines gives rise to the existence of miscibility windows in the temperature–composition slices at constant pressure. An indicative example is shown in the temperature–composition phase diagram of Fig. 12. At low temperatures, the mixture separates into two liquid phases and remains immiscible up to an UCST. At much higher temperatures, the mixture exhibits a small region of liquid–liquid separation (as is more clearly seen in the inset) and two vapor–liquid equilibrium regions. Between the low- and hightemperature phase separation regions the mixture is totally miscible. This type of behavior is very relevant in solvent design studies where the aim is to determine the temperature and pressure window and the amount of solvent to be used for an effective separation or reaction within a process. Taking advantage, for instance, of the

3.2.4. Type IV phase behavior Size asymmetry, accompanied by differing chemistry of the components, leads to mixtures that exhibit phase behavior of type IV. Methane + hex-1-ene and a number of polymer–solvent mixtures are typical examples of systems exhibiting type IV phase behavior [3]. Here, a type IV model mixture with parameters ˛1,1 /˛2,2 = 0.43 and ˛1,2 /˛2,2 = 0.707852 is studied. The pressure–temperature projection of the phase surface for this mixture is shown in Fig. 11. Type IV phase behavior combines features of both type II and type V (type V will be discussed in the next section). Two separate three-phase lines are observed in this case: the ﬁrst occurs at low temperatures and ends at an UCEP, while the second is

Fig. 10. Two insets from Fig. 9.

A. Giovanoglou et al. / Fluid Phase Equilibria 275 (2009) 95–104

Fig. 12. Isobaric temperature–composition phase diagram for a binary mixture with parameters ˛1,1 /˛2,2 = 0.43 and ˛1,2 /˛2,2 = 0.707852 at the reduced pressure slice of Fig. 11: P1 = 0.00242546. The black circles denote critical points.

characteristics of type IV phase behavior, the same solvent could be used to guarantee a single homogeneous liquid phase within a reaction unit and liquid–liquid phase separation inside a subsequent separation process. As an interesting example of this, Yelash and Kraska [2] have obtained the pressure–temperature diagrams for mixtures of small spherical molecules and chain molecules and have represented solvent and anti-solvent processes in polymerization as paths in the phase diagram. To construct the phase diagram shown in Fig. 12, 700 P–T ﬂash calculations were required. 3.2.5. Type V phase behavior Methane + n-hexane is a typical mixture exhibiting type V phase behavior [12], and in general mixtures of methane with n-alkanes of higher molecular weight than n-hexane also exhibit this type of phase behavior. The mixture phase behavior is characterized by the presence of a region of liquid–liquid immiscibility close to the crit-

103

Fig. 14. Isothermal pressure–composition phase diagram for a binary mixture with parameters ˛1,1 /˛2,2 = 0.38 and ˛1,2 /˛2,2 = 0.759 at the reduced temperature slice of Fig. 13, T1 = 0.036490.

ical point of the light component. Studies of the transition between types I and V can be formulated in terms of process and solvent optimisation in polymerization processes [1]. Here, a model mixture with parameters ˛1,1 /˛2,2 = 0.38 and ˛1,2 /˛2,2 = 0.759 is studied. The pressure–temperature projection of the phase surface for this mixture is shown in Fig. 13. As mentioned, in this case the mixture exhibits a small region of liquid–liquid separation close to the critical point of component 1. A three-phase line is also observed and is bounded by upper and lower critical end points as seen more clearly in the inset of the ﬁgure. A stable critical line extends from the critical point of component 2 to the LCEP while a second much smaller critical line starts at the critical point of component 1 and extends towards the UCEP. In Fig. 14 the pressure–composition phase diagram corresponding to a temperature T1 which lies between the lower and the upper critical end points of Fig. 13 is shown. The phase diagram consists of a large two-phase region where a liquid phase is in equilibrium with the vapor phase (L1V), a liquid–liquid equilibrium (L1L2) region which is bounded at high pressure by a critical point, and a second vapor–liquid equilibrium region (L2V) also limited by a critical point (better seen in the inset). In total, 600 point calculations were required to construct the phase diagram of Fig. 14. The use of the material stability tests allows the L2V phase region to be detected, while the use of the limits of material stability as bounds guarantees the successful solution of the phase equilibrium calculation. 4. Conclusions

Fig. 13. Pressure–temperature projection of the diagram for a type V mixture (˛1,1 /˛2,2 = 0.38 and ˛1,2 /˛2,2 = 0.759). The continuous curves correspond to the vapor pressure curves of components 1 and 2 and the three-phase line. The threephase line lies between an LCEP and an UCEP denoted by black triangles. The dashed curves denote critical lines and the black circles pure component critical points. The dotted line denotes the temperature projection at which the phase diagram in Fig. 14 was constructed.

A prototype implementation of the algorithm proposed in Part I of this work [27] for P–T ﬂash calculations and for the automated construction of isothermal pressure–composition and isobaric temperature–composition phase diagrams is presented. The performance of the algorithm is thoroughly tested on a large number of P–T ﬂash calculations for non-azeotropic mixtures. The algorithm proves to be reliable despite the complexity of the phase behavior studied. The mechanical and material stability tests employed in step II of the algorithm are successful in identifying the number and type of existing phases for difﬁcult points encountered at near-critical behavior, when liquid–liquid separation is observed or close to three-phase lines. In step III, the solution of the two-phase equilibrium between all pairs of existing phases is obtained reliably by applying the restricted

104

A. Giovanoglou et al. / Fluid Phase Equilibria 275 (2009) 95–104

necessary conditions for equilibrium. Liquid–liquid equilibrium regions and small two-phase regions close to critical points and three-phase lines are probably the most challenging cases when binary non-azeotropic mixtures are studied. Even in such cases the algorithm proves very reliable. The average CPU time on an Intel Pentium III 600MHz processor is 1.5 s and never exceeds 2 s. Given the prototypical nature of the implementation in GAMS, this consistency in the speed of calculations is promising for a high-performance implementation of the approach. Results on the automated construction of the pressure– composition and temperature–composition phase diagrams is presented through a review and analysis of the different types of phase behavior seen in binary mixtures exhibiting non-azeotropic ﬂuid phase behavior in the classiﬁcation of Scott and van Konynenburg [4,12]. The augmented van der Waals equation of state is used to study the ﬁve main types of ﬂuid phase behavior characteristic of these mixtures (types I–V). Pressure–compositions and temperature–compositions phase diagrams are constructed to highlight the main features of each non-azeotropic type. The loci of the limits of mechanical and material stability (readily calculated with the algorithm) are also shown in order to demonstrate the complexity of the problem and to analyze the performance of the algorithm in difﬁcult regions of phase space. Positive azeotropic behavior is associated with types II and III, giving rise to types II-A and III-HA. This type of behavior will be discussed in a forthcoming publication. It is also the intention to extend the phase equilibrium algorithm to more sophisticated equations of state such as the statistical associating ﬂuid theory (SAFT) [33,34] and SAFT-VR [35–38]. Acknowledgments We thank D.N. Theodorou for very useful insights during the PhD examination of A. Giovanoglou. Funding to the Molecular Systems Engineering group from the Engineering and Physical Sciences Research Council (EPSRC) of the UK (grants GR/N35991 and EP/E016340), the Joint Research Equipment Initiative (JREI) (GR/M94427), and the Royal Society-Wolfson Foundation refurbishment grant is gratefully acknowledged. Appendix A. Supplementary Data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.ﬂuid.2008.08.018.

References [1] A. Giovanoglou, C.S. Adjiman, A. Galindo, G. Jackson, Edited by A. Kraslawski, I. Turunen, Proceedings of the 13th European Symposium on Computer Aided Process Engineering, Computer-Aided Chemical Engineering, vol. 14, Elsevier, 2003, pp. 137–142. [2] L.V. Yelash, T. Kraska, Phys. Chem. Chem. Phys. 1 (1999) 4315. [3] J.S. Rowlinson, F.L. Swinton, Liquids and Liquid Mixtures, 3rd ed., Butterworths, 1982. [4] R.L. Scott, P.H. van Konynenburg, Discuss. Faraday Soc. 49 (1970) 87. [5] D. Furman, R.B. Grifﬁths, Phys. Rev. A 17 (1978) 1139. [6] J.A. Barker, W. Fock, Discuss. Faraday Soc. 15 (1953) 188. [7] G. Jackson, Mol. Phys. 72 (1991) 1365. [8] L.A. Davies, G. Jackson, L.F. Rull, Phys. Rev. Lett. 82 (1999) 5285. [9] L.A. Davies, G. Jackson, L.F. Rull, Phys. Rev. E 61 (2000) 2245. [10] L.V. Yelash, T. Kraska, Ber. Bunsen-Ges. Phys. Chem. 102 (1998) 213. [11] R.L. Scott, Phys. Chem. Chem. Phys. 1 (1999) 4225. [12] P.H. van Konynenburg, R.L. Scott, Philos. Trans. A 298 (1980) 495. [13] V.A. Mazur, L.Z. Boshkov, V.G. Murakhovsky, Phys. Lett. A 104 (1984) 415. [14] L.Z. Boshkov, Dokl Akad. Nauk SSSR 294 (1987) 901. [15] U.K. Deiters, I.L. Pegg, J. Chem. Phys. 90 (1989) 6632. [16] A. van Pelt, C.J. Peters, J. de Swaan Arons, J. Chem. Phys. 95 (1991) 7569. [17] T. Kraska, U.K. Deiters, J. Chem. Phys. 96 (1992) 539. [18] A. Bolz, U.K. Deiters, C.J. Peters, T.W. de Loos, Pure Appl. Chem. 70 (1998) 2233. [19] B.J. Alder, C.E. Hecht, J. Chem. Phys. 50 (1969) 2032. [20] K.N. Marsh, M.L. McGlashan, C. Warr, Trans. Faraday Soc. 66 (1970) 2453. [21] N.F. Carnahan, K.E. Starling, AIChE J. 18 (1972) 1184. [22] G. Jackson, J.S. Rowlinson, C.A. Leng, J. Chem. Soc.: Faraday Trans. 1, 82 (1986) 3461. [23] N.F. Carnahan, K.E. Starling, J. Chem. Phys. 51 (1969) 635. ´ [24] J. Kolafa, I. Nezbeda, J. Pavlcek, W.R. Smith, Fluid Phase Equilib. 146 (1998) 103. ´ [25] J. Kolafa, I. Nezbeda, J. Pavlcek, W.R. Smith, Phys. Chem. Chem. Phys. 1 (1999) 4233. [26] J.-L. Wang, R.J. Sadus, Fluid Phase Equilib. 214 (2003) 67. [27] A. Giovanoglou, A. Galindo, G. Jackson, C.S. Adjiman, Fluid Phase Equilib. 275 (2009) 79. [28] A. Brooke, D. Kendrick, A. Meeraus, R. Raman, GAMS Language Guide, GAMS Development Corporation, Washington, DC, USA, http://www.gams.com, 1997. [29] A.R. Imre, T. Kraska, J. Chem. Phys. 122 (2005) 064507. [30] M. Cismondi, M.L. Michelsen, J. Supercrit. Fluids 39 (2007) 287. [31] PSE Ltd. gPROMS Introductory User Guide, London, UK, http://www. psenterprise.com, 2004. [32] J.W. Tester, M. Modell, Thermodynamics and Its Applications, 3rd ed., Prentice Hall, 1997. [33] W.G. Chapman, K.E. Gubbins, G. Jackson, M. Radosz, Fluid Phase Equilib. 52 (1989) 31. [34] W.G. Chapman, K.E. Gubbins, G. Jackson, M. Radosz, Ind. Eng. Chem. Res. 29 (1990) 1709. [35] A. Gil-Villegas, A. Galindo, P.J. Whitehead, S.J. Mills, G. Jackson, A.N. Burgess, J. Chem. Phys. 106 (1997) 4168. [36] A. Galindo, L.A. Davies, A. Gil-Villegas, G. Jackson, Mol. Phys. 93 (1998) 241. [37] L.A. Davies, A. Gil-Villegas, G. Jackson, Int. J. Thermophys. 19 (1998) 675. [38] L.A. Davies, A. Gil-Villegas, G. Jackson, J. Chem. Phys. 111 (1999) 8659.