Calculations of Complex Phase Equilibrium by ... - ACS Publications

Article pubs.acs.org/IECR

Calculations of Complex Phase Equilibrium by Semigrand Canonical Ensemble Yiping Tang* Honeywell Process Solutions, London, ON N6A 6K2, Canada

Ensheng Zhao Honeywell Process Solutions, Calgary, AB T2H 0C2, Canada ABSTRACT: The recently developed method of semigrand canonical (sGC) ensemble (Tang, Y. J. Chem. Phys. 2012, 136, 034505) was employed to calculate phase equilibria for binary mixtures. The calculations cover vapor−liquid equilibria (VLE), liquid−liquid equilibria (LLE), solid−liquid equilibria (SLE), and solid−solid equilibria (SSE). The thermodynamic models in this study include cubic and noncubic equations of state, the activity coefficient model, and the alloy-based model. On the basis of extensive calculations and comparisons with classical methods, we conclude that the sGC method is universal for calculating complex phase equilibria. In particular, it is very robust in cases of multiple vdW loops with and without intersections, which are often difficult and require special treatments for classical methods. The phase stability characteristics of the sGC method are also addressed.

1. INTRODUCTION Phase equilibrium calculations are essential in the oil and gas industry for the design and simulation of distillation columns, pressure vessels, heat exchangers, and so on. A successful simulation depends on both the rationality of the employed models and the reliable and efficient computation algorithm. There are many models for describing phase behavior, and these models can be classified into cubic equations of state (EOS), noncubic EOS, and activity coefficient models (ACMs). Representative examples are the Peng−Robinson (PR) EOS, the perturbed-chain statistical associating fluid theory (PCSAFT), and the van Laar ACM. In some cases, one must model solid solutions, and the model is, in general, highly empirical. For computation algorithms, one must address the issues of both phase stability and phase equilibrium. Decades ago, Michelson1 successfully developed the minimum tangent plane distance (TPD) criterion for phase stability and the successive substitution (SS) method for phase equilibrium calculation.2 The SS method or, more generally, equation-solving methods make use of the equality of the chemical potentials of each component in different phases. The method requires an initialization and the updating of fugacity coefficients in the computation process. Another commonly mentioned approach is the equal-area (EA) method, which makes use of equal areas between the van der Waals (vdW) loop above and below the equilibrium line3−6 to extract equilibrium points, noting that the vdW loop exists in all analytical equations of state and even in computer simulations at low temperatures. Both methods have advantages and disadvantages. The SS method is computationally very efficient and quite popular in industrial simulations. However, a trivial solution could occur and proper initialization is needed. In addition, the physics of the results is not guaranteed, which often requires an analysis of the global minimization of the Gibbs free energy of the system. The EA © 2013 American Chemical Society

method has no problem with nonconvergence. However, the justification of its solution is subject to the complexity of the vdW loops because these loops could be multiple and overlap each other. Therefore, neither of the two methods can warrant their solution with full confidence. In general, computational difficulties arise from the complexity of the phase behavior at the transition, such as from vapor−liquid equilibrium (VLE) to vapor−liquid−liquid equilibrium (VLLE) and liquid−liquid equilibrium (LLE) at some temperatures and pressures. To solve the problems, others have proposed many new methods.7−11 These solutions were developed by some mathematical manipulations that depend on the system and are not expected to be generally applicable. Recently, we proposed a theory to explore phase transitions within the framework of the grand canonical (GC) ensemble for pure fluids12 and the semigrand canonical (sGC) ensemble for binary mixtures.13 We found that, in principle, conventional intensive thermodynamic properties, namely, the Helmholtz free energy density for pure fluids and the Gibbs free energy for mixtures, are related to extensive properties, namely, the volume and number of molecules in the system, respectively. To explore their influence on thermodynamic properties, we merged two subsystems into a new subsystem and found that the corresponding free energy is changed or transformed at unstable states. By carrying out this transform successively, we eventually obtained the stable free energy. In the curve of stable Helmholtz free energy density or Gibbs free energy, there are a number of tangent lines whose end points are physically at phase equilibrium. The GC and sGC transforms show all Received: Revised: Accepted: Published: 9690

March 26, 2013 May 23, 2013 June 7, 2013 June 7, 2013 dx.doi.org/10.1021/ie400958v | Ind. Eng. Chem. Res. 2013, 52, 9690−9697

