The Combinatorial Term for COSMO-Based Activity Coefficient Models

Jan 27, 2011 - Rafael de P. Soares*. Departamento de ... Rua Engenheiro Luis Englert, s/n, Bairro Farroupilha, CEP 90040-040, Porto Alegre, RS, Brazil...
0 downloads 0 Views 781KB Size
Download PDF
ARTICLE pubs.acs.org/IECR

The Combinatorial Term for COSMO-Based Activity Coefficient Models Rafael de P. Soares* Departamento de Engenharia Química, Escola de Engenharia, Universidade Federal do Rio Grande do Sul, Rua Engenheiro Luis Englert, s/n, Bairro Farroupilha, CEP 90040-040, Porto Alegre, RS, Brazil ABSTRACT: COSMO-based activity coefficient models, similarly to the UNIQUAC method, rely on two contributions: combinatorial and residual. The combinatorial term should account for differences in size, shape, and free volume. The residual term should account for electrostatics, dispersion, and hydrogen-bonding effects. The recent literature shows a great effort in improving the residual term, while less attention is given to the combinatorial contribution. This is partly justified by the fact that the residual contribution can be orders of magnitude larger than the combinatorial one. Nevertheless, once the activity coefficient of a substance in mixture is given by the product of these two contributions, both are important. In this work, the combinatorial expressions currently in use in COSMO-based models are reviewed. It is also shown that the model performance can be much improved by assuming a combinatorial contribution similar to that present in the modified UNIFAC model.

’ INTRODUCTION Currently, the most successful predictive models for activity coefficients are group contributions methods such as UNIFAC and modified UNIFAC.1 These methods are an extension of the UNIQUAC model, taking advantage of the fact that the free energy of molecules in solution is (to a considerable degree) additive.2 Despite the undoubted merits of these models, they need large sets of experimental data for the derivation of the subgroup and group interaction parameters and suffer from the inability to distinguish between isomers. An alternative approach is to use COSMO-based methods. COSMO is one of the many quantum mechanical continuum solvation models available today3 and stands for COnductor-like Screening MOdel. COSMO-based models allow the prediction of thermophysical properties (almost) without any experimental data. A few input parameters are necessary and need to be estimated only once for all substances.4,5 COSMO for real solvents (COSMO-RS), introduced in 1995 by Klamt,2 was the first model in this category. The required COSMO calculation has been implemented in many quantum chemistry software packages such as Gaussian6,7 TURBOMOLE,8 MOPAC,9 DMol3,10 and GAMESS.11 In 2002, Lin and Sandler12 proposed a COSMO-RS variation called COSMO-SAC (segment activity coefficient). Later, Grensemann and Gmehling1 developed another modification called COSMO-RS(Ol). It can be expected that these models will become a valuable supplement for the traditional group contribution methods in physical chemistry and chemical engineering. For instance, by examining predictions for the infinite dilution activity coefficient (IDAC), Putnam et al.13 concluded that COSMO-based models are quite close to reaching the point of becoming a practical chemical engineering utility. The authors found that COSMORS can be superior to UNIFAC variations for some aqueous mixtures. Similar conclusions were also obtained by Gerber and Soares14 when using COSMO-SAC in a more recent work. Using COSMO-RS(Ol), Grensemann and Gmehling1 have raised some partly negative remarks. The authors have found r 2011 American Chemical Society

that, in all cases tested, the modified UNIFAC (Do) method provides the most reliable results not only for vapor-liquid equilibria, but also for IDACs and excess enthalpies. Further, they have found that the new method, COSMO-RS(Ol), could not reduce the isomer weakness present in group contribution methods. In spite of some more and other less encouraging results, it is clear that COSMO-based models still need improvement to replace the well-known group activity coefficient methods. In the following section the current combinatorial expressions present in several COSMO-based versions are reviewed. Then the need for a better representation of combinatorial effects is justified. The use of a simple empirically modified Staverman-Guggenheim combinatorial term is suggested. By comparison with experimental IDAC data, it is concluded that the model performance of COSMO-based implementations can be much improved by using the suggested expression.

’ COSMO-BASED MODELS In COSMO-based models, in a very similar way to UNIQUAC, the activity coefficient is computed as the result of two contributions:12 comb ln γi ¼ ln γres i þ ln γi

