1496
Ind. Eng. Chem. Res. 2011, 50, 1496–1503
First-Principles Prediction of Vapor-Liquid-Liquid Equilibrium from the PR+COSMOSAC Equation of State Chieh-Ming Hsieh and Shiang-Tai Lin* Department of Chemical Engineering, National Taiwan UniVersity Taipei, Taiwan 10617, Taiwan
We present the prediction of vapor-liquid-liquid equilibrium (VLLE) from a newly developed PR+COSMOSAC equation of state (EOS). Unlike semiempirical thermodynamic models whose parameters vary with the data used in regression and reparameterization are usually necessary for describing VLLE, this theoretically sound model with its parameters determined from first-principles calculations is not biased against descriptions of either VLE or LLE and can thus be used to predict the occurrence conditions of VLLE. In this model, the molecular interaction and size parameters in the Peng-Robinson EOS are determined directly from a solvation model based on first principle solvation calculations. This new EOS contains neither species dependent parameter nor binary interaction parameters and can be used to predict vapor pressure, liquid density, critical properties of pure substances, and VLE and LLE of mixtures. Without adjustment to the model parameters, this method can predict VLLE with accuracy similar to that of the method modified UNIFAC model. This model is particular useful for the design of new processes involving chemicals whose interaction parameters are not available due to the lack of experimental data. 1. Introduction The knowledge of thermodynamic properties and phase equilibria of pure and mixture fluids is of great importance in the design and optimization of chemical processes.1,2 For example, the information regarding temperature, pressure, and composition in vapor-liquid equilibrium (VLE) and liquid-liquid equilibrium (LLE) is crucial for the design of distillation and extraction processes.3-6 For some special processes, such as extractive distillation, the capability of describing the more complex phase behaviors such as the vapor-liquid-liquid equilibrium (VLLE) becomes necessary.7-11 While the condition at which VLLE occurs can be obtained from experimental measurements, it is often more economical and time-saving to estimate such data using a suitable thermodynamic model. Many thermodynamic models have been proposed and are proven useful for the correlation and/or prediction of the phase behaviors of fluids. Since van der Waals first proposed the twoparameter cubic equations of state (EOS),5 such an EOS (e.g., the SRK12 and the PR EOS13) has been widely used in the industry for process design, simulation, and optimization because of its high accuracy and ease of use. However, these EOS are usually less accurate for associating fluids and liquid mixtures. Some of these issues are address in more recent, none cubic types of EOS such as SAFT,14,15 PC-SAFT,16 and CPA.17 Another class of thermodynamic models focuses on the modeling of nonideality, or the excess free energy (Gex), in the liquid phase. Examples such as the Wilson,18 NRTL,19 UNIQUAC,20 UNIFAC,3 modified UNIFAC,21-23 COSMO-RS,24,25 and COSMO-SAC26,27 model have been implemented in many process simulators28 for the correlation or prediction of low pressure VLE (with assumption of gas phase to be ideal) and LLE of mixtures. Parameters in liquid models usually depend on the experimental data used in the regression. For example, parameters obtained from fitting to VLE experimental data cannot describe LLE well,29-31 or vice versa. For the same reason, the original UNIFAC parameters usually fail to describe LLE, and an additional parameter set, UNIFAC-LLE parameter * To whom correspondence should be addressed. E-mail stlin@ ntu.edu.tw.
table,29 was developed specifically for LLE.32 Therefore, these models should be used with caution as the origin of parameters may affect the design strategy for new processes.1 Some of these problems are removed in the modified UNIFAC model,21-23 whose parameter set was optimized against thousands of experimental data (including both VLE and LLE). Unfortunately, there are still many missing parameters in the parameter matrix for modified UNIFAC33 which severely limit its applications to new processes. A third class is the combination of a cubic equation of state (EOS) with a liquid model through the use of a Gex-based mixing rule.34-39 This approach has the advantages of both the EOS (accurate description for pressure-volume-temperature relation) and liquid models (accurate description for liquid phase nonideality). This kind of approach can successfully describe VLE for a wide range of temperature and pressure. EscobedoAlvarado and Sandler40 and Matsuda et al.41 have demonstrated the prediction of high pressure LLE using parameters determined from low pressure LLE data. Nevertheless, there has been no study regarding the description of LLE and VLE using the same set of parameters. In this work, we present the prediction of VLLE using a recently developed PR+COSMOSAC equation of state. In this method, the energy parameter a(T, x) and volume parameter b(x) in the Peng-Robinson EOS are obtained from first principle solvation calculations. This method has been shown to provide satisfying predictions for a wide variety of thermodynamic properties, VLE, and LLE.42-46 In this work, we show that, with its parameters obtained from first-principle calculations, this method is not biased against the description VLE or LLE and can be more reliable for producing the correct VLLE phase behaviors of mixtures. Furthermore, there is no issue of missing parameters in this method because no species dependent parameter is required. This approach combines the merits from both the EOS, accurate for low density systems, and the quantum-mechanically based solvation model, accurate for liquid systems, and therefore allows for the simultaneously description of VLE, LLE, and VLLE. We believe that this method is ideal
10.1021/ie100781a 2011 American Chemical Society Published on Web 06/16/2010
Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011
for the initial design of new processes, where some or all the needed experimental data are not immediately available. 2. Theory The Peng-Robinson equation of state13 has been proven to be a powerful equation that accurately describes both the vapor and the liquid phase. In this equation, the P-V-T (pressure-volume-temperature) of a fluid are related to one another as follows P)
a RT V _-b V _ (V _ + b) + b(V _ - b)
(1)
where a(T) is a temperature-dependent energy parameter and b is a covolume parameter. Conventionally, these two parameters must be determined from the critical properties (Tc and Pc) and the acentric factor (ω) of the chemical substance. For mixtures, a mixing rule, such as the van der Waals one-fluid mixing rule, Wong-Sandler mixing rule,36,37 or the MHV1 mixing rule,38 is necessary to describe the concentration dependence of these two parameters. Recently, Lin et al.45 and Hsieh and Lin42-44,46 propose a new approach to determine these two parameters, a(T, x) and b(x), from solvation charging free energy. Taking the advantage of advances in quantum mechanical calculations, the solvation free energy can be obtained with input of only the atom connectivity. Therefore this model is, in principle, applicable to any kind of chemical species and their mixtures, especially when no experiment of critical properties, acentric factor, or binary interaction parameter is available. These two parameters are determined from a(T, _x) )
b(x_) ∆G _ *chg(T, _x) CPR
b(x_) )
∑xb
where xi is the mole fraction of substance i in mixture; bi is the volume parameter of species i and equals the molecular cavity volume (multiplied by the Avogadro’s number) in the solvation calculation; CPR is a EOS-dependent constant and found to be -0.623 for the PR EOS;43 ∆G _ *chg is the total solvation charging 45 free energy free energy and is calculated from
∑ x ∆G_ i
*chg x) i/S (T, _
(4)
i
where the subscript i/S denotes solute i dissolved in solvent S. *chg [Equation 4 implies ∆G _ *chg(T) ) ∆G _ i/i (T) for pure fluids.] Once the a(T,x) and b(x) in the PR EOS are determined from eqs 2-4, the fugacity for a species in the fluid can be determined from Appendix I. jf i(T, P, _x) ∆G*chg Pbi i/S Pb ) - ln 1 - ln z ln + xiP RT zRT zRT - Pb (5)
(
*chg (eq 6) + ∆G _ *chg _ *chg _ i/i i/S ) ∆G i/S (eq 6) - ∆G
aiCPR bi
(7)
where ∆G _ *chg _ *chg i/S (eq 6) and ∆G i/i (eq 6) indicate the corresponding charging free energies calculated from eq 6; ai and bi are calculated from Tc, Pc, and ω as in the original PR EOS.13 This approach (eq 7) is referred to as PR+COSMOSAC+TcPcω. It should be noted that the only difference between PR+COSMOSAC (eq 6) and PR+COSMOSAC+TcPcω (eq 7) is the equation used *chg in the determination of ∆G _ i/S . 3. Computational Details
i
∆G _ *chg(T, _x) )
The first three terms on the right of eq 6 are determined from the results of COSMO calculation.48 The last term, i.e., the dispersion contribution, is a semiempirical function of temperature and the exposed surface area of molecule i.44 Details on the calculation of these contributions can be found elsewhere43,44,46 and are briefly summarized in Appendix II. It should be noted that the 4 energy components to the solvation charging free energy are of comparable importance (see Appendix III for the examples of their contributions). The success of this method relies on the accurate modeling of all these free energy components. Since the solvation charging free energy is obtained from the results of COSMO calculation48 and the COSMOSACmodel,26,27 thisapproachisrefereedtoasthePR+COSMOSAC EOS. While the PR+COSMOSAC does not require input of any experimental data, Hsieh and Lin43 showed that the prediction accuracy, in particular the mixture VLE, can be significantly improved when experimental critical properties and acentric factor are used. This is because that the vapor pressures of the pure fluids can be better described. In this case, the solvation charging energy are calculated from
(2) (3)
i i
1497
)
The solvation charging free energy is determined from the summation of four contributions: the ideal solvation (is), charge averaging correction (cc), restoring (res), and dispersion (dsp) contribution,47 ∆G _ *chg x) ) ∆G _ *is _ *cc + ∆G _ *res x) + ∆G _ *dsp (T) i/S (T, _ i + ∆G i i/S (T, _ i (6)
To use the Peng-Robinson equation of state, the energy parameter a(T, x) and volume parameter b(x) of the system are determined first. These two parameters are determined from the results of quantum mechanical (QM) COSMO solvation calculations. This is the most time-consuming step; however, a large database of DFT/COSMO calculation results maintained by Liu’s group at the Virginia Polytechnic Institute and State University (denoted as the VT database49-51) is freely available and is used in this work. Once the QM data are available, the calculation of a(T, x) and b(x) takes a fraction of a second using a modern personal computer. Temperature and composition independent properties (∆Gi*is and ∆Gi*cc) are first determined (eqs A2-1 and A2-3) for each component. The volume parameter bi is determined from the molecular cavity volume in solvation calculation results. For a given temperature, the composition independent dispersion term (∆Gi*dsp) is calculated (eqs A210-A2-12). For a given mixture composition (x), the composition dependent terms (∆Gi/S*res(T, x), ∆Gi/S*chg(T, x), and ∆G*chg(T, x)) can then be obtained (eqs A2-6-A2-9, eq 6 (PR+COSMOSAC) or eq 7 (PR+COSMOSAC+TcPcω), and eq 4). Finally, the interaction parameter a(T, x) is obtained from eq 2 and b(x) from eq 3. Once parameters a(T, x) and b(x) are determined, the vapor-liquid and liquid-liquid equilibrium calculation can be done through the standard procedure.5,6 For isobaric VLLE calculations, the procedure suggested by Iwakabe and Kosuge30 is used. This calculation consists of four steps: (1) For a given pressure P, the VLE over the whole composition range is calculated and the temperatures and compositions at the two turning points are recorded. (The
0.3377 0.2341 0.9709 0.3826 0.9707 0.4955 7.25% 0.9437 0.9854 0.4571 0.9253 0.4638 0.9638 5.64% 331.96 344.69 362.97 347.32 365.87 360.34 0.15% 0.1301 0.1463 0.2873 0.6453 0.3273 0.7599 0.5846 1.73% 0.3754 0.2311 0.2256 0.9843 0.4322 0.9791 0.5566 5.24% a
0.3478 0.7590 0.6201 0.9498 0.6393 0.9590 346.67 365.92 360.47
0.4236 0.9481 0.7019
0.3130 0.2240 0.9880 344.25
Component 1 is water in all eight cases. b AAD ) 1/M ∑ i |xical - xiexp| × 100% and AARD ) 1/M ∑ i |Tical - Tiexp|/Tiexp × 100%.
0.9101 0.9739 0.9895 0.5290 0.9735 0.6585 0.9848 1.65% 327.90 329.77 344.03 363.58 347.19 367.91 361.05 0.25%
y1
0.2883 0.1434 0.8886 342.72
0.2102 0.0777 0.0607 0.1736 0.4140 0.1600 0.5952 0.4513 13.18% 0.3121 0.2718 0.1941 0.1996 0.9839 0.3159 0.6693 0.5619 11.09% 0.9365 0.9467 0.9758 0.9909 0.5212 0.9833 0.9763 0.9832 7.88% 355.14 332.48 327.05 349.83 371.91 348.02 382.15 374.27 2.54% 0.2850 0.1340 0.2260 0.3000 0.9630 0.8980
Tcal (K) y1 x1
x1 Tcal (K) y1 x1
x1 Texp (K)
comp 2
344.05 326.85
xII1 xI1 Tcal (K)
modified UNIFAC
y1 x1
x1
II II
PR+COSMOSAC
I II I a
Table 1. Comparison of Predicted VLLE Point at Atmospheric Pressure from Different Models
This section presents the application of PR+COSMOSAC for the prediction of VLLE. There are a total of 9 binary systems considered in this work, all of which are highly nonideal. The predicted results form PR+COSMOSAC, PR+COSMOSAC+TcPcω, and a widely used group contribution method, the modified UNIFAC model,21 are compared to experimental data. 4.1. Isobaric VLLE at Atmospheric Pressure. The predictions of VLLE temperature and composition in three phases at atmospheric pressure for 8 binary mixtures at are summarized in Table 1. The corresponding phase diagrams of isobaric VLE and LLE are illustrated in Figure 1. As shown in Figure 1, the PR+COSMOSAC, without use of any experimental data, captures the general features of VLLE, although it is least accurate among the three approaches considered in this study. The poor accuracy is a result of its inaccuracy in predicting the pure fluid vapor pressure or normal boiling temperature.43,44 (A comparison of predicted normal boiling temperature and liquid molar volume from various approaches is provided in Appendix IV.) Therefore, when the critical properties and acentric factor are used (eq 7), i.e., the PR+COSMOSAC+TcPcω model,43 the prediction accuracy is significantly improved. The accuracy in VLLE predictions from the modified UNIFAC model (required the use of tens of group interaction parameters and the use of experimental vapor pressures) is similar to that from PR+COSMOSAC+TcPcω. While PR+COSMOSAC gives the largest deviations, it is the only one approach that can successfully predict VLLE point for all eight systems. The modified UNIFAC model fails to describe the LLE (or miscible gap) for water + acrolein [Figure 1f] and PR+COSMOSAC+TcPcω predicts a too low upper critical solution temperature for water + acrylonitrile [Figure 1h]. Furthermore, the turning point (as defined by Iwakabe and Kosuge30) was not found in the VLE calculations for water + acrolein using modified UNIFAC and for water + acrolein using PR+COSMOSAC+TcPcω. This confirms that VLLE cannot be predicted using these two methods. In general, accurate VLLE can be predicted if a model is accurate on both VLE and LLE. This can be seen in the case of water + ethyl acetate using either PR+COSMOSAC+TcPcω or modified UNIFAC [Figure 1e]. Furthermore, good description on LLE helps in getting better liquid phase compositions;
I
PR+COSMOSAC+TcPcω
4. Results and Discussions
acrylonitrile acrolein methyl acetate ethyl acetate 2-methyl-1-propanol 2-butanone 1-butanol 2-butanol AAD/AARDb
turning points are the kinks in the dew point curves in the Pxy diagram.30 If VLLE exists, two turning points can be found in the VLE calculations.) (2) The LLE calculation is then performed at the average temperature T from the two turning points. The compositions (xI, xII) of the two partially miscible liquids are recorded. (3) Given pressure P and the liquid phase composition (xI or xII), the temperature (TI or TII) and vapor phase composition (yI or yII) are determined from isobaric VLE calculations. (4) If the temperature in all these phases and vapor phase composition from the two isobaric calculations are the same (|T - TI| + |T - TII| < 10-3 K and |yI - yII| < 10-6), the VLLE calculation is converged. Otherwise, the average temperature (0.5 × (TI + TII)) is used as the new initial guess for T and steps 2-4 are repeated until convergence. The VLLE calculations usually converge within 10 iterations.
0.1793 0.3248 0.6728 0.3515 0.7504 0.6063 0.82%
Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011
experimental data57,58
1498
Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011
1499
Figure 1. Comparison of VLLE from experiments and predictions for (a) water (1) + 2-methyl-1-propanol (2),59,60 (b) water (1) + 1-butanol (2),57,59 (c) water (1) + 2-butanol (2),57,59 (d) water (1) + methyl acetate (2),59,60 (e) water (1) + ethyl acetate (2),59,60 (f) water (1) + acrolein (2),59,60 (g) water (1) + 2-butanone (2),58,59 (h) water (1) + acrylonitrile (2).59,60 The open squares and triangles are experimental VLE and LLE. The dashed lines, solid lines, and dotted lines are results from PR+COSMOSAC, PR+COSMOSAC+TcPcω, and modified UNIFAC, respectively.
1500
Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011
Table 2. Comparison of Predicted VLLE Point at Different Pressures for Water (1) + Methyl Propionate (2) PR + COSMOSAC
experimental data P (kPa)
Texp (K)
x1I
x1II
y1
Tcal (K)
x1I
35.0 43.7 53.3 65.2 79.2 95.3 113.8 AAD/AARD
318.24 323.24 328.24 333.15 338.15 343.15 348.15
0.1107 0.1162 0.1241 0.1506 0.1675 0.1770 0.1933
0.9890 0.9891 0.9892 0.9893 0.9894 0.9899 0.9891
0.1226 0.2648 0.2701 0.2805 0.2951 0.2733 0.2593
317.83 323.55 328.86 334.44 340.02 345.51 350.95 0.41%
0.1445 0.1546 0.1644 0.1752 0.1866 0.1982 0.2103 2.78%
x1II
PR + COSMOSAC + TcPcω y1
0.9934 0.1101 0.9930 0.1186 0.9927 0.1267 0.9923 0.1354 0.9919 0.1443 0.9915 0.1533 0.9910 0.1623 0.30% 11.64%
whereas good description on VLE helps in getting better equilibrium temperature and gas phase composition. For example, in Figure 1b, PR+COSMOSAC+TcPcω has a better description of LLE, but modified UNIFAC has better description for VLE. 4.2. Isobaric VLLE at Different Pressures. In this section, the VLLE at different pressures for water + methyl propionate are studied and the predicted results are summarized in Table 2 and shown in Figure 2. The accuracy of VLLE predictions are very similar at different pressures. The predicted results are all acceptable from these three approaches. The average errors in the equilibrium temperature and compositions are similar from the modified UNIFAC (vapor pressures are used in VLE calculations) and the PR+COSMOSAC+TcPcω (critical properties and acentric factor are used). The error from PR+COSMOSAC is slightly higher; however, no experimental data is used in any of these calculations. Therefore, PR+COSMOSAC can serve as a complementary approach that provides an acceptable ab initio prediction for VLLE with the molecular structure as the only input.
Modified UNIFAC
Tcal (K)
x1I
x1II
y1
Tcal (K)
x1I
x1II
y1
318.27 323.53 328.39 333.49 338.57 343.57 348.52 0.09%
0.1597 0.1710 0.1821 0.1944 0.2074 0.2208 0.2347 4.72%
0.9930 0.9927 0.9923 0.9919 0.9914 0.9909 0.9904 0.25%
0.2362 0.2492 0.2614 0.2741 0.2868 0.2994 0.3117 3.30%
318.50 323.70 328.52 333.59 338.67 343.68 348.65 0.13%
0.1424 0.1508 0.1589 0.1680 0.1775 0.1874 0.1976 2.05%
0.9894 0.9894 0.9895 0.9895 0.9894 0.9894 0.9893 0.03%
0.2765 0.2877 0.2981 0.3093 0.3205 0.3316 0.3428 5.72%
against a large set of experimental data, the PR+COSMOSAC EOS contains only a few (about 33) nonspecies dependent, universal parameters. There is no issue of missing parameters if a new chemical species is involved. The time-consuming QM calculations has to be done only once for each chemical species and can be stored in a database. Once the database is established (e.g., the VT COSMO database), the time need for phase equilibrium calculations using PR+COSMOSAC is similar to that using modified UNIFAC on a modern personal computer. We consider the PR+COSMOSAC EOS as an ideal complementary method for the development and design of new processes. Acknowledgment This work is dedicated to Prof. Cheng-Ching Yu, who has been a dear friend, and wonderful colleague, and great teacher to the authors. The authors would like to thank the financial support from Grant NSC 98-2221-E-002-087-MY3 by the National Science Council of Taiwan and computation resources from the National Center for High-Performance Computing of Taiwan.
5. Conclusion Combining first-principles solvation calculations and PengRobinson EOS, the PR+COSMOSAC EOS has been shown to be able to provide reasonable prediction on vapor pressures, liquid densities, and critical properties for pure fluids and VLE and LLE for mixtures with a single model and a single set of parameters. Here this method is applied to study complex vapor-liquid-liquid equilibrium of binary mixtures. The occurrences of VLLE in nine studied systems are successfully predicted without use of experimental data or readjustment of any parameter. When the pure fluid critical properties and acentric factor are used, the PR+COSMOSAC+TcPcω provides similar accuracy to the modified UNIFAC model. Unlike the modified UNIFAC whose parameter matrix was optimized
Appendix I: Expression of Fugacity from Solvation Charging Free Energy Lin et al.45 have shown that the fugacity coefficient φi can be *sol determined from the solvation free energy ∆G _ i/S and the compressibility factor z as ln φi ) ln
jf i(T, P, _x) ∆G _ *sol x) i/S (T, P, _ ) - ln z (A1-1) xiP RT
where the solvation free energy is the sum of two contributions: the cavity formation and charging contributions, i.e., *chg ∆G _ *sol _ *cav _ i/S i/S ) ∆G i/S + ∆G
(A1-2)
*cav Hsieh and Lin43,44 showed that ∆G _ i/S can be obtained from
∆G*cav Pbi i/S Pb ) -ln 1 + RT zRT zRT - Pb
(
)
(A1-3)
where b and bi are the volume parameter of mixture and species *sol i, respectively. Then, substituting ∆G _ i/S in eq A1-1 with eqs A1-2 and A1-3, one has Figure 2. Comparison of VLLE from experiments61 and predictions for water (1) + methyl propionate (2) at pressures ranging from 35 to 113.8 kPa (numbers are equilibrium pressure of each point). The open squares are experimental VLLE points. The dashed lines, solid lines, and dotted lines are results from PR+COSMOSAC, PR+COSMOSAC+TcPcω, and modified UNIFAC, respectively.
*chg jf i(T, P, _x) ∆G _ i/S (T, P, _x) Pb ) - ln 1 ln φi ) ln + xiP RT zRT Pbi - ln z (A1-4) zRT - Pb
(
)
Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011 Table A1. Values of Parameters in the PR+COSMOSAC EOS Universal Parameters
n
σm )
parameter
value
CPR aeff (Å2) cHH (kcal · Å4/mol/e2) cAA (kcal · Å4/mol/e2) cOO (kcal · Å4/mol/e2) cHA (kcal · Å4/mol/e2) cHO (kcal · Å4/mol/e2) cAO (kcal · Å4/mol/e2) fpol fdecay Adsp,HB (J/mol) Bdsp,HB (K) Cdsp,HB (K) Adsp,RING (J/mol/K) Bdsp,RING (J/mol)
-0.623 7.5 1757.9468 1121.4047 1757.9468 2462.3206 933.4108 2057.9712 0.6916 3.57 -465876.8150 -429.5556 -141.8436 -0.9181 -365.0667
n
rn2reff2
∑r
2 n
Adsp,i (J/mol/K/Å2)
Bdsp,i (J/mol/Å2)
H C N O F Cl
1.30 2.00 1.83 1.72 1.72 2.05
0.1694 0.1694 0.4045 0.2701 0.1806 0.1566
-191.4602 -191.4602 -207.9411 -178.0767 -125.7842 -201.7754
∑ x n p (σ)
dmn2 rn2 + reff2
)
)
(A2-2)
∆G*res/ RT
∆G*dsp/ RT
∆G*chg/ RT
-8.681 -8.089 -7.316 -6.981 -5.123 -7.155 -5.792 -5.448 -8.336 -6.774
4.210 4.466 4.792 4.847 2.815 4.663 3.124 2.983 -0.416 4.734
-3.909 -4.580 -5.116 -5.587 -5.238 -5.108 -5.105 -5.488 -1.656 -5.481
-8.381 -8.203 -7.640 -7.721 -7.546 -7.600 -7.773 -7.952 -10.407 -7.520
Appendix II: Determination of the Four Components in Solvation Charging Free Energy The ideal solvation contribution accounts for the energy difference of molecule i in ideal conductor and ideal gas phase, i.e.
(A2-1)
where ECOSMO and EIG i i are obtained from the results of quantum mechanical solvation and geometry optimization calculation, respectively. In COSMO-based models,24-27,52 the interactions between molecules are considered as the interaction between the screening charges on the molecular cavity surface. The distribution of these charges on the molecular surface are originally determined from COSMO solvation calculation.48 An important assumption in these models is that the segment interaction only existing between contact surface segments of the same size aeff, so a charge averaging process is needed to eliminate interactions between segments. In this work, the charge averaging equation is
(A2-3)
p(σ) )
i
∑xn
∑ x A p (σ) i i i
i i i
(∆G*is + ∆G*cc)/ RT
COSMO ∆G _ *is - EIG i ) Ei i
(
exp -fdecay
rn + reff2
where fpol ) 0.691647 is the polarization factor. The dielectric energy is defined as Ediel(q_*) ) 1/2∑VφVq*V where φV is the electrostatic potential due to the solute at position V and q*V is the screening charge at some position V on the cavity surface. After the charge averaging process, σ-profile p(σ), the probability of finding a surface segment with screening charge density σ, can be determined from its definition, i.e.
Table A2. Comparison of Free Energy Components in the PR+COSMOSAC EOS
350.45 325.84 330.09 350.26 381.04 352.79 390.88 372.70 373.15 352.80
+ reff2
(
dmn 2
∆G _ *cc ) fpol1/2[Ediel(q_) - Ediel(q_*)] i
Ri (Å)
acrylonitrile acrolein methyl acetate ethyl acetate 2-methyl-1-propanol 2-butanone 1-butanol 2-butanol water methyl propionate
+ reff2
exp -fdecay
1501
where σ and σ* are charge density after and before charge averaging process; σn ) qn/an is the charge density of segment (an is the surface area of the segment m); rn ) (an/π)1/2 is the radius of segment n; reff ) (aeff/π)1/2 is the radius of a standard surface segment; the parameter fdecay (set to 3.57) was introduced to balance the different unit (between Bohr and angstroms);53 dmn is the distance between segments m and n. The charge averaging contribution considers the energy shift due to this charge averaging process,
atom type
Tb/K
2 n
n
Atom Specific Parameters
compound
rn2reff2
∑ σ*r
2
)
i
(A2-4)
∑xA
i i
i i
i
i
where ni, the number of surface segment for molecule i, is the ratio of molecular surface area (Ai) and surface area of standard segment (aeff); pi(σ) is the σ-profile of substance i and specific property of pure substance. It should be noted that different quantum mechanical calculation for COSMO file will generate different σ-profiles and need different sets of parameters.54-56 On the basis of the approach proposed by Hsieh and Lin,46 the σ-profile is separated into four components for the propose of well describing hydrogen-bonding interaction, hydro pi(σ) ) pnhb (σ) + pamino (σ) + pother (σ) i (σ) + pi i i (A2-5)
where nhb, hydro, amine, and other are components of σ-profile collects non-hydrogen-bonding segments, segments on the hydroxyl groups (e.g. water, alcohol, etc.), segments on the amino groups of primary and secondary amines, and segments on O, N, and F atoms which do not connect to hydrogen atom (e.g. nitro, ether, etc.), respectively. Once the σ-profile is generated, the restoring contribution can be calculated from ∆G _ *res i/S ) ni RT
nhb,hydro, amino,other
∑ ∑ p (σ )ln Γ (σ ) s i
s
s m
s S
s m
(A2-6)
σm
where the segment activity coefficient Γ(σ) is obtained from ln ΓSt(σmt) )
{
- ln
nhb,hydro, amino,other
∑∑ s
σn
pSs(σns) exp
[
-∆W(σmt, σns) + RT ln ΓSt(σns) RT
}
]
(A2-7)
1502
Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011
Table A3. Comparison of Normal Boiling Temperatures and Corresponding Liquid Molar Volumes from Experiments and Predictions experimenta
PR+COSMOSAC
PR+COSMOSAC+TcPcω
b
method name
Tb (K)
VL (m3/kmol)
Tb(K)
VL(m3/kmol)
Tb(K)
VL(m3/kmol)
water acrylonitirle acrolein methyl acetate ethyl acetate 2-methyl-1-propanol 2-butanone 1-butanol 2-butanol methyl propionate overall deviation
373.15 350.45 325.84 330.09 350.26 381.04 352.75 390.81 372.65 352.6
1.81 × 10-2 7.18 × 10-2 6.98 × 10-2 8.37 × 10-2 1.06 × 10-1 1.03 × 10-1 9.79 × 10-2 1.03 × 10-1 1.02 × 10-1 1.05 × 10-1
397.38 357.65 330.74 325.65 352.13 374.71 347.91 389.59 377.29 348.92 1.76%
1.83 × 10-2 5.93 × 10-2 5.99 × 10-2 7.25 × 10-2 8.88 × 10-2 8.53 × 10-2 8.03 × 10-2 8.51 × 10-2 8.51 × 10-2 8.89 × 10-2 14.71%
375 348.93 327.14 330.76 350.56 380.42 352.7 390.11 372 352.48 0.22%
2.25 × 10-2 1.20 × 10-1 8.06 × 10-2 8.74 × 10-2 1.12 × 10-1 1.03 × 10-1 1.05 × 10-1 1.03 × 10-1 1.04 × 10-1 1.08 × 10-1 12.94%
a The experimental normal boiling temperatures are taken from VLE data and liquid molar volumes are from the DIPPR database. properties, the PR+COSMOSAC+TcPcω predicts the same results as those from the original Peng-Robinson EOS.
where the superscript t and s can be nhb, hydro, amino, or other. The segment exchange energy ∆W is determined from ∆W(σmt, σsn) ) fpol
{
0.3aeff3/2 t (σm + σsn)2 - chb(σmt, σsn)(σmt - σsn)2 2ε0 (A2-8)
where ε0, the permittivity of vacuum, is 8.8542 × 10-12 (A2 · s4/ kg · m3); chb(σmt , σsn) is independent of temperature and is cHH if s ) t ) hydro and σmt · σsn < 0 cAA if s ) t ) amino and σmt · σsn < 0 cOO if s ) t ) other and σmt · σsn < 0 t s chb(σm, σn) ) cHA if s ) hydro, t ) amino, and σmt · σsn < 0 cHO if s ) hydro, t ) other, and σmt · σsn < 0 cAO if s ) amino, t ) other, and σmt · σsn < 0 0 otherwise (A2-9)
The dispersion contribution is a function of temperature and exposed surface area of molecule i, ∆G*dsp (T) ) i
∑ S (A k
dsp,kT
*dsp *dsp + Bdsp,k) + GHB (T) + GRING (T)
k
(A2-10) where Sk is the total exposed surface area of atom type k, Adsp,k and Bdsp,k are the dispersion parameters of atom type k, and *dsp *dsp GHB and GRING are the empirical corrections for hydrogenbonding and ring containing (cyclic or aromatic) molecules. The expressions of these two terms are *dsp GHB (T) )
{
}
Adsp,HB 1 NHBH 1 + exp[-(T - Bdsp,HB)/Cdsp,HB] (A2-11)
*dsp GRING (T) ) NAR(Adsp,RINGT + Bdsp,RING)
(A2-12)
where NHBH and NAR are the number of hydrogen-bonding donors (the hydrogen atoms connecting to either nitrogen, oxygen, or fluorine atom) and the number of atoms involved in any ring-structure, respectively. The values of parameters are directly taken from previous work46 and listed in Table A1. Appendix III: Comparison of the Components of Solvation Charging Free Energies Table A2 lists the solvation charging free energy and its components (ideal solvation plus charge averaging correction,
b
For pure fluid
restoring, and dispersion for the chemical species considered in this work. The charging free energy is determined at the normal boiling temperature of each chemical species. It can be seen that the free energy components are of the same order of magnitude. Therefore, accurate prediction of properties requires an accurate modeling of all free energy components. Appendix IV: Comparison of Experimental and Predicted Properties for Pure Substances The normal boiling temperature can be determined by solving the temperature at which the vapor pressure is 1 atm () 101325 Pa). The comparison of experimental and predicted normal boiling temperatures and the corresponding liquid molar volumes for all ten substances considered in this work are summarized in Table A3. Literature Cited (1) Cadoret, L.; Yu, C. C.; Huang, H. P.; Lee, M. J. Effects of Physical Properties Estimation on Process Design: A Case Study Using Aspenplus. Asia-Pac. J. Chem. Eng. 2009, 4, 729–734. (2) Wang, S. J.; Yu, C. C.; Huang, H. P. Plant-Wide Design and Control of Dmc Synthesis Process Via Reactive Distillation and Thermally Coupled Extractive Distillation. Comput. Chem. Eng. 2010, 34, 361–373. (3) Fredenslund, A.; Gmehling, J.; Michelsen, M. L.; Rasmussen, P.; Prausnitz, J. M. Computerized Design of Multicomponent DistillationColumns Using Unifac Group Contribution Method for Calculation of Activity-Coefficients. Ind. Eng. Chem. Process Design DeV. 1977, 16, 450– 462. (4) Banerjee, T.; Sahoo, R. K.; Rath, S. S.; Kumar, R.; Khanna, A. Multicomponent Liquid-Liquid Equilibria Prediction for Aromatic Extraction Systems Using Cosmo-Rs. Ind. Eng. Chem. Res. 2007, 46, 1292–1304. (5) Sandler, S. I. Chemical and Engineering Thermodynamics, 3rd ed.; John Weily & Sons: New York, 1999. (6) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria. 3rd ed.; Pearson Education Taiwan Ltd.: Taipei, 2004. (7) Higler, A.; Chande, R.; Taylor, R.; Baur, R.; Krishna, R. Nonequilibrium Modeling of Three-Phase Distillation. Comput. Chem. Eng. 2004, 28, 2021–2036. (8) Steyer, F.; Qi, Z. W.; Sundmacher, K. Synthesis of Cylohexanol by Three-Phase Reactive Distillation: Influence of Kinetics on Phase Equilibria. Chem. Eng. Sci. 2002, 57, 1511–1520. (9) Gumus, Z. H.; Ciric, A. R. Reactive Distillation Column Design with Vapor/Liquid/Liquid Equilibria. Comput. Chem. Eng. 1997, 21, S983– S988. (10) Bekiaris, N.; Morari, M. Multiple Steady States in Distillation: Infinity/Infinity Predictions, Extensions, and Implications for Design, Synthesis, and Simulation. Ind. Eng. Chem. Res. 1996, 35, 4264–4280. (11) Bekiaris, N.; Meski, G. A.; Morari, M. Multiple Steady States in Heterogeneous Azeotropic Distillation. Ind. Eng. Chem. Res. 1996, 35, 207– 227. (12) Soave, G. Equilibrium Constants from a Modified Redlich-Kwong Equation of State. Chem. Eng. Sci. 1972, 27, 1197–1203.
Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011 (13) Peng, D.; Robinson, D. B. New 2-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59–64. (14) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. New Reference Equation of State for Association Liquids. Ind. Eng. Chem. Res. 1990, 29, 1709–1721. (15) Chapman, W. G.; Jackson, G.; Gubbins, K. E. Phase-Equilibria of Associating Fluids Chain Molecules with Multiple Bonding Sites. Mol. Phys. 1988, 65, 1057–1079. (16) Gross, J.; Sadowski, G. Perturbed-Chain Saft: An Equation of State Based on a Perturbation Theory for Chain Molecules. Ind. Eng. Chem. Res. 2001, 40, 1244–1260. (17) Kontogeorgis, G. M.; Voutsas, E. C.; Yakoumis, I. V.; Tassios, D. P. An Equation of State for Associating Fluids. Ind. Eng. Chem. Res. 1996, 35, 4310–4318. (18) Wilson, G. M. Vapor-Liquid Equilibrium 0.11. New Expression for Excess Free Energy of Mixing. J. Am. Chem. Soc. 1964, 86, 127–130. (19) Renon, H.; Prausnitz, J. M. Local Compositions in Thermodynamic Excess Functions for Liquid Mixtures. AIChE J. 1968, 14, 135–144. (20) Abrams, D. S.; Prausnitz, J. M. Statistical Thermodynamics of Liquid-Mixtures - New Expression for Excess Gibbs Energy of Partly or Completely Miscible Systems. AIChE J. 1975, 21, 116–128. (21) Gmehling, J.; Lohmann, J.; Jakob, A.; Li, J. D.; Joh, R. A Modified Unifac (Dortmund) Model. 3. Revision and Extension. Ind. Eng. Chem. Res. 1998, 37, 4876–4882. (22) Gmehling, J.; Li, J. D.; Schiller, M. A Modified Unifac Model 0.2. Present Parameter Matrix and Results for Different Thermodynamic Properties. Ind. Eng. Chem. Res. 1993, 32, 178–193. (23) Weidlich, U.; Gmehling, J. A Modified Unifac Model 0.1. Prediction of Vle, He, and Gamma-Infinity. Ind. Eng. Chem. Res. 1987, 26, 1372– 1381. (24) Klamt, A. Conductor-Like Screening Model for Real Solvents - a New Approach to the Quantitative Calculation of Solvation Phenomena. J. Phys. Chem. 1995, 99, 2224–2235. (25) Klamt, A.; Jonas, V.; Burger, T.; Lohrenz, J. C. W. Refinement and Parametrization of Cosmo-Rs. J. Phys. Chem. A 1998, 102, 5074–5085. (26) Lin, S. T.; Sandler, S. I. A Priori Phase Equilibrium Prediction from a Segment Contribution Solvation Model. Ind. Eng. Chem. Res. 2002, 41, 899–913. (27) Wang, S.; Sandler, S. I.; Chen, C. C. Refinement of Cosmo-Sac and the Applications. Ind. Eng. Chem. Res. 2007, 46, 7275–7288. (28) Aspen Plus Unit Operation Models, Aspen Technology Inc: Burlington, MA, 2007. (29) Magnussen, T.; Rasmussen, P.; Fredenslund, A. Unifac Parameter Table for Prediction of Liquid-Liquid Equilibria. Ind. Eng. Chem. Process Des. DeV. 1981, 20, 331–339. (30) Iwakabe, K.; Kosuge, H. A Correlation Method for Isobaric VaporLiquid and Vapor-Liquid-Liquid Equilibria Data of Binary Systems. Fluid Phase Equilib. 2008, 266, 202–210. (31) Skjoldjorgensen, S.; Rasmussen, P.; Fredenslund, A. On the Concentration-Dependence of the Uniquac Unifac Models. Chem. Eng. Sci. 1982, 37, 99–111. (32) Gupte, P. A.; Danner, R. P. Prediction of Liquid-Liquid Equilibria with Unifac - a Critical Evaluation. Ind. Eng. Chem. Res. 1987, 26, 2036– 2042. (33) Jakob, A.; Grensemann, H.; Lohmann, J.; Gmehling, J. Further Development of Modified Unifac (Dortmund): Revision and Extension 5. Ind. Eng. Chem. Res. 2006, 45, 7924–7933. (34) Vidal, J. Mixing Rules and Excess Properties in Cubic Equations of State. Chem. Eng. Sci. 1978, 33, 787–791. (35) Huron, M. J.; Vidal, J. New Mixing Rules in Simple Equations of State for Representing Vapor-Liquid-Equilibria of Strongly Non-Ideal Mixtures. Fluid Phase Equilib. 1979, 3, 255–271. (36) Wong, D. S. H.; Orbey, H.; Sandler, S. I. Equation of State Mixing Rule for Nonideal Mixtures Using Available Activity-Coefficient Model Parameters and That Allows Extrapolation over Large Ranges of Temperature and Pressure. Ind. Eng. Chem. Res. 1992, 31, 2033–2039. (37) Wong, D. S. H.; Sandler, S. I. A Theoretically Correct Mixing Rule for Cubic Equations of State. AIChE J. 1992, 38, 671–680. (38) Michelsen, M. L. A Modified Huron-Vidal Mixing Rule for Cubic Equations of State. Fluid Phase Equilib. 1990, 60, 213–219. (39) Boukouvalas, C.; Spiliotis, N.; Coutsikos, P.; Tzouvaras, N.; Tassios, D. Prediction of Vapor-Liquid-Equilibrium with the Lcvm Model - a Linear Combination of the Vidal and Michelsen Mixing Rules
1503
Coupled with the Original Unifac and the T-Mpr Equation of State. Fluid Phase Equilib. 1994, 92, 75–106. (40) Escobedo-Alvarado, G. N.; Sandler, S. I. Study of Eos-G(Ex) Mixing Rules for Liquid-Liquid Equilibria. AIChE J. 1998, 44, 1178–1187. (41) Matsuda, H.; Kurihara, K.; Ochi, K.; Kojima, K. Prediction of Liquid-Liquid Equilibria at High Pressure for Binary Systems Using EosG(E) Models: Methanol Plus Hydrocarbon Systems. Fluid Phase Equilib. 2002, 203, 269–284. (42) Hsieh, C. M.; Lin, S. T. Prediction of 1-Octanol-Water Partition Coefficient and Infinite Dilution Activity Coefficient in Water from the Pr Plus Cosmosac Model. Fluid Phase Equilib. 2009, 285, 8–14. (43) Hsieh, C. M.; Lin, S. T. First Principle Predictions of Vapor-Liquid Equilibria for Pure and Mixture Fluids from the Combined Use of Cubic Equations of State and Solvation Calculations. Ind. Eng. Chem. Res. 2009, 48, 3197–3205. (44) Hsieh, C. M.; Lin, S. T. Determination of Cubic Equation of State Parameters for Pure Fluids from First Principle Solvation Calculations. AIChE J. 2008, 54, 2174–2181. (45) Lin, S. T.; Hsieh, C. M.; Lee, M. T. Solvation and Chemical Engineering Thermodynamics. J. Chin. Inst. Chem. Eng. 2007, 38, 467– 476. (46) Hsieh, C. M.; Lin, S. T. Prediction of Liquid-Liquid Equilibrium from the Peng-Robinson+Cosmosac Equation of State. Chem. Eng. Sci. 2010, 65, 1955–1963. (47) Lin, S. T.; Chang, J.; Wang, S.; Goddard, W. A.; Sandler, S. I. Prediction of Vapor Pressures and Enthalpies of Vaporization Using a Cosmo Solvation Model. J. Phys. Chem. A 2004, 108, 7429–7439. (48) Klamt, A.; Schuurmann, G. Cosmo - a New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Energy and Its Gradient. J. Chem. Soc.-Perkin Trans. 2 1993, 799–805. (49) St Cholakov, G.; Stateva, R. P.; Shacham, M.; Brauner, N. Prediction of Properties in Homologous Series with a Shortcut Qs2pr Method. Aiche J. 2007, 53, 150–159. (50) Mullins, E.; Oldland, R.; Liu, Y. A.; Wang, S.; Sandler, S. I.; Chen, C. C.; Zwolak, M.; Seavey, K. C. Sigma-Profile Database for Using CosmoBased Thermodynamic Methods. Ind. Eng. Chem. Res. 2006, 45, 4389– 4415. (51) Mullins, E.; Liu, Y. A.; Ghaderi, A.; Fast, S. D. Sigma Profile Database for Predicting Solid Solubility in Pure and Mixed Solvent Mixtures for Organic Pharmacological Compounds with Cosmo-Based Thermodynamic Methods. Ind. Eng. Chem. Res. 2008, 47, 1707–1725. (52) Grensemann, H.; Gmehling, J. Performance of a Conductor-Like Screening Model for Real Solvents Model in Comparison to Classical Group Contribution Methods. Ind. Eng. Chem. Res. 2005, 44, 1610–1624. (53) Lin, S. T.; Sandler, S. I. A Priori Phase Equilibrium Prediction from a Segment Contribution Solvation Model (Vol 41, Pg 903, 2004). Ind. Eng. Chem. Res. 2004, 43, 1322–1322. (54) Wang, S.; Stubbs, J. M.; Siepmann, J. I.; Sandler, S. I. Effects of Conformational Distributions on Sigma Profiles in Cosmo Theories. J. Phys. Chem. A 2005, 109, 11285–11294. (55) Mu, T. C.; Rarey, J.; Gmehling, J. Performance of Cosmo-Rs with Sigma Profiles from Different Model Chemistries. Ind. Eng. Chem. Res. 2007, 46, 6612–6629. (56) Wang, S.; Lin, S. T.; Watanasiri, S.; Chen, C. C. Use of Gamess/ Cosmo Program in Support of Cosmo-Sac Model Applications in Phase Equilibrium Prediction Calculations. Fluid Phase Equilib. 2009, 276, 37– 45. (57) Iwakabe, K.; Kosuge, H. Isobaric Vapor-Liquid-Liquid Equilibria with a Newly Developed Still. Fluid Phase Equilib. 2001, 192, 171–186. (58) Moon, H. M.; Ochi, K.; Kojima, K. Vapor-Liquid-Equilibria for the Ethyl Methyl Ketone + Water - System with Limited Miscibility. J. Chem. Eng. Data 1995, 40, 468–471. (59) Sørensen, J. M.; Arlt, W. Liquid-Liquid Equilibrium Data Collection; Dechema: Frankfurt, 1979; Vol. 1. (60) Gmehling, J.; Onken, U.; Arlt, W. Vapor-Liquid Equilibrium Data Collection; Dechema: Frankfurt, 1977. (61) Hsieh, C. T.; Lee, M. J.; Lin, H. M. Vapor-Liquid -Liquid Equilibria for Aqueous Systems with Methyl Acetate, Methyl Propionate, and Methanol. Ind. Eng. Chem. Res. 2008, 47, 7927–7933.
ReceiVed for reView April 1, 2010 ReVised manuscript receiVed May 29, 2010 Accepted June 3, 2010 IE100781A