Industrial & Engineering Chemistry Research

Article

1 ln Ni + 1 ⎧ ⎫ ⎡ g (x + y ) + g (x − y ) ⎤ cx ⎪ ⎪ i ⎢ i ⎥Ni + 1⎬ dy − β exp⎨ ⎪ ⎪ ⎥⎦ 0 2 ⎩ ⎢⎣ ⎭

known characteristics of classical thermodynamics. They yield the same stability conditions as required by the TPD method (see the Appendix) and reproduce the classical relations of equal chemical potentials between two equilibrium phases and equal area to split the vdW loop at phase equilibrium. Despite these well-known relations, the corresponding SS and EA methods are only byproducts of the GC or sGC transform, which are necessary but insufficient conditions for phase equilibrium calculations. The GC and sGC methods have been well tested for Lennard-Jones fluids and mixtures modeled by the first-order mean spherical approximation (FMSA). The purpose of this work was to apply the sGC transform to phase equilibrium calculations for realistic mixtures and to compare its performance with that of the two classical methods. To demonstrate its universal capability, we chose several highly nonideal mixtures for testing, of which some are well discussed in the literature. Because the same philosophy was used in developing the GC and sGC methods, we believe that the conclusions drawn in this work are applicable to the GC method for pure fluids.

βgi + 1(x) = −

∫

(3)

in which cx = min(x, 1 − x) and Ni+1 is the number of molecules at stage i + 1 (Ni+1 = 2Ni). At the limit Ni → ∞, eq 3 reduces to the sGC equation ⎡ g (x + y ) + g (x − y ) ⎤ i ⎥ gi + 1(x) = min ⎢ i ⎥⎦ 0 ≤ y ≤ cx⎢ 2 ⎣

(4)

In eq 4, the transform of gi(x) is performed in a straightforward manner by finding the minimum of gi(x) with composition fluctuations in the range of [0, min(x, 1 − x)], noting that the fluctuations are of equal amplitude y around x. As an illustration, we plot in Figure 1 the profiles of gi(x) and in

2. SUMMARY OF THE SGC METHOD The sGC ensemble was originally developed to simulate the physical properties and phase behavior of molecules at the microscopic level, where temperature T; pressure P; and the difference between two chemical potentials Δμ21 = μ1 − μ2, with μ1 for component 1 and μ2 for component 2, are fixed.14,15 We reformulated the ensemble as double integrals from the canonical ensemble. Taking the chemical potential from sGC partition function and the Helmohltz free energy from the GC partition function, we obtained the following relation between μ2 and Δμ12 βμ2 = −

1 ln( N

∫0

1

exp{ − Nβ[g (x) − Δμ21x]} dx)

(1) Figure 1. sGC transform of the Gibbs free energy for Margules model with reduced parameters A12/RT = 2.2 and A21/RT = 1.0. The solid curve, three dashed curves, and solid line are for MFT, three transform stages, and the solution at stage 10, respectively.

where the reduced Gibbs free energy is given by β g (x ) = −

1 ln[ N

∫v

∞

exp( − Nβ {[Pv + a(x , v , T )]}) dv]

min

(2)

In eqs 1 and 2, x is the composition of component 1; N is the total number of molecules in the system; and v and vmin are the unit of volume and its minimum value, respectively. In addition to the composition, g(x) depends implicitly on temperature T and pressure P, which are omitted here for simplicity. Equation 2 can be taken as a new definition of Gibbs free energy for mixtures. This equation has been shown13 to exhibit a number of unusual thermodynamic behaviors. The most significant is that the intensive property g(x) is a function of the total number of molecules in the system N, which is an extensive property. To explore the effects of this equation on thermodynamics, a process was designed for g(x) to evolve for the system from small but macroscopic subsystems. We assumed that there are many identical subsystems and combined them successively in pairs. We found that each combination generates a new or transformed g(x) within the framework of the sGC ensemble. At the limit of the transform, g(x) is stable and exhibits a phase transition. Details about the combination and the corresponding transform can be found in ref 13. The conclusion is that, for the transform from stage i to stage i + 1, we have

Figure 2. Corresponding derivative of the Gibbs free energy for Figure 1. 9691

dx.doi.org/10.1021/ie400958v | Ind. Eng. Chem. Res. 2013, 52, 9690−9697