ð1Þ

It is well-known that the residual contribution can be orders of magnitude larger than the combinatorial one. However, as long as the activity coefficient is given by the product of the contributions, a poor performance in the prediction of the combinatorial part will directly impact the model quality. For descriptions of the currently available formulations for the residual term, check Received: October 13, 2010 Accepted: January 6, 2011 Revised: December 21, 2010 Published: January 27, 2011 3060

dx.doi.org/10.1021/ie102087p | Ind. Eng. Chem. Res. 2011, 50, 3060–3063

Industrial & Engineering Chemistry Research Grensemann and Gmehling,1 Klamt et al.,4 Wang et al.,15and Hsieh et al.16 Combinatorial Term. The combinatorial contribution should take into account differences in size, shape, and free volume. Free-volume effects are usually neglected, because they are more pronounced only for polymer/solvent systems and can only be evaluated if the molar volume of the substances is known.17 Once the COSMO computations need to define a cavity comprising the molecule, differences in size and shape can be considered without any additional information other than that already determined. In the COSMO-RS formulation by Klamt et al.4 a very simple expression for the combinatorial term is suggested:   AS ln γcomb ¼ λ ln ð2Þ i Ai where λ is a parameter estimated by the authors as 0.14, AS is the mean cavity surface area of all components, and Ai is the cavity surface area of substance i. According to Lin and Sandler,12 and also followed by the several COSMO-SAC implementations currently available,14-16 a better fit is possible if the combinatorial term is computed by the Staverman-Guggenheim (SG) equation:     z φi φi q þ 1 ln γcomb ¼ ln φ þ 1 φ ln ð3Þ i i i i 2 θi θi P where P φ i  r i / j r j x j is the normalized volume fraction; θi  q i / j q j x j is the normalized surface-area fraction; z is the coordination number, usually taken as 10; x i is the mole fraction; r i = V i/r and q i = A i/q are the normalized volume and surface area, respectively; A i is the cavity surface area and V i is the cavity volume; q and r are universal parameters of the model. As well pointed out by Gerber and Soares,14 it is easy to check that φi and θi, as defined above, are independent of the parameters q and r. Thus, the model is independent of r and only the quotient z/q is relevant and not the parameters z and q independently. Finally, if the typical value of 10 is assumed for z, eq 3 has q as its single adjustable parameter. Grensemann and Gmehling 1 also use eq 3 for the combinatorial contribution, but with φ i and θi empirically modified by composition-dependent exponents. For θi , the exponent # is given by δ[1 - (q min,i/q max,i )], with P q min,i = min(q i ,q i ), # # q max,i = max(q i ,q i ), q i = (1 - x i ) j6¼i q j /x j , and δ is an adjustable parameter. The authors have shown that the new expression performed better than the original one, when comparing infinite dilution activity coefficients of binary mixtures of n-heptane/alkane. A potential drawback of the proposed exponentials is that they are defined in terms of discontinuous functions of the composition and it is not clear how these empirical modifications will perform for multicomponent mixtures. On the other hand, according to Zhong,18 one of the most important improvements in eq 3 was obtained by replacing the volume fraction φi in the first portion of the equation (the Flory-Huggins term) by X pj p r j xj ð4Þ φi0  r i i xi = j

ARTICLE

Table 1. Resulting Optimized Parameters and Comparison with Literature Values q/Å2 eq 3, parameters from ref 15 eq 5, calibrated in this work

79.53 124

p

ln AAD

-

0.31

0.637

0.036

The resulting expression, which is used in modified UNIFAC implementations with a single value for pi = p, is given by     z φi φi 0 0 comb ln γi ¼ ln φi þ 1 - φi - qi ln þ1ð5Þ 2 θi θi Equation 5 has already proven to be very effective in describing the combinatorial contribution when the molecule surface area and volume are computed by group contribution methods. In this work the performance of eq 5, when the molecule surface area and volume come from the COSMO cavity, is investigated. A component-independent exponent p is considered, likewise adopted in modified UNIFAC. Thus, there are only two parameters to be calibrated: q and p.

