Ind. Eng. Chem. Res. 2001, 40, 2149-2159
2149
A Model-Free Approach Data Treatment of Vapor-Liquid Equilibrium Data in Ternary Systems. 2. Applications Elizabeth Lam,†,‡ Andre´ s Mejı´a,† Hugo Segura,*,† Jaime Wisniak,§ and Sonia Loras| Department of Chemical Engineering, Universidad de Concepcio´ n, P.O. Box 160-C, Correo 3, Concepcio´ n, Chile, Department of Chemical Engineering, Ben-Gurion University of the Negev, Beer-Sheva, Israel, and Departamento de Ingenierı´a Quı´mica, Facultad de Quı´mica, Universitat de Vale` ncia, 46100 Burjassot, Valencia, Spain
In the first part of this work, a model-free technique for treating ternary systems was presented. The procedure avoids the use of empirical activity coefficient models by generating numerically vapor-phase concentrations, activity coefficients, and excess Gibbs energy values directly from [T, P, x] data. Here we illustrate the use of the method as a complementary tool for modeling and assessing multicomponent vapor-liquid equilibrium data. Guidelines are given for fitting bubble surfaces and for treating azeotropic systems. In addition, reasonable approximations are outlined for handling ternary systems for which excess volume or enthalpy data are not available. Examples are presented for experimental isothermal and isobaric systems. Introduction In the first part of this work,1 we have presented two numerical procedures, based on a model-free approach, for treating vapor-liquid equilibria (VLE) of ternary systems. Both methods are based on the Barker2 equation, written as a nonlinear partial differential equation on the excess Gibbs energy (GE). The first method, called mixed approach (MA), solves the partial differential equation using a finite difference technique, where the GE derivatives are approximated by single and central differences. The second method, called second-order approach (SOA), is also based on a finite difference technique, but it uses central differences for approximating the GE derivatives. In both approaches the discretization of the partial differential equation proceeds over a triangular grid of N × (N + 1)/2 nodes that represents the concentrations of a ternary system. Thus, the partial differential equation yields a set of nonlinear algebraic equations that must be solved in every node of the discretization grid. Numerical stability analysis of the methods indicates that MA is useful for obtaining the GE function of nonazeotropic ternary systems, while SOA is effective for treating azeotropic systems, although it demands more computational effort. MA and SOA are both applicable when the bubble surface of a ternary system is smoothed by a function that depends on the liquid-phase composition, xi. In addition, for isothermal systems, it is necessary to know the excess volume of the liquid phase (V ˜ E) as a function of the composition, while for isobaric systems, it is required to know the excess enthalpy (H ˜ E). The fitting of bubble surfaces is a central problem in the treatment of VLE by means of model-free approaches, because this is the input information for calculating equilibrium properties. To obtain accurate * Corresponding author. E-mail:
[email protected]. † Universidad de Concepcio ´ n. ‡ Present address: Department of Chemical Engineering, Universidad Cato´lica del Norte, Avenida Angamos 0610, Antofagasta, Chile. § Ben-Gurion University of the Negev. | Universitat de Vale ` ncia.
results, the fit must represent the data within the range of experimental errors. In addition, as shown in our previous paper,1 the curvature of the bubble surface in the vicinity of an azeotrope is constrained by geometrical relations. When the curvature constraints are violated, no GE value is able to represent the azeotrope, which should then be qualified as nonphysical. The satisfaction of curvature constraints is a convergence requirement for every model-free approach. Consequently, the fit must also match these constraints when the VLE data under consideration are azeotropic. In this second part, guidelines are given for applying the model-free techniques to ternary systems. Special attention is given to the fit of bubble surfaces and to the selection of the approach (MA or SOA). In addition, reasonable approximations are outlined for handling ternary systems for which excess volumes or enthalpies are not available. Fit of Bubble Points in Ternary Systems The basic relation for solving the GE function is the nonlinear Barker differential equation, given by1
{
[( ) ] { [( ) [( ) ]} { ] [( ) ]} [( ) ] [( ) ]}
x1P01 ∂gE P) exp gE + (1 - x1) Φ1 ∂x1 x2
∂gE ∂x2
x1
∂gE ∂x2
E
x1
+ ψ x1 -
x2P02 ∂gE + exp gE - x1 Φ2 ∂x1
+ ψx2
ψx1 + (1 - x2)
x2
∂g ∂x1
x2
x1
+ ψ x2
+
+ ψx1 - x2
+
x2
x3P03 exp gE Φ3
∂gE ∂x2
x1
+ ψx2
(1)
where P is the vapor pressure, gE is the dimensionless excess Gibbs energy (G ˜ E/RT), Φi is a fugacity correction of the phases defined as
[
φˆ i V ˜ Li Φi ) exp (P - P0i ) φi RT
10.1021/ie000753u CCC: $20.00 © 2001 American Chemical Society Published on Web 04/06/2001
]
(2)
2150
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
and ψxi is the following function of the excess volume and enthalpy of the liquid phase: E
( )
H ˜ ∂T ψxi ) 2 ∂x RT i
E
( )
V ˜ ∂P RT ∂xi xj*i
(3)
xj*i
When T, P, ψx1, and ψx2 are known functions of the liquid-phase mole fractions, then both MA and SOA are able to generate gE information, as explained in depth in the first part of this work.1 For isothermal systems, the basic information required is a fit of the bubble-point pressures with the composition xi, as well as data on volumes of mixing. For isobaric systems, the calculation requires a fit of the bubble-point temperatures with composition and mixing enthalpy data. In this section we propose fitting functions for pressure or temperature in ternary systems. Isothermal Systems. Activity coefficients γ may be directly calculated from experimental data using the definition3
γi )
yiPΦi
(4)
xiP0i
Equation 4 yields the following relation for ternary systems:
of the Redlich-Kister5 equation: 3
xixj(Aij + Bij{xi - xj} + Cij{xi - xj}2 + ...) + ∑ j>i x1x2x3(Q0 + x1Q11 + x2Q12 + ∑ [Qk1{x1 - x2}k + kg2
Ω)
Qk2{x1 - x3}k + Qk3{x2 - x3}k]) (10) where Aij, Bij, Cij, etc., are binary constants, which can be obtained by fitting the data of the constituent binaries. Parameters Q0, Q11, Q12, etc., correspond to ternary constants that, in turn, can be calculated from ternary data. For isothermal systems, eq 3 (used in the discretization of the Barker relations) becomes
gE )
yi
3
3
xi ln γi ) ln P + ∑xi ln + ∑xi ln ∑ x i)1 i)1 i)1 i
Φi P0i
∆ ln P ) (5)
∑ i)1
xi ln
P0i
[ ]
∆H ˜ 0i T0i xi 1) 0 T i)1 RT i 3
)
P
∑
3
xi ∑ i)1
so that 3
ln P )
∆H ˜ 0i
-
1
3
∑xi
Ti)1
RT0i
∆H ˜ 0i R
(12)
3
xi ln P0i + gE - ω - ∑xi ln Φi ∑ i)1 i)1
(6)
where ω corresponds to the weighted volatility function4
ω)
xi ln ∑ x i)1
(7)
i
Equation 6 yields the following expression for the mixing logarithm of the bubble pressure (∆ ln P): 3
xi ln ∑ i)1
˜ 0i are the boiling temperatures and where T0i and ∆H heats of vaporization of the pure components at the pressure of the system. From eq 12 we have 3
yi
3
∆ ln P ) ln P -
(11)
xj*i
If no experimental data are available for the mixing volumes of the ternary system, we can assume that ψxi ≈ 0 because the ratio V ˜ E/RT is almost zero ( 0 4∆V ˜ Px1x1 + G
(20)
while the temperature curvature of a isobaric binary azeotrope is constrained by
∆H ˜ T +G ˜ yV1y1 > 0 Az x1x1 T
-4
(21)
When an azeotropic-like state violates the constraining relations given by eqs 18-21, then no model-free approach is able to calculate the corresponding value of gE. It should be pointed out that this type of nonphysical azeotrope may be predicted by eqs 9 and 14, when applied to azeotropic systems characterized by flat bubble surfaces. For such cases, the model-free approaches discussed here will either generate a noisy gE function in the vicinity of the problematic azeotrope11 or, simply, will not converge for the ternary system. Consequently, it is always recommended to test the fit of bubble surfaces before applying a model-free technique. When the fit violates the physical constraints, eqs 18-21 are introduced as penalty functions, and the parameters of eq 10 should be redetermined using constrained optimization techniques. Physical constraints affect exclusively positive deviation and saddle azeotropes,1 which may be found by the simultaneous solution of the relations
∂T ∂P ) )0 ∂x1 ∂x1 ∂P ∂T ) )0 ∂x2 ∂x2
(22)
applied to eq 9 or eq 14. Once the azeotropic concentra-
2152
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
tions have been calculated, it is possible to determine the deviation from ideal behavior by inspecting the signs of the eigenvalues λ1 and λ2 of the matrices
[
Px1x1 Px1x2 Px1x2 Px2x2
[
]
Tx1x1 Tx1x2 Tx1x2 Tx2x2
Discussion, Applications, and Examples
isothermal case
]
isobaric case
(23)
Positive deviation azeotropes are characterized by two negative eigenvalues, while saddle azeotropes are characterized by eigenvalues with mixed signs. The different properties that appear in eqs 18-21 can be calculated as follows:
∆V ˜ )
(ZV - ZL)RT RT ≈ P Pf0 P
∆H ˜ ) -H ˜E+
∑xi∆H˜ Vi ≈ ∑xi∆H˜ Vi
(24)
In eq 24, Z is the compressibility factor, ∆H ˜ V is the heat of vaporization of the pure component, and H ˜ E is the excess enthalpy of the liquid phase. All of the properties in eq 24 must be evaluated at azeotropic conditions, as determined from eq 22. In addition, the Gibbs energy function for the vapor phase is given by
G ˜V )
∑yi{RT ln[yiφˆ i] + G˜ X,V i } ≈ Pf0 RT∑yi{RT ln yi + G ˜ X,V i }
(25)
corresponds to a reference state of Gibbs where G ˜ X,V i energy for the pure component, which depends on pressure and temperature, and the φˆ i correction can be calculated from an equation of state. The composition derivatives of eq 25 yield the elements of the Hessian matrix of the vapor-phase Gibbs energy which, for ideal gas behavior, is given by the following relation for a ternary azeotrope:
[
y1 + y3 1 y1y3 y3 H(G ˜ V) ) RT y2 + y3 1 y3 y2y3
]
(26) yi)xiAz
Equations 20 and 21 should be also considered for constraining the fit of the bubble curve in binary systems because, otherwise, no solution may be obtained for the ternary system. We have found some cases in which the bubble curve of binaries is not well correlated by such a constrained fit. This situation was detected especially in highly nonideal binary azeotropic systems. For such problematic binaries, we recommend the following strategy: the Wilson model12 may be used for performing a previous smoothing of the bubble curve. The Wilson model parameters can be calculated using the objective function Nd
OF )
problematic binaries is used for fitting the pertinent binary parameters of eq 10.
2 2 exptl [|Pexptl - Pcalc - ycalc ∑ j j | + |yj j | ] j)1
(27)
Then, instead of using directly the binary experimental data, the Wilson prediction of the bubble curve of these
Systems That Exhibit Complex Topology in the Bubble Surface. The capabilities of the procedures proposed here will be illustrated first by the following fictitious isothermal ternary system, where its vapor pressure is assumed to be represented by the following empirical relation:
P/kPa ) 103x1 + 102x2 + 100x3 + 3
xixj(10 - 40{xi - xj} + 10{xi - xj}2) ∑ j>i
(28)
Equation 28 corresponds to a case where all of the binaries are polyazeotropic, each one exhibiting two azeotropes. Furthermore, the system exhibits three ternary azeotropes, as shown in Figure 1a. Each one of these ternary azeotropes satisfies the constraining relation for the curvature of the bubble pressure surface and thus has physical significance. In this example, we will assume that the vapor phase is ideal (Φi ) 1) and that ψxi ≈ 0. Equation 28 can be used to generate VLE P-x data for all of the binaries. These binary data are useful for analyzing the capability of common gE models (previously fitted to binaries) for predicting the topology of the ternary system. Wisniak et al.13 have indicated that the NRTL model12 is adequate for fitting the data of binary polyazeotropic systems. Table 1 presents the pertinent parameters, correlation statistics for the binaries, and prediction statistics for the ternary system. Although from Table 1 we could conclude that the NRTL model gives a good representation of the data, Figure 1b reveals that it fails in yielding the correct topology of the ternary azeotropes. We can conclude then that adequate prediction statistics are not enough to guarantee that the model predicts the correct topological details of the data. The example given here corresponds to a family of systems whose thermodynamic behavior is difficult to predict from binaries, and thus it constitutes a good test for data analysis using a model-free approach. As mentioned, azeotropes are present in all of the constituent binaries of the ternary system. Consequently, the MA is not recommendable for treating these data. Let us consider the SOA. Figure 1c shows how the composition of the vapor phase converges to a fixed value (within an error 10-4) for N ) 100. Therefore, the properties of the system will be generated using N ) 101 (this is usually a reasonable value for many systems treated with the SOA and can be handled efficiently by a personal computer). Figure 1d shows the Gibbs energy surface of the system and that it combines ranges of positive and negative deviations from ideal behavior. Such a topology agrees well with the polyazeotropic scenario of the binary systems,13 particularly when the vapor phase can be considered ideal. Figure 1d shows the calculated phase diagram for the ternary system and the correct prediction of the topology of each ternary azeotrope (Figure 1a). The example given here demonstrates that the suggested model-free approach is able to generate the excess energy, the activity coefficients, and the vaporphase mole fractions of the ternary system, a difficult
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2153
Figure 1. (a) Isobars for the fictitious system predicted by eq 28: (A) maximum pressure azeotrope; (B) saddle azeotrope; (C) minimum pressure azeotrope. (b) Isobars for the system in Figure 1a, as predicted from binaries using the NRTL model: (A) saddle azeotrope; (B) saddle azeotrope; (C) minimum pressure azeotrope. (c) Vapor-phase mole fraction convergence of the SOA for the system in eq 28. (d) Excess Gibbs energy surface calculated for the system in eq 28. Data generated using the SOA with N ) 101: (‚‚‚) ternary mixture; (s) binaries. (e) Phase diagram calculated from the SOA for the system in eq 28. (f) Vapor-phase residuals for the system in eq 28, according calc to the NRTL model prediction. yi* corresponds to the theoretical value obtained from the model-free approach. yi is the NRTL prediction for vapor-phase mole fractions.
2154
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
Table 1. NRTL Parameters, Correlation and Prediction Statistics for the VLE Data Generated from Eq 28a system
Aij/(J mol-1)
Aji/(J mol-1)
Ri,j
∆P/%b
1+2 1+3 2+3 1 + 2 + 3c
-3268.64 6671.30 6645.81
-3252.26 -3259.39 6670.12
0.3 0.3 0.3
0.35 0.38 0.37 1.14
a
Twenty equally spaced concentration points have been considered for each binary, and 175 equally spaced concentration d points have been considered for the ternary. b ∆P ) 100/Nd∑N i | calc exptl Pexptl P |/P , where N is the number of data points. d i i i c Prediction from binaries.
task when using conventional gE models with parameters fitted to binaries. The required input data are a correlation function for the bubble surface P(x1,x2) (such as eq 28), which can be determined from VLE measurements in a static equilibrium cell. Because the starting point of the model-free approach is the Gibbs-Duhem equation, the method will automatically generate an accurate set of properties that is fully consistent with the P(x1,x2) fit equation. Let us assume that eq 28 has been generated from accurate isothermal P-x data measured in a static cell and that the vapor-phase concentration has not been measured. The method proposed here will generate the corresponding theoretical concentrations of the vapor phase. These may be used for testing a particular gE model by calculating the residuals based on the predicted and theoretical values of yi. This application is illustrated in Figure 1f, where the predictive capability of the NRTL model is tested. It is seen that residuals are not randomly distributed and that they are autocorrelated and biased to overpredict the composition of the vapor phase. We conclude then that the NRTL model gives an inaccurate prediction of the data of this system. Treatment of Simple Isothermal Ternary Systems. The treatment of nonazeotropic systems may be performed using the MA discussed before. A representative example is the ternary mixture 1-heptene (1) + heptane (2) + octane (3), reported by Kudryavtseva et al.14,15 at 328.15 K. The data of the ternary system and its constituent binaries have been smoothed according to eq 17, and the corresponding parameters and correlation statistics are presented in Table 2. The smoothing relation shows that the system exhibits positive deviation from ideal behavior and that no azeotrope is present. Vapor-phase fugacity corrections have been neglected (Φi ) 1), and ψxi is assumed to be zero. Figure 2a illustrates the convergence of the composition of the vapor phase for both model-free approaches. From the figure we conclude that the MA requires a discretization grid with N ∼ 500. A broader grid can be used when considering the SOA, yielding the same precision. However, the MA runs faster and requires less physical memory than the SOA. Another important conclusion from Figure 2a is that both approaches converge to the same equilibrium values. Figure 2b shows the Gibbs energy surface for the ternary system and that each constituent binary deviates positively from ideal behavior. In addition, the GE surface exhibits a saddle geometry, a fact that cannot be predicted easily from the deviation observed in binaries. When saddle points appear in the GE surface, unexpected local negative deviations are possible for the ternary system.
Table 2. Parameters and Correlation Statistics for the Vapor Pressures of the System 1-Heptene (1) + Heptane (2) + Octane (3) at 328.15 K (Data of Kudryavtseva et al.14,15) Binary Parameters for Eq 10 system/i + j
Aij
1+2 1+3 2+3
Bij
∆P/%a σP/kPab
Cij
0.099 355 0.000 000 0.000 000 0.896 178 -0.356 716 0.217 630 0.631 803 -0.199 507 0.000 000
0.05 0.26 0.14
0.01 0.06 0.03
Ternary Parameters for Eq 10 system
Q0
Q11
Q12
∆P/%a
1+2+3
-2.354 20
2.572 07
1.858 40
0.72
Pure-Component Vapor Pressures Used in Eq
a
σP/kPab 0.19 9c
P01/kPa
P02/kPa
P03/kPa
27.49
23.09
8.43
exptl exptl d ∆P ) 100/Nd∑N - Pcalc , i |Pi i |/Pi of data points. b Standard deviation.
Percentage deviation, where Nd is the number c Taken from the data of Kudryavtseva et al.14,15
As indicated above, calculation of the equilibrium properties of an isothermal system using the model-free approach requires fitting the bubble pressure surface with the composition of the liquid phase. The method is then capable of predicting the theoretical compositions of the vapor phase, which are consistent with the fitted P(x1,x2) surface. The data of Kudryavtseva et al. include experimental vapor-phase mole fractions, and these can be compared to the predicted ones as a measure of consistency. This comparison is shown in Figure 2c as vapor-phase mole fraction residuals. It is seen that, although the vapor-phase residuals satisfy the consistency criteria [100|δy| e 1] required by the test of Van Ness et al.,16 as modified by Fredenslund et al.,17 they do not scatter randomly about the zero value line. Treatment of Complex Isothermal VLE Data for Ternary Systems. In this example we consider the data of the ternary system acetone (1) + chloroform (2) + methanol (3), reported by Tamir et al.18 at 298.15 K. The system exhibits azeotropic behavior for each constituent binary and combines positive and negative deviation from ideal behavior. In addition, association plays an important role in the thermodynamic behavior of the system, because of the presence of an alcohol and associating functional groups in the other components. According to the chemical theory,19 positive deviations from ideal behavior may be explained in terms of selfassociation, while negative deviations may be attributed to cross-association between the components of a mixture. In the ternary in question, we can expect a complex association situation, because competitive association regimes (self- and cross-association) can take place depending on the concentration of each component. The binary acetone (1) + chloroform (2) has a minimum pressure azeotrope, and the deviation from ideality may be explained in terms of cross-association of its components. The remaining binaries, acetone (1) + methanol (2) and chloroform (2) + methanol (3), exhibit maximum pressure azeotropes, which, in turn, may be explained in terms of the extensive self-association of the alcohol. Practice has shown that associating systems are difficult to correlate and assess by means of common gE models, because they are based on a theoretical physical basis instead of on chemical interactions. The data of the ternary system and their constituent binaries have been smoothed using eq 9. It should be
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2155 Table 3. Parameters and Correlation Statistics for the Vapor Pressures of the System Acetone (1) + Chloroform (2) + Methanol (3) at 298.15 K (Data of Tamir et al.18) Binary Parameters for Eq 10 binary system parameters
1+2
1+3
2 + 3a
Aij Bij Cij Dij Eij Fij Gij Hij Iij Jij Kij ∆P/%b σP/kPac
-0.940 484 0.653 261 0.294 477 -0.609 311
0.889 966 -0.381 584 0.417 723 -0.319 304
1.487 360 -0.283 513 0.720 015 -0.074 742 -0.337 711 1.662 946 3.914 123 -3.298 105 -6.566 601 3.236 125 4.911 296 0.63 0.21
Q00 Q11 Q12 Q21 Q22 Q23 ∆P/%b σP/kPac
0.42 0.13
0.18 0.06
Ternary Parameters for Eq 10 -7.039 304 Q31 -4.408 728 5.074 548 Q32 9.240 015 8.564 383 Q33 0.000 000 -2.845 334 Q41 20.600 073 1.808 301 Q42 -24.769 853 0.000 000 Q43 0.000 000 1.25 0.42
Pure-Component Vapor Pressures Used in Eq 9d P01/kPa
P02/kPa
P03/kPa
30.73
26.46
17.02
a
Obtained from the prediction of the Wilson model fitted to exptl d experimental data. b Percentage deviation, ∆P ) 100/Nd∑N i |Pi exptl c Standard - Pcalc |/P , where N is the number of data points. d i i deviation. d Taken from the data of Tamir et al.18
Figure 2. (a) Convergence of vapor-phase mole fractions for the system 1-heptene (1) + heptane (2) + octane (3) at 328.15 K. (b) Excess Gibbs energy surface for the system 1-heptene (1) + heptane (2) + octane (3) at 328.15 K. Data generated using the mixed approach with N ) 601: (‚‚‚) ternary mixture; (s) binaries. (c) Vapor-phase mole fraction residuals for the system 1-heptene (1) + heptane (2) + octane (3) at 328.15 K. δyi ) yexptl - yi*, i where yi* corresponds to the theoretical value obtained from the model-free approach.
noted that the fit of bubble pressure surfaces with eq 9 is tedious work when treating noisy VLE data. No satisfactory constrained fit of the data (eq 20) was found for the binary system chloroform (2) + methanol (3). Consequently, it was first treated with the Wilson model, and the predicted vapor pressure curve was used for determining the binary parameters A23, B23, etc., of eq 17. All of the remaining parameters were determined directly from experimental data. Table 3 gives the pertinent parameters and correlation statistics, from which we can see that, in comparison with previous examples, a larger number of parameters are needed for fitting the bubble curve. In addition, the statistics in Table 3 suggest that the ternary data are not completely compatible with the binary data; however, the larger deviations from the fit are probably due to the presence of pressure outliers in the experimental data of the ternary system. The isobars of the fit appear in Figure 3a, where it is seen that the system exhibits azeotropes for the binaries (A-C) and a saddle, almost equimolar, ternary azeotrope (D). Because all of the binaries exhibit azeotropic behavior, we select the SOA for calculating the equilibrium properties. For integration purposes, we have assumed that ψxi ≈ 0 and calculated the fugacity corrections according to eq 2. The fugacity coefficients were estimated from the method of Hayden and O’Connell20 using the association and solvation parameters suggested by Prausnitz et al.21 Pure-component liquid volumes were approximated using Rackett’s correlation,22 and critical properties were taken from DIPPR.23
2156
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
Figure 3. (a) Isobars for the system acetone (1) + chloroform (2) + methanol (3) at 298.15 K predicted by eq 9: (A) minimum pressure azeotrope; (B and C) maximum pressure azeotropes; (D) ternary saddle azeotrope. (b) Excess Gibbs energy surface for the system acetone (1) + chloroform (2) + methanol (3) at 298.15 K. Data generated using the SOA with N ) 101: (‚‚‚) ternary mixture; (s) binaries. (c) Phase diagram calculated from the SOA for the system acetone (1) + chloroform (2) + methanol (3) at 298.15 K. (d) Residual plot for the ternary system acetone (1) + chloroform (2) + methanol (3) at 298.15 K: (O) vapor-phase mole fraction residuals; (b) vapor-pressure residuals.
Figure 3b presents the excess Gibbs energy surface for the system and the complex topology present. The system combines positive and negative deviations from ideal behavior that reflect in inflections on the GE surface. Figure 3c displays the calculated phase diagram and the binary and ternary azeotropes. The residuals of the vapor-phase mole fraction and vapor pressures for the constituent binaries and ternary systems are presented in parts d and f of Figure 3, respectively. From Figure 3d it is possible to observe that the vaporphase mole fraction residuals appear outside the range of consistency for all of the binaries, although vaporpressure residuals are reasonably distributed for systems 1 + 2 and 1 + 3. For the case of the system 2 + 3, vapor-pressure residuals are distorted as a consequence of the previous smoothing with the Wilson model. Results for binaries are in good agreement with the
point-to-point inconsistency reported by DECHEMA Chemistry Data Series24 for the quoted systems. From Figure 3f, it is possible to observe that the vapor-phase mole fraction residuals appear outside the acceptable range, while vapor-pressure residuals show a reasonable distribution about the zero line. In addition, although the residuals of the ternary system scatter about the zero line for acetone, they do not do so for chloroform as its concentration increases. By their nature, the data of Tamir et al. are very difficult to model, and no satisfactory fit was found in this work with common gE models. However, the present analysis, where no model has been considered, suggests that this ternary system and their constituent binaries are not very accurate. Treatment of Complex Isobaric VLE Data for Ternary Systems. In this illustration we consider the isobaric data reported by Loras et al.25 for the system
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2157
Figure 4. (a) Isothermals for the system 2-methyl-2-propanol (1) + methyl 1,1-dimethylpropyl ether (2) + 2,2,4-trimethylpentane (3) at 101.3 kPa predicted by eq 14: (A and B) minimum temperature azeotropes; (C) ternary saddle azeotrope. (b) Excess Gibbs energy surface for the system 2-methyl-2-propanol (1) + methyl 1,1-dimethylpropyl ether (2) + 2,2,4-trimethylpentane (3) at 101.3 kPa. Data generated using the SOA with N ) 101: (‚‚‚) ternary mixture; (s) binaries. (c) Phase diagram calculated from the SOA for the system 2-methyl-2propanol (1) + methyl 1,1-dimethylpropyl ether (2) + 2,2,4-trimethylpentane (3) at 101.3 kPa. (d) Vapor-phase mole fraction residuals for the system 2-methyl-2-propanol (1) + methyl 1,1-dimethylpropyl ether (2) + 2,2,4-trimethylpentane (3) at 101.3 kPa.
2-methyl-2-propanol (1) + methyl 1,1-dimethylpropyl ether (2) + 2,2,4-trimethylpentane (3) at 101.3 kPa. This ternary system exhibits minimum temperature azeotropes for the binaries 2-methyl-2-propanol (1) + methyl 1,1-dimethylpropyl ether (2), and 2-methyl-2-propanol (1) + 2,2,4-trimethylpentane (3). In addition, it presents a ternary saddle azeotrope, which is clearly seen in the contour plot of isothermals shown in Figure 4a. The bubble temperatures of the ternary system were fitted using eq 14. The pertinent parameters and correlation statistics appear in Table 4. The system 1 + 3 is azeotropic; thus, integration has been performed using the SOA. The excess enthalpy needed in eq 3 was estimated from the modified UNIFAC group contribution method,6 and the temperature derivative was calculated from the temperature fit given by eq 14, with the parameters shown in Table 4. Fugacity corrections
were calculated from refs 20-22 using the critical constants reported in ref 23. It should be realized that there are no major differences between the treatment of isothermal and isobaric data by means of the model-free approaches proposed here. However, in the case of isobaric systems, the vapor pressures of pure components depend on temperature and vary from node to node in the discretization grid. The dependence of pure-component vapor pressures on temperature can be handled using the Antoine parameters given in Table 4. Figure 4b depicts the calculated GE surface using N ) 101 and the fact that the system presents positive deviations from ideal behavior. In addition, Figure 4c shows the predicted phase diagram for the ternary system. Once again, the model-free approach allows us to analyze the residuals of the vapor phase as a measure
2158
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001
Table 4. Parameters and Correlation Statistics for the Bubble Temperatures of the System 2-Methyl-2-propanol (1) + Methyl 1,1-Dimethylpropyl Ether (2) and 2-Methyl-2-propanol (1) + 2,2,4-Trimethylpentane (3) at 101.3 kPa (Data of Loras et al. 24) Binary Parameters for Eq 18 binary system parameters
1+2
1+3
2+3
Aij Bij Cij Dij Eij Fij Gij Hij ∆T/Ka σT/Kb
-0.604 679 0.160 413 -0.280 516
-1.405 914 0.776 311 -1.300 477 0.324 381 0.423 339 0.007 307 -1.743 080 1.217 445 0.10 0.17
-0.165 053 0.071 784 -0.040 651
0.07 0.10
0.06 0.08
Ternary Parameters for Eq 18 -2.030 889 Q23 1.589 448 3.587 903 Q31 0.977 611 0.624 367 Q32 5.405 007 -1.225 610 Q33 -0.761 250 -2.980 535 0.06 0.08
Q00 Q11 Q12 Q21 Q22 ∆T/Ka σT/Kb
Pure-Component Vapor-Pressure Data and Additional Constants Needed in Eq 23c component
Ai
Bi
Ci
∆H0i d/(J mol-1)
1 2 3
14.8533 14.3501 13.8022
2649.89 3111.28 2999.10
96.69 39.52 45.96
41560 32658 32447
a Average absolute deviation, ∆T ) 100/N ∑Nd|Texptl - Tcalc|, d i i i where Nd is the number of data points. b Standard deviation. c The Antoine parameters are taken from the data of Loras et al.;24 ln P0i /kPa ) Ai - Bi/T/K - Ci. d Heat of vaporization calculated from the Antoine equation at the temperature of the system.
of the reliability of the data (Figure 4d). The figure shows an excellent distribution of the residuals. In fact, the data reported by Loras et al. satisfy the constraining relations for isobaric azeotropes given by eqs 18-21. In addition, the vapor-phase residuals scatter about the zero line within the acceptable criteria of consistency of a point-to-point test. Therefore, we can conclude that the reported data for the ternary are accurate and compatible with their constituent binaries. Conclusions In this work we have applied model-free techniques for calculating equilibrium properties of ternary systems. It is shown that the proposed approach is a valuable tool for analyzing modeling problems and for assessing the accuracy of the data. From a modeling viewpoint, it is possible to perform a detailed study about the compatibility of a ternary system with its constituent binaries. In contrast, when assuming a GE model to be valid, ternary data are usually predicted from binaries. In this last approach it is difficult to explain the deviations of the predictions, which may originate on inaccurate data, on model weakness, or with a combination of both aspects. Lenka and Anderko26 have addressed the problem of predicting lowpressure ternary data from binary contributions using an equation of state (EOS) approach to phase equilibrium. Ternary equilibrium has been predicted using a van der Waals type EOS, with composition-dependent mixing rules fitted to binaries. Composition-dependent
mixing rules have similar fitting capabilities when compared to traditional GE models for correlating binary VLE data. However, as pointed out by Lenka and Anderko, there are some cases in which it is not possible to represent simultaneously binary and ternary data using the same set of parameters. Another example is provided by the parameters of GE models reported by DECHEMA Chemistry Data Series,24 where specific nonbinary parameters are used for fitting ternaries. In most cases, when multicomponent systems are modeled, only binary interactions have to be considered because the effects of additional interactions are usually negligible in the low-pressure range. Such a simplification has a good background on results related to dispersion forces, although experience shows that polar or association forces require better theories. An accurate description of the behavior of complex ternary mixtures is important because it is necessary to account for more than one unlike interaction. Model-free techniques offer a potential alternative for treating and analyzing these data. Consistency procedures are well established for treating binary data; however, weaker approaches are available for assessing ternary systems. Known consistency approaches for ternary data are the method of McDermott and Ellis,27 which is a kind of integral test, the method of Wisniak,28 which is not based on the GibbsDuhem relation, and the method of Fonseca and Lobo,29 which is based on the Barker equation and on the Redlich-Kister GE model.5 The weaknesses of integral and GE model based consistency tests have been reviewed in depth elsewhere.30 The model-free approaches discussed in this work offer an additional tool that is not limited by modeling pitfalls and, in addition, are rigorously based on the Gibbs-Duhem equation for predicting equilibrium properties. Particularly, the analysis of vapor-phase concentration residuals is promising for testing the accuracy of the data. In addition, the analysis of the constraining relations offers a good opportunity for testing azeotropic data. Acknowledgment This work was partially financed by FONDECYT, Santiago, Chile (Projects 1990402 and 7990065), and by MEC, Valencia, Spain (Project PB96-0788). E.L. was sponsored by CONICYT, Chile. A.M. acknowledges a grant from SF, Concepio´n, Chile. Nomenclature Aij, Bij, Cij, ... ) binary parameters for eq 10 C ) number of components g ) dimensionless Gibbs energy (G ˜ /RT) G ) Gibbs energy H ) enthalpy N ) grid size P ) absolute pressure Q0, Q11, Q12, ... ) ternary parameters for eq 10 R ) universal gas constant T ) absolute temperature V ) volume x, y ) mole fractions of the liquid and vapor phases Z ) compressibility factor Greek Letters ∆ ) mixing or vaporization property γ ) activity coefficient δ ) residual
Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2159 Φ ) fugacity correction, defined in eq 2 ψ ) function defined in eq 3 φˆ ) fugacity coefficient Superscripts ∼ ) molar property Az ) azeotropic composition 0 ) pure-component saturation property E ) excess property exptl ) experimental value calc ) calculated value L ) liquid phase V ) vapor phase X ) reference state in eq 25 * ) theoretical value
Literature Cited (1) Lam, E.; Mejı´a, A.; Segura, H.; Wisniak, J.; Loras, S. A Model-Free Approach Data Treatment of Vapor-Liquid Equilibrium Data in Ternary Systems. 1. Theory and Numerical Procedures. Ind. Eng. Chem. Res. 2000, 40, 2134. (2) Barker, J. A. Determination of Activity Coefficients from Total Pressure Measurements. Aust. J. Chem. 1953, 6, 207-210. (3) Van Ness, H. C.; Abbott, M. M. Classical Thermodynamics of Nonelectrolyte Solutions With Applications to Phase Equilibria; McGraw-Hill: New York, 1982. (4) Malesinski, W. Azeotropy and Other Theoretical Problems of Vapour-Liquid Equilibrium; Interscience Publishers: London, 1965. (5) Redlich, O.; Kister, A. T. Algebraic Representation of Thermodynamic Properties and the Classification of Solutions. Ind. Eng. Chem. 1948, 40, 345-348. (6) Larsen, B.; Rasmussen, P.; Fredenslund, Aa. A Modified UNIFAC Group-Contribution Method for Prediction of Phase Equilibria and Heats of Mixing. Ind. Eng. Chem. 1987, 26, 22742286. (7) Gmehling, J.; Li, J.; Schiller, M. A. Modified UNIFAC Model. 2. Present Parameter Matrix and Results for Different Thermodynamic Properties. Ind. Eng. Chem. Res. 1993, 32, 178193. (8) Wisniak, J.; Tamir, A. Correlation of the Boiling Point of Mixtures. Chem. Eng. Sci. 1976, 31, 631-635. (9) Tamir, A. New Correlations for Fitting Multicomponent Vapor-Liquid Equilibria Data and Prediction of Azeotropic Behavior. Chem. Eng. Sci. 1981, 36, 1453-1465. (10) Sayegh, S. G.; Vera, J. H. Model-Free Methods for Vapor Liquid Equilibria Calculations. Binary Systems. Chem. Eng. Sci. 1980, 35, 2247-2256. (11) Wisniak, J.; Apelblat, A.; Segura, H. Application of ModelFree Computation Techniques and Limiting Conditions for Azeotropy. An Additional Assessment of Experimental Data. Chem. Eng. Sci. 1997, 52, 4393-4402. (12) Walas, S. M. Phase Equilibria in Chemical Engineering; Butterworth Publishers: Boston, 1985. (13) Wisniak, J.; Segura, H.; Reich, R. Polyazeotropy in Binary Systems. Ind. Chem. Res. 1996, 35, 3742-3758. (14) Kudryavtseva, L. S.; Kirss, H.; Eisen, O. Vapor-Liquid Equilibrium in 1-Heptene-Heptane, 1-Heptene-Octane, Hep-
tane-Octane, Benzene-Thiophene, Benzene- Heptane, and Thiophene-Heptane at 55°. Eesti NSV Tead. Akad. Toim., Keem., Geol. 1971, 20, 292-296. (15) Kudryavtseva, L. S.; Kirss, H.; Eisen, O. Calculation and Study of Vapor-Liquid Equilibriums in the Heptene-HeptaneOctane and Benzene-Thiophene-Heptane Ternary Systems. Eesti NSV Tead. Akad. Toim., Keem., Geol. 1972, 21, 19-23. (16) Van Ness, H. C.; Byer, S. M.; Gibbs, R. E. Vapor-Liquid Equilibrium: Part I. An Appraisal of Data Reduction Methods. AIChE J. 1973, 19, 238-244. (17) Fredenslund, Aa.; Gmehling, J.; Rasmussen, P. VaporLiquid Equilibria Using UNIFAC; Elsevier: Amsterdam, The Netherlands, 1977. (18) Tamir, A.; Apelblat, A.; Wagner, M. An Evaluation of Thermodynamic Analyses of the Vapor-Liquid Equilibria in the Ternary System Acetone-Chloroform-Methanol and its Binaries. Fluid Phase Equilib. 1981, 6, 113-139. (19) Prausnitz, J. M.; Lichtenthaler, R. N.; Gomes de Azevedo, E. Molecular Thermodynamics of Fluid-Phase Equilibria; PrenticeHall: New York, 1999. (20) Hayden, J.; O’Connell, J. A Generalized Method for Predicting Second Virial Coefficients. Ind. Eng. Chem. Process Des. Dev. 1975, 14, 209-216. (21) Prausnitz, J. M.; Anderson, T.; Grens, E.; Eckert, C.; Hsieh, R.; O’Connell, J. Computer Calculations for Multicomponent Vapor-Liquid and Liquid-Liquid Equilibria; Prentice-Hall: New York, 1980. (22) Rackett, H. G. Equation of State for Saturated Liquids. J. Chem. Eng. Data 1970, 15, 514-517. (23) Daubert, T. E.; Danner, R. P. Physical and Thermodynamic Properties of Pure Chemicals. Data Compilation; Taylor & Francis: Bristol, PA, 1989. (24) Gmehling, J.; Onken U.; Arlt, W.; Grenzheuser, P.; Kolbe, B.; Rarey, J.; Weidlich, U. Vapor-Liquid Equilibrium Data Collection. DECHEMA Chemistry Data Series; DECHEMA: Frankfurt, 1977; Vol. I, parts 1-18. (25) Loras, S.; Aucejo, A.; Mun˜oz, R.; de la Torre, J. Isobaric Vapor-Liquid Equilibrium for Binary and Ternary Mixtures of 2-Methyl-2-propanol + Methyl 1,1-dimethylpropyl ether + 2,2,4Trimethylpentane. Fluid Phase Equilib. 2000, accepted for publication. (26) Lenka, M.; Anderko, A. On the Composition-Dependent Interaction Parameters in Equations of State. Chem. Eng. Commun. 1991, 107, 173-188. (27) McDermott, C.; Ellis, S. R. M. A Multicomponent Consistency Test. Chem. Eng. Sci. 1965, 20, 293-296. (28) Wisniak, J. A New Test for the Thermodynamic Consistency of Vapor-Liquid Equilibrium. Ind. Eng. Chem. Res. 1993, 32, 1531-1533. (29) Fonseca, I. M. A.; Lobo L. Q. Error Analysis in Barker’s Method: Extension to Ternary Systems. Fluid Phase Equilib. 1999, 154, 205-211. (30) Wisniak, J.; Apelblat, A.; Segura, H. An Assessment of Thermodynamic Consistency Tests for Vapor-Liquid Equilibrium Data. Phys. Chem. Liq. 1997b, 35, 1-58.
Received for review August 15, 2000 Revised manuscript received January 4, 2001 Accepted January 10, 2001 IE000753U