Article

Table 1. Comparsion between the SS Method of Heidemann and Michelsen7 and Our Implementation of the SS Algorithm P (MPa) 0.0367 (200 mmHg) 0.0213 (160 mmHg) 0.0207 (155 mmHg) ≤0.0153 (≤115 mmHg)

Heidemann and Michelsen7

this work

dozens of iterations hundreds of iterations oscillates among a few sets of phase compositions switches between full vapor and full liquid and eventually becomes chaotic

Figure 2 the corresponding profiles of gi′(x), defined as the derivative of gi(x) with respect to x, in the sGC transform. The two figures show the evolutions of gi(x) and g′i (x) with the successive transforms. gi(x) starts from the orginal form of gi(x) at stage 0, or g0(x), usually referred to as the Gibbs free energy of mean field theory (MFT), to that at stage 1, from that at stage 1 to that at stage 2, and so on until stage 10, where g10(x) is almost stationary. In the process of transformation, one can observe a substantional increment of oscillation frequency and a decrement of oscillation amplitude, resulting in two straight lines in Figures 1 and 2. It is noted that g′(x) is much more sensitive than g(x) to the transform, and subsequent discussions focus on g′(x) for illustration. In addition, eq 4 is also capable of checking phase stability, which has not previously been mentioned. The condition of phase stability was found to be identical to that from the TPD method, and its details are left to the Appendix.

similar level of iterations 225/213 iterations for inputs z = 0.15/0.8 at convergence nonconvergent even with a maximum of 1000 iterations behaves similarly; chaotic calculation continues even below the azeotrope point (0.0148 MPa)

mixtures, which are irrelevant to our calculation algorithm, can be found in the original works and are simply mentioned here. In our calculations, the grid of composition space was finely divided into 1000 intervals over the composition range [0, 1] for both the sGC and EA methods, and the tolerance of convergence was 1.0 × 10−12 for every point. For flash calculations by the SS method, we used 100 composition inputs within composition range [0, 1] and a tolerance of convergence of 1.0 × 10−10.

4. EXAMPLES OF PHASE EQUILIBRIUM CALCULATIONS 4.1. Pyridine + Acetic Acid. In applying the SS method to phase equilibrium calculation, Heidemann and Michelsen7 found that the convergence is highly dependent on the nonideality of the system. For binary mixtures, the convergence of the method is determined mathematically by the eigenvalues of the derivative matrix of the fugacity for an EOS or of the activity coefficient for an ACM. The absolute value of each eigenvalue should be less than unity for convergence. To verify this, Heidemann and Michelsen7 performed a VLE calculation at T = 353.15 K by the SS method for pyridine + acetic acid. The mixture was described by the ideal gas equation for the vapor phase and by the van Laar model for the liquid phase, showing an azeotropic phase envelope. They found that the convergence condition was violated in the composition range [0.18, 0.70 ] by the van Laar model, which severely downgraded the performance of the SS method. For the purpose of validating our SS method, our results, obtained using the flash computation described earlier, are summarized and compared with theirs in Table 1. In Figure 3, we show the results in detail, noting particularly the horizontal line for nonconvergence. This line concides with that found by Heidemann and Michelsen,7 and all calculations below it became chaotic. With these comparsions, we are confident that the SS method implemented in this work behaves very similarly to the literature SS method and that subsequent calculations in later cases are reliable. In contrast to the unstable performance of the SS method, the sGC method, the results of which are plotted as solid curves in Figure 3, is much different. It passes from higher pressure with a single phase envelope, through middle pressure with double phase envelopes, and to lower pressue with a single phase without any problems. In paticular, there are no symptoms of any calculation issues around the nonconvergent line, whereas the line is critical for the SS method. As one example, Figure 4 shows the calculation of the sGC transform at P = 0.018 MPa, where two VL equilibria coexist. To illustate the process of the sGC transform, we show two dashed lines at stage 3 and the final solid lines at stage 29. These lines were obtained successively by eq 4, without any prior considerations of the existence of vdW loops. The EA method, althought providing the same answer as the sGC method, has to check the number of vdW loops in advance and solves two equations of equal