’ COMPUTATIONAL ASPECTS All tests in this work were performed by modifying the computational program known as JCOSMO, previously implemented by our group and described elsewhere.14 The tests were performed using a single processor of an Intel Core 2 Quad CPU with 3.2 GB of memory running Ubuntu Linux 10.04. The molecule cavity was determined using MOPAC as described by Gerber and Soares,14 but using RM1, instead of AM1, which tends to improve the performance for organic molecules and biomolecules.19 MOPAC is a semiempirical quantum chemistry program which differs from Gaussian 03, TURBOMOLE, and DMol3 in the level of sophistication. The approximations considered in semiempirical methods enable it to reduce the computing time by orders of magnitude. For instance, MOPAC can perform a COSMO computation with geometry optimization in 0.12 s for WATER. If GAMESS is used, which compares to DMol3 in complexity, for a similar computation using the configuration flags as described by Wang et al.20 plus geometry optimization, the computation time is 61.7 s. It should be noted that, although Gerber and Soares14 have shown that the COSMO-SAC performance can be degraded when based on MOPAC computations, this is mainly due to the residual contribution. For the combinatorial contribution only the cavity is taken into account and not the apparent surface charges. Thus, the quantum chemistry package used is of less importance. ’ RESULTS AND DISCUSSION Experimental Database. A convenient way to investigate the nonenergetic term alone is by considering alkane mixtures of differing size and shape. The apparent surface charges of alkanes determined by COSMO computations are essentially zero. Thus, the residual contribution becomes negligible and the combinatorial term can be investigated alone. In fact, the resolution necessary for the precise determination of enthalpic effects in alkane-alkane systems is usually within the experimental error.21,22 However, paraffin mixtures are not exactly athermal and Kato et al.21 have obtained significant enthalpic effects for these systems. On the other hand, Briard et al.23 have shown 3061

dx.doi.org/10.1021/ie102087p |Ind. Eng. Chem. Res. 2011, 50, 3060–3063

Industrial & Engineering Chemistry Research

ARTICLE

Figure 1. Experimental IDAC logarithm versus model IDAC logarithm, 146 experimental points.

recently by calorimetry that the mixing of C20-C36 n-alkanes in n-heptane can be safely considered an athermal process. None of the currently available COSMO-based formulations is able to reproduce thermal effects in alkane mixtures, and the usual assumption of predominant entropic effects1,12,24 will be considered here. The comparisons were based on the prediction of experimental infinite dilution activity coefficients (IDACs) taken from Parcher et al.,22 which contains IDAC data for C4-C8 linear and branched alkanes in C22-C36 linear alkanes. The database was complemented by some data present in the collections from Kato25 and Kontogeorgis et al.26 From the last one, we have taken data for a large solute (n-octadecane) in small solvents (C6-C7). The resulting database contains a total of 146 IDAC points of linear and branched alkanes ranging from C4 to C36, which assures a variety of dissimilarities in size and shape. Parameter Optimization. In order to calibrate the two parameters q and p of the modified SG, eq 5, the IDAC logarithm of the absolute average deviation (ln AAD) was used: NP 1 X exp ln AAD ¼ jln γi - ln γcalc i j NP i¼1

ð6Þ

where NP is the total number of experimental points. The parameters were optimized to minimize the ln AAD using a direct search method.27 All collected experimental data were considered in the parameter optimization. Once there are only two parameters to be optimized, there is no need for different fit and test sets. With the use of eq 3 and q = 79.53 Å2, as suggested by Wang et al.,15 the ln AAD value was 0.31. The resulting optimized values for eq 5 with q = 124 Å2 and p = 0.637 lead to a ln AAD value as small as 0.036. These results are summarized in Table 1, and the fit improvement can be seen in Figure 1. It should be noted that very similar fits were possible either using atoms radii as suggested by Klamt et al.4 or using multiples around 1.2 of the radii determined by Bondi.28 Of course, different radii led to different paramenters. The values reported above were obtained with radii 18% larger than those determined by Bondi. Model Predictions and Discussion. As mentioned before, it was possible to reduce the ln AAD from 0.31 to 0.036 by using eq 5 and adjusting only two universal parameters. A better visual idea of the improvement is possible with the aid of Figure 1. From Figure 1 it is also possible to see that eq 3 always underestimate the activity coefficient. It is clear that when there is