3. NUMERICAL IMPLEMENTATION OF THE SGC METHOD Because eq 4 requires only Gibbs free energy, a basic property for any MFT model, its implementation is very straightforward. In this work, we apply the following procedure: (1) Take g0(x) from the given model and make a number of evenly distributed discrete points in the composition range [0, 1]. (2) Apply eq 4 to each point and obtain g1(x). (3) If g1(x) = g0(x) for all points, then the mixture is stable, and no phase transition occurs. In this case, the calculation can be terminated. (4) Otherwise, continue the process from stage i to stage i + 1. If phase stability is the only concern, a point x with gi+1(x) ≠ gi(x) indicates that phase stability has been lost. (5) If |gi+1(x) − gi(x)| meets the criterion of convergence, where we define gi+1(x) as g∞(x) for convenience, then the process can be terminated. (6) Determine equilibrium points by finding two end points of a linear tangent line in g∞(x). Note that there could exist multiple tangent lines and each of them represents one phase equilibrium. As a comparison, we also performed the calculations by the SS and EA methods. The two methods require chemical potentials or the derivatives of g(x) and a proper initialization. For the SS method, we followed the algorithm of Michelsen and Mollerup,16 where the system is initialized by Wilson’s method,7 and performed flash calculations to obtain phase equilibrium from a number of composition inputs. To clarify the true behavior of the SS method, neither analysis of the phase stability nor minimization of the Gibbs free energy was applied during the calculation. For the EA method, we followed the procedures of Tang et al.6 and assumed that each vdW loop represents one phase equilibrium. In the following discussion, the details of the cited models and the corresponding parameters and performances for our 9692



Article

is about the hydrocarbon side as shown in Figure 5, noting that the inset provides a broader view of the phase diagram. The

Figure 3. Phase diagram calculated by the SS (dotted curve), EA (dashed curve), and sGC (solid curve) methods for pyridine + acetic acid at T = 353.15 K, as modeled by the ideal gas equation for the vapor phase and by the van Laar ACM for the liquid phase. The dotted horizontal line is the convergence limit for the SS method.

Figure 5. Phase diagram calculated by the SS (dotted curve), EA (dashed curve), and sGC (solid curve) methods for n-butane + water at T = 344.26 K, as modeled by the Peng−Robinson EOS. The inset is a broader view of the diagram.

figure also shows the VLLE line, which was found to be at P = 0.864 MPa, very close to the value of P = 0.863 MPa reported by Hanif et al.,4 as well as the experimental result of P = 0.865 MPa. Our test for the SS method, running from low to high pressures, revealed that the method worked well for VLE1 but missed VLE2, likely because of the K-value initialization in Wilson’s method.7 More disappointingly, VLE1 spread to higher pressures until P = 0.913 MPa for the appearance of LLE. The true VLLE line was missed, and this problem was also found by Hanif et al.4 Compared with the earlier solution for pyridine + acetic acid, the SS calculation was numerically stable but yielded a partially unphysical answer for the present mixture. When applying the EA method for the solution, we assumed that each vdW loop corresponds to one phase equilibrium. Such an assumption was used to calculate the azeotropes in Figure 4 and found to work successfully. For the present system, we found that the EA method worked well for VLE1 and VLE2. However, it missed the VLLE line and failed to calculate LLE at high pressures. The phase curve on the right side of Figure 5 seems to deviate from the physical solution to a greater extent with increasing pressure. To understand this unusual behavior, we plot g′(x) in Figure 6 at P = 0.90 MPa, as well as the two additional values of P = 0.80 and 0.85 MPa. It can be seen that, at P = 0.80 MPa, there is only one vdW loop, for which the horizontal line for phase equilibrium can be solved from the EA relation A = B. For P = 0.85 MPa, two vdW loops can be seen in Figure 6b, and the two corresponding horizontal lines can be obtained by solving the two equations A = B and C = D. However, for the higher pressure of 0.90 MPa, the EA method was found to be problematic. Two vdW loops, and thus two phase equilibria, remain, as shown in inset c by two dashed horizontal lines. However, the two sets of compositions at equilibrium overlap each other, which is physically impossible. To solve this problem, Hanif et al.11 developed an algorithm by assuming that only one horizontal line, or a single phase equilibrium, exists and that the line can be located by solving the equation A + C = B + D, as shown in inset d of Figure 6. This algorithm was successful in calculating

Figure 4. sGC transform for pyridine + acetic acid at T = 343.15 K and P = 0.018 MPa, where the solid curves, dashed curves, and solid horizontal lines are for MFT and transforms at stages 3 and 10, respectively.

areas A = B and C = D for the two vdW loops individually. When the two loops approach each other closely, one must be cautious because the EA solution could be unphysical, as shown later in this article. It is somewhat unexpected that, when the system is closer to the azeotropic point, the sGC method takes fewer steps for convergence. One explanation might be that, near the azeotropic point, the vdW loop is much flatter and thus the sGC transform is much faster. 4.2. n-Butane + Water Mixture. The n-butane + water mixture was studied in detail by Hanif et al.4 with one new EA method. The details of their EA method can be found in the relevant publications.6−8,11 Our calculation was based on the Peng−Robinson EOS with the standard vdW mixing rule and the binary parameter k12 = 0.5. The vdW mixing rule is quite good for water solubility in a hydrocarbon phase, but gives hydrocarbon concentrations in the water phase that are a few orders of magnitude lower than real values.17 Nevertheless, we focus on the calculation algorithm of phase equilibrium, which 9693



Article

changed tremendously during the process. In the beginning, there were only two vdW loops, and the loops exhibited enormous oscillations with the progress of the sGC transform. In the meantime, the oscillation amplitudes were suppressed tremendously. At stage 12, the profile was visually horizontal, and its two intersection points with the original g′(x) function represent the two points of phase equilibrium. There is only one horizontal line at this pressure, whereas the original EA method gives two horizontal lines, or two phase equilibria. We note that Figure 7 shows much stronger oscillation than either Figure 4 or Figure 9 for the next mixture considered, which might indicate that the n-butane + water mixture is highly nonideal because the solubility of n-butane in water is as low as 10−12. 4.3. Water + 1-Pentanol. The water + 1-pentanol system was modeled by PC-SAFT with associating interactions,18 and its phase diagram was calculated within the temperature range from 298.15 to 423.15 K at the fixed pressure of P = 0.1013 MPa. As shown in the inset of Figure 8, this mixture displays a

Figure 6. Derivative of the Gibbs free energy for n-butane + water at T = 344.26 K and P = 0.90 MPa, as well as (inset a) 0.80 and (inset b) 0.85 MPa. (Insets c,d) Local views for two EA methods. Solid curves and horizontal lines are for the MFT and sGC solutions, respectively.

VLE, VLLE and LLE but missed VLE2, as one consequence of the assumption of a single phase equilibrium. One can deduce that there is a transition from inset b to inset c of Figure 6 or a VLLE line that brings the double phase equilibria to a single phase equilibrium. Unfortunately, such a transition is unknown in advance, and one has to resort to a plot of g(x) and g′(x), as instructed by the new algorithm. It is anticipated that the EA method, which is typical of other mathematical manipulations, requires the initialization of the compositions and that the initialization and subsequent calculation could become much more complicated when more vdW loops are present. In contrast to the two earlier classical methods, the implementation of the sGC method is straightforward, and it experiences no discrepancies in phase equilibrium calculations for different pressures or different phase behaviors, such as those at P = 0.80, 0.85, and 0.90 MPa in Figure 6. The sGC method simply follows the transform in eq 4 successively until gi+1(x) is stable. It does not require g(x) and g′(x) to be plotted or the composition to be initialized. As one example, we chose the sGC transform at P = 0.90 MPa in Figure 7 to illustrate a few steps in the process. One can see that the g′(x) profile

Figure 8. Phase diagram calculated by the SS (dotted), EA (dashed), and sGC (solid) methods for water + 1-pentanol at P = 0.1013 MPa, as modeled by PC-SAFT. The inset is a broader view of the diagram.

variety of phase behaviors, including LLE, double VL equilibria, and single VLE with increasing temperature, which is somewhat similar to Figure 5 with increasing pressure at a fixed temperature. Our calculations show that the SS method is numerically very unstable in the temperature range of 367.5− 372.5 K. Two distorted curves were found, and consequently, there was no way to determine the phase equilibrium. Such a failure, as noted in previous cases, occurs near the region for a phase transition from LLE to double VL equilibria. The EA method, after assuming double VL equilibria for double vdW loops, gives correct results above the VLLE line. However, this assumption is misleading for phase equilibrium below the VLLE line, where an unphysical VLLE line is created. These observations are quite similar to those obtained for the nbutane + water mixture at high pressure. Although a resolution can be obtained in the same way as by Hanif et al.,11 some detailed graphic analysis of the vdW loops is needed. The sGC method was found to be as robust as before, and no abnormal calculations were required for the region of multiple phase transitions. The VLLE line was found at T = 369.51 K, and its formation is shown in Figure 9 with the sGC transform. Essentially, the system proceeds from two vdW loops and shows a rapid increment in the oscillation frequency and