no size or shape effect (ln γ ≈ 0; the combinatorial contribution is negligible) both versions give reasonable predictions. It is important to note that the deviations given by eq 3 (currently considered in most COSMO-SAC implementations) will impact all systems in the absolute same way, regardless of whether the residual term is more pronounced, as long as there are size or shape effects involved. Regarding the simple relation present in eq 2, used in the COSMO-RS implementation by Klamt et al.,4 we have found a major issue. Let us consider two molecules with different surface areas. From the simple relation in eq 2 it is easy to verify that of two dissimilar molecules will necessarily have the ln γcomb i opposite signs while it is known that they both should be negative. From our tests, we have also found that the λ value, determined by the authors to be 0.14, is suitable only for the activity of large solutes dissolved in small solvents. In order to obtain reasonable values of the activity of a small solute dissolved in a big solvent, the λ value needs to be something around -0.14.

’ CONCLUSIONS COSMO-based models, likewise UNIQUAC and UNIFAC, consider a combinatorial contribution and a residual contribution. The recent literature shows a great effort in improving the residual term, while less attention is given to the combinatorial contribution. This is partly justified by the fact that the residual contribution can be orders of magnitude larger than the combinatorial one. In this work, the need for a better combinatorial term for COSMO-based models was justified. First, the combinatorial models currently in use in COSMO-based models were reviewed. The expression suggested by Klamt et al.4 to be used with their COSMO-RS implementation was found to be suitable only for computing the activity of large solutes dissolved in small solvent molecules. The Staverman-Guggenheim (SG) equation suggested by Lin and Sandler,12 and also followed by the several COSMO-SAC implementations currently available,14-16 always underestimates the activity coefficient. It was also shown that an empirically modified SG term, depending on only two universal parameters, can much improve the model performance. This modified form is identical to that used by modified UNIFAC. A comparison was made in terms of the deviation with respect to experimental data for the logarithm of infinite dilution activity coefficients of linear and branched alkanes ranging from C4 to C36. While the deviation for the SG term currently in use in most COSMO-SAC implementations 3062

dx.doi.org/10.1021/ie102087p |Ind. Eng. Chem. Res. 2011, 50, 3060–3063

Industrial & Engineering Chemistry Research was 0.31, a deviation as low as 0.036 was possible with the modified form.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Tel.: þ55 51 33083528. Fax: þ55 51 33083277.

’ ACKNOWLEDGMENT This work was partially supported by ARD-2010 FAPERGS. ’ LIST OF SYMBOLS δ = adjustable parameter of the exponent given by Grensemann and Gmehling1 ¥ γi = infinite dilution activity coefficient = combinatorial contribution to the activity coefficient γcomb i = residual contribution to the activity coefficient γres i γi = activity coefficient of the substance i in solution λ = adjustable parameter of eq 2 φi = normalized volume fraction θi = normalized surface-area fraction Ai = total cavity surface area of substance i, Å2 AS = mean solution cavity surface area, Å2 p = empirical exponent for modified Staverman-Guggenheim equation q = area normalization parameter, Å2 qi = normalized surface-area parameter r = volume normalization parameter, Å3 ri = normalized volume parameter Vi = total cavity volume of the substance i, Å3 xi = mole fraction of component i in solution z = coordination number ’ REFERENCES (1) 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. (2) 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. (3) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999–3094. (4) Klamt, A.; Jonas, V.; Burger, T.; Lohrenz, J. C. W. Refinement and Parametrization of COSMO-RS. J. Phys. Chem. A 1998, 102, 5074–5085. (5) Klamt, A.; Eckert, F. COSMO-RS: a novel and efficient method for the a priori prediction of thermophysical data of liquids. Fluid Phase Equilib. 2000, 172, 43–72. (6) Truong, T. N.; Stefanovich, E. V. A new method for incorporating solvent effect into the classical, ab initio molecular orbital and density functional theory frameworks for arbitrary shape cavity. Chem. Phys. Lett. 1995, 240, 253–260. (7) Barone, V.; Cossi, M. Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a Conductor Solvent Model. J. Phys. Chem. A 1998, 102, 1995–2001. (8) Schafer, A.; Klamt, A.; Sattel, D.; Lohrenz, J. C. W.; Eckert, F. COSMO Implementation in TURBOMOLE: Extension of an Efficient Quantum Chemical Code towards Liquid Systems. Phys. Chem. Chem. Phys. 2000, 2, 2187–2193. (9) Stewart, J. J. P. MOPAC 2009, Stewart Computational Chemistry, Version 10.251L, http://OpenMOPAC.net. (10) Andzelm, J.; Kolmel, C.; Klamt, A. Incorporation of solvent effects into density functional calculations of molecular energies and geometries. J. Chem. Phys. 1995, 103, 9312–9320.