Figure 7. sGC transform for n-butane + water at T = 344.26 K and P = 0.90 MPa. The solid curve, three dashed curves, and solid horizontal line are for MFT and transforms at three stages and at stage 12, respectively. 9694



Article

bound to the phase equilibrium. When the temperature continues to drop to 600 K, as shown in Figure 10d, the bcc segment is replaced by an fcc segment, and remarkably, this gives rise to three vdW loops. However, our sGC solution shows again that only one phase equilibrium exists, with the two compositions of 0.001556 for the fcc solid and 0.9986 for the liquid. At the lowest temperature in our test of T = 300 K, the three segments in Figure 10e represent the fcc, hcp, and bcc solids, and the liquid no longer participates. The two equilibrium values are 5.21 × 10−6 for the fcc solid and 1 × 10−8 for the bcc solid, noting that the hcp segment remains. Compared with earlier examples, one can see that the fluctuations in panels a and b of Figure 10 are moderate, quite similar to those in Figures 2 and 4, and are very strong in panels c−e of Figure 10, which are more like in Figure 7. This suggests that this alloy is highly nonideal at low temperatures. In fact, the alloy model developed through some polynomialtype fittings from segment to segment, quite different from the patterns observed in our classical thermodynamics examples. Figure 6f shows the full phase envelope from solid to liquid. In these calculations, just as simply as in the cases discussed in sections 4.1−4.3, no special treatments are needed for the number and type of phase coexistences. The complexity of g(x) is instinctively taken into account by the sGC transform or, in other words, by fluctuations of the mixture composition. The two classical methods, as byproducts of the sGC transform, can face complexities such as the determination of the type and number of phases, the initialization of the composition, and numerical jumps among multiple segments of g(x). These problems are, in particular, serious around multiple phase transitions from one type to another, as illustrated well by the examples discussed in this work.

Figure 9. sGC transform for water + 1-pentanol at P = 0.1013 MPa and T = 369.51 K. The solid curve, two dashed curves, and solid horizontal line are for MFT and transforms at stages 2, 4, and 12, respectively.

decrement in the amplitude with the progress of the sGC transform. Compared with Figure 7, the present mixture has relatively weaker oscillations, which can be observed by comparing their profiles at stage i = 4. This indicates that the mixture is mildly nonideal, which is supported by the fact that the solubility of 1-pentanol in water is at the level of 10−3, in contrast to 10−12 for the solubility of n-butane in water. 4.4. Ca + Na Mixture. Alloys have very complex phase behaviors and, with changing temperature, could become liquids or solids with three types of structures, namely, facecentered cubic (fcc), body-centered cubic (bcc), and hexagonal close packed (hcp). Therefore, it can be deduced that there could have as many as 10 types of phase envelopes, in contrast to three types (VVE, VLE, LLE) in classical thermodynamics. In addition, these phase envelopes can coexist, which makes the calculation of phase equilibrium much more complicated. Therefore, the classical methods can face significant difficulties, and a very tedious analysis of the phase behavior has to be made in advance. As one illustrative example to show the universality of the sGC method, we chose the Ca + Na mixture, which was modeled by Zhang et al.19 for g(x). g(x) is highly temperature-dependent and thus creates a variety of phase behavior. These phase behaviors are shown in Figure 10a−e with decreasing temperatures. Panels a and b of Figure 10 at two high temperatures are quite similar to Figures 2 and 4 presented earlier, where g′(x) is relatively simple and yields a single LLE and double LL equilibria, noting that double LL equilibria are seldom found in classical thermodynamics. When the temperature is decreased to 800 K, the left side of liquid g′(x) is overtaken by two segments, namely, bcc and hcp, and the mixture exhibits two vdW loops, which can easily lead to the conclusion of the existence of double phase equilibria as found in Figure 10c. The sGC calculation, however, shows the existence of a single phase equilibrium, with equilibrium compositions of 0.0173 for the bcc solid and 0.9886 for the liquid. One interesting point is that the hcp segment looks irrelevant to the phase equilibrium between bcc and liquid. However, this appearance is superficial, because if the hcp curve were moved up or down, the fluctuation pattern of the compositions would be different and would give different phase equilibrium results. Therefore, this fact supports our view that the vdW loop, despite its artificial appearance, is physically