ARTICLE

(11) Baldridge, K.; Klamt, A. First principles implementation of solvent effects without outlying charge error. J. Chem. Phys. 1997, 106, 6622–6633. (12) 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. (13) Putnam, R.; Taylor, R.; Klamt, A.; Eckert, F.; Schiller, M. Prediction of Infinite Dilution Activity Coefficients Using COSMORS. Ind. Eng. Chem. Res. 2003, 42, 3635–3641. (14) Gerber, R. P.; Soares, R. de P. Prediction of Infinite-Dilution Activity Coefficients Using UNIFAC and COSMO-SAC Variants. Ind. Eng. Chem. Res. 2010, 49, 7488–7496. (15) Wang, S.; Sandler, S. I.; Chen, C.-C. Refinement of COSMOSAC and the Applications. Ind. Eng. Chem. Res. 2007, 46, 7275–7288. (16) Hsieh, C.-M.; Sandler, S. I.; Lin, S.-T. Fluid Phase Equilib. 2010, 297, 90–97. (17) Kouskoumvekaki, I. Fluid Phase Equilib. 2002, 202, 325–335. (18) Zhong, C. Improvement of predictive accuracy of the UNIFAC model for vapor-liquid equilibria of polymer solutions. Fluid Phase Equilib. 1996, 123, 97–106. (19) Rocha, G. B.; Freire, R. O.; Simas, A. M.; Stewart, J. J. P. RM1: a reparameterization of AM1 for H, C, N, O, P, S, F, Cl, Br, and I. J. Comput. Chem. 2006, 27, 1101–1111. (20) 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. (21) Kato, S.; Hoshino, D.; Noritomi, H.; Nagahama, K. Determination of Infinite-Dilution Partial Molar Excess Entropies and Enthalpies from the Infinite-Dilution Activity Coefficient Data of Alkane Solutes Diluted in Longer-Chain-Alkane Solvents. Ind. Eng. Chem. Res. 2003, 42, 4927–4938. (22) Parcher, J. F.; Weiner, P. H.; Hussey, C. L.; Westlake, T. N. Specific retention volumes and limiting activity coefficients of C4-C8 alkane solutes in C22-C36 n-alkane solvents. J. Chem. Eng. Data 1975, 20, 145–151. (23) Briard, A.-J.; Bouroukba, M.; Petitjean, D.; Dirand, M. Variations of Dissolution Enthalpies of Pure n-Alkanes in Heptane as a Function of Carbon Chain Length. J. Chem. Eng. Data 2003, 48, 1574– 1577. (24) Kontogeorgis, G. Improved models for the prediction of activity coefficients in nearly athermal mixtures Part II. A theoreticallybased GE-model based on the van der Waals partition function. Fluid Phase Equilib. 1997, 127, 103–121. (25) Kato, S. Infinite dilution activity coefficients of n-alkane solutes, butane to decane, in nalkane solvents, heptane to hexatriacontane. Fluid Phase Equilib. 2002, 194-197, 641–652. (26) Kontogeorgis, G. M.; Nikolopoulos, G. I.; Fredenslund, A.; Tassios, D. P. Improved models for the prediction of activity coefficients in nearly athermal mixtures Part II. A theoretically-based GE-model based on the van der Waals partition function. Fluid Phase Equilib. 1997, 127, 103–121. (27) Wright, M. H. Direct Search Methods: Once Scorned, Now Respectable. In Numerical Analysis 1995, Proceedings of the 1995 Dundee Biennial Conference in Numerical Analysis; Griffiths, D. F., Watson, G. A., Eds.; Addison Wesley Longman: Harlow, U.K., 1996; pp 191-208. (28) Bondi, A. van der Waals Volumes and Radii. J. Phys. Chem. 1964, 68, 441–451.

3063

dx.doi.org/10.1021/ie102087p |Ind. Eng. Chem. Res. 2011, 50, 3060–3063