5. CONCLUSIONS This article has illustrated the application of the sGC method for complex phase equilibrium calculations of binary mixtures. To demonstrate the capability of the sGC method, we selected several highly nonideal mixtures, namely, pyridine + acetic acid, n-butane + water, water + 1-pentanol, and Ca + Na metal, for phase diagram calculations. These mixtures were modeled by the van Laar ACM, the PR EOS, PC-SAFT, and an alloy-based model, respectively. The phase behaviors of these mixtures include VLE, LLE, SLE, SSE, and their coexistences. Two classical methods, namely, the SS and EA methods, were also studied for comparison. It was demonstrated that the sGC method is robust and universal for phase equilibrium calculations, whereas the two classical methods suffer a number of problems around multiple phase transitions. These problems include numerical instability, unphysical phase transitions, and the absence of phase equilibrium, which are caused by initialization; the mathematical algorithm; and, more fundamentally, the coexistence of multiple vdW loops. The sGC method, which was developed from our theory of composition fluctuations, completely avoids these problems. As another support for the sGC method, we demonstrated (see the Appendix) that its phase stability criterion is identical to that of the TPD method. In practical applications, the computation time is an important factor of consideration. The computational cost of the sGC method is generally higher than those of the SS and EA methods, because of its consideration of a wide range of composition fluctuations. However, the level of cost depends highly on the case, as well as the discrete number of 9695



Article

Figure 10. (a−e) sGC transforms and (f) phase diagram for Na + Ca, as modeled by polynomial correlation.19 The solid curves and horizontal lines are for the MFT and sGC solutions, respectively. The dashed curves are for transforms at (a,b) stage 3 and (c−e) stage 5.

the sGC method can, in principle, be applied to the transform of the Gibbs free energy by making fluctuations of, for example, two composition variables for a three-component mixture. This requires further investigation because of the factor of time cost, which is proportional to N2, or the square of the number of discrete points. Nevertheless, it can act as one tool to check the newly developed flash algorithm, which, when reduced to binary mixtures, should match the sGC method in terms of phase equilibrium calculations.

compositions. For example, using the numerical parameters in section 3, the cost of the sGC method is ∼12 times higher than that of the SS method and 2 times higher than that of the EA method for the phase diagram calculation in Figure 5, which was modeled by the cubic PR EOS. For the noncubic PC-SAFT EOS in Figure 8, the costs of the three methods are comparable because the main computational burden is shifted to calculating the large number of property derivatives and solving the association term of the SAFT model. Note that, in these comparisons, the correctness of the physical solution was not taken into account. This fact, as well as the facts that neither initialization nor property derivatives are needed, suggests that the sGC method is more favorable for noncubic EOS. Overall, the sGC method is the best for finding the correct solution when complex phase behaviors are encountered. It is particularly useful to build up the phase diagram when time is not a significant concern. Regarding multiple components,

■

APPENDIX: INVESTIGATION OF PHASE STABILITY BY THE SGC METHOD

Equation 4 can also be used to study phase stability, which was not mentioned previously in our work. To ensure phase stability of a binary mixture at composition z, Michelsen1 confirmed that the tangent plane distance (TPD) function 9696


Industrial & Engineering Chemistry Research TPD(x) = g0(x) − xμ1(z) − (1 − x)μ2 (z) ≥ 0

Article

(11) Quinones-Cisnerosa, S. E.; Deiters, U. K. An efficient algorithm for the calculation of phase envelopes of fluid mixtures. Fluid Phase Equilib. 2012, 329, 22. (12) Tang, Y. A new grand canonical ensemble method to calculate first-order phase transitions. J. Chem. Phys. 2011, 134, 224508. (13) Tang, Y. A new method of semigrand canonical ensemble to calculate first-order phase transitions for binary mixtures. J. Chem. Phys. 2012, 136, 034505. (14) Kofke, D. A. Semigrand canonical Monte Carlo simulation; integration along coexistence lines. Adv. Chem. Phys. 1999, 105, 405. (15) Errington, J. R.; Shen, V. K. Direct evaluation of multicomponent phase equilibria using flat-histogram methods. J. Chem. Phys. 2005, 123, 164103. (16) Michelsen, M. L.; Mollerup, J. M. Thermodynamic Models: Fundamentals and Computational Aspects; Tie-Line Publications: Holte, Denmark, 1992. (17) Economou, I. G.; Tsonopoulos, C. Associating models and mixing rules in equation of state for water/hydrocarbon mixtures. Chem. Eng. Sci. 1997, 52, 511. (18) Gross, J.; Sadowski, G. Application of the perturbed-chain SAFT equation of state to associating systems. Ind. Eng. Chem. Res. 2002, 41, 5510. (19) Zhang, S.; Shin, D.; Lin, Z. K. Thermodynamic modeling of the Ca−Li−Na system. CALPHAD: Comput. Coupling Phase Diagrams Thermochem. 2003, 27, 235.

(A1)

must be satisfied for all values of composition x in the range [0, 1]. Application of eq 4 to eq A1 gives, at the first stage of the transform ⎡ g (x + y ) + g (x − y ) ⎤ 0 ⎥ g1(x) = min ⎢ 0 ⎥⎦ 0 < y < cx⎢ 2 ⎣ =

g0(x + ym ) + g0(x − ym ) 2

≥ xμ1(z) + (1 − x)μ2 (z) (A2)

and therefore g1(z) = g0(z) = zμ1(z) + (1 − z)μ2 (z)

(A3)

where ym is the fluctuation magnitude required to give the minimum of eq 4. It can be deduced that the relations A2 and A3 exist for the whole process of the sGC transform, and thus, g0(z) = g1(z) = ...... = g∞(z), or in other words, the phase is stable. Conversely, if gi+1(z) ≠ gi(z) at one stage i, then phase stability at z is broken. Therefore, phase stability can be simply determined during the process of sGC transformation, prior to phase equilibrium. Physically, this indicates that the combination of subsystems could break phase stability.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: 1-5196406515. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS The authors thank John Wasson for useful suggestions in the revision of this article. REFERENCES

(1) Michelsen, M. L. The isothermal flash problem, Part I. Stability. Fluid Phase Equilib. 1982, 9, 1. (2) Michelsen, M. L. The isothermal flash problem, Part II. Phasesplit calculation. Fluid Phase Equilib. 1982, 9, 21. (3) Hanif, N. S. M.; Shyu, G, S.; Hall, K. R.; Eubank, P. T. Areas of J. C. Maxwell for pure-component fluid-phase equilibria from equations of state. Ind. Eng. Chem. Res. 1996, 35, 2431. (4) Hanif, N. S. M.; Shyu, G. S.; Hall, K. R.; Eubank, P. T. Calculation of multiphase equilibria using the equal area rule with application to hydrocarbon/ water mixtures. Fluid Phase Equilib. 1996, 126, 53. (5) Eubank, P. T.; Hall, K. R. Equal area rule and algorithm for determining phase compositions. AIChE J. 1995, 41, 924. (6) Tang, Y.; Stephenson, G.; Zhao, E.; Agrawal, M.; Saha, S. Calculations of Complex Phase Equilibrium by Equal-Area Method. AIChE J. 2012, 58, 591. (7) Heidemann, R. A.; Michelsen, M. L. Instability of Successive Substitution. Ind. Eng. Chem. Res. 1995, 34, 958. (8) Zhu, Y.; Wen, H.; Xu, Z. Global stability analysis and phase equilibrium calculations at high pressures using the enhanced simulated annealing algorithm. Chem. Eng. Sci. 2000, 55, 3451. (9) Rossi, C. C. R. S.; Cardozo-Filho, L.; Guirardellob, R. Gibbs free energy minimization for the calculation of chemical and phase equilibrium using linear programming. Fluid Phase Equilib. 2009, 278, 117. (10) Giovanoglou, A.; Galindo, A.; Jackson, G.; Adjiman, C. S. Fluid phase stability and equilibrium calculations in binary mixtures Part I: Theoretical development for non-azeotropic mixtures. Fluid Phase Equilib. 2009, 275, 79. 9697


Calculations of Complex Phase Equilibrium by ... - ACS Publications

Recommend Documents