Using Molecular Orbital Calculations To Describe the Phase Behavior

Ind. Eng. Chem. Res. 1998, 37, 2917-2928

2917

Using Molecular Orbital Calculations To Describe the Phase Behavior of Cross-associating Mixtures Jeffrey P. Wolbach and Stanley I. Sandler* Center for Molecular and Engineering Thermodynamics, Department of Chemical Engineering, University of Delaware, Newark, Delaware 19716

In previous studies, we have used molecular orbital calculations to determine the thermodynamic changes of dimerization for a number of pure and cross-associated species. We have shown that by using these results in a physical equation of state, the statistical associating fluid theory (SAFT), we are able to accurately model the phase behavior of pure self-associating compounds and binary mixtures of a self-associating compound and a nonassociating compound with a reduction in the number of adjustable parameters. In this study, we consider the phase behavior of binary mixtures in which cross-associated species may occur. To determine the equation of state parameters describing cross association, we introduce a mixing rule based on the results of our molecular orbital calculations. We show that using information derived from our quantummechanical calculations results in correlations of mixture vapor-liquid equilibrium data with fewer adjustable parameters and no loss of accuracy, indeed frequently with improved accuracy, compared to the original SAFT model. Introduction Molecular association dramatically affects the phase behavior of pure fluids and mixtures. Equation of state models that perform well for associating fluids and their mixtures must explicitly account for the association that occurs. These equations can be separated into two classes: chemical equations that postulate the formation of distinct molecular aggregates in solution and physical equations that model the association as being the result of a strong, specific molecular interaction. Both classes of equations suffer from the use of a large number of ill-defined adjustable parameters in order to obtain good fits of experimental phase behavior data. If one were able to determine the equation of state parameters describing the association a priori, such models would be more useful, and possibly even completely predictive. In previous studies we used ab initio molecular orbital calculations to determine the thermodynamic parameters of association for pure hydrogen-bonded dimers of small organic compounds and have compared the results of our calculations to experimental estimates of the thermodynamic parameters where they are available (Wolbach and Sandler, 1997a,b). We then used the results of our calculations to determine parameters in a physical equation of state, the statistical associating fluid theory (SAFT). Further, we showed that using the information obtained from molecular orbital calculations allows a reduction in the number of adjustable parameters in the SAFT equation of state for pure associating compounds (Wolbach and Sandler, 1997c) as well as for binary mixtures of an associating compound and a diluent (Wolbach and Sandler, 1997d) with little or no loss in accuracy. In this work we are concerned with mixtures in which cross-associated species may form. We consider first binary mixtures of self-associating compounds, and then * To whom correspondence should be addressed. E-mail: [email protected].

mixtures of a self-associating compound and a non-selfassociating compound in which cross-associated species may form. An example of the second type of system is a mixture of methanol and acetone. We again compare the results obtained using our quantum-mechanically determined parameters with those using the original SAFT model. Description of the Molecular Orbital Calculations We have performed ab initio molecular orbital calculations to determine the enthalpy, entropy, Gibbs free energy, and heat capacity changes for a number of hydrogen-bonding dimerization “reactions”. These calculations are for the association of isolated pairs of molecules and are representative of vapor-phase dimerization. The calculations were performed using the Gaussian suite of computational chemistry programs (Gaussian, Inc., 1993, 1995), and the results were converted to thermochemical properties with standard methods of statistical mechanics (McQuarrie, 1976). We have used two different calculational methods: the computationally inexpensive Hartree-Fock (H-F) method and the more rigorous density functional theory method. The H-F calculations were done using the small 6-31g(d,p) basis set, while the density functional theory calculations used the Becke3LYP (B3LYP) (Becke, 1993) functional and the larger 6-31++g(2d,p) basis set. Details of these calculations and a further discussion of these two calculational methods were presented earlier (Wolbach and Sandler, 1997a). As the B3LYP calculations involve a more rigorous calculational procedure, one would expect the results to be more accurate. One goal of this work was to determine whether it was necessary to go to this level of computational intensity and complexity to obtain parameters leading to accurate correlations of mixture vapor-liquid equilibrium (VLE) data, or if the less rigorous H-F calculations were sufficient. In previous studies we have found that the results of the H-F

S0888-5885(97)00781-1 CCC: $15.00 © 1998 American Chemical Society Published on Web 06/02/1998

2918 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 Table 1. Thermodynamic Parameters of Dimerization Determined by Molecular Orbital Calculations H-F

B3LYP

compound

∆H (kcal/mol)

∆S (cal/(mol‚K))

∆CP (298 K) (cal/(mol‚K))

Keq (298 K)

∆H (kcal/mol)


∆CP0 (298 K) (cal/(mol‚K))

Keq (298 K)

water methanol formic acid acetic acid propionic acid

-3.929 -3.877 -13.557 -13.924 -13.720

-18.324 -19.937 -34.515 -34.052 -34.311

2.48 3.07 2.22 1.26 3.44

7.48 × 10-2 3.05 × 10-2 248 582 362

-2.993 -3.130 -13.497 -13.958 -13.742

-21.087 -22.379 -37.290 -36.097 -37.205

1.97 2.75 1.39 1.46 1.50

3.85 × 10-3 2.53 × 10-3 55.5 220 87.5

0

calculations were sufficient to accurately model pure hydrogen-bonding fluids (Wolbach and Sandler, 1997c) and binary mixtures of an associating fluid and a diluent (Wolbach and Sandler, 1997d). The results of our calculations for dimers of the selfassociating compounds considered in this study are presented in Table 1. The methanol and water dimers are linear, while the carboxylic acid dimers are cyclic, as these compounds have been experimentally observed to form only such dimers in the vapor phase. Previously, we have shown that both calculational methods led to results that were in agreement with data derived from experiment. Although the results for the methanol and water dimers are quite different from the two calculational methods, there are also large discrepancies in the reported experimental estimates for these weakly associating compounds, so it is not possible to decide which of the calculational methods is more accurate. In a previous study (Wolbach and Sandler, 1997b) we also considered cross dimers between self-associating compounds. We concluded from the results that the enthalpy and entropy changes upon association for a cross dimer of two self-associating compounds could be approximated by the algebraic average of the enthalpy and entropy changes of the two self-dimers. That is, if compounds i and j can both self-associate, for the crossassociation reaction i + j f ij

for the properties of associating and nonassociating fluids over a wide range of molecular sizes (Huang and Radosz, 1990; Economou and Donohue, 1991, 1992). In this model pure fluids are treated as chains of equal-sized spherical segments interacting with a squarewell potential. The equation of state is expressed in terms of the residual Helmholtz free energy and contains three parameters for nonassociating compounds: the number of segments (m), the volume of an individual segment (v∞), and the depth of the square-well potential (u0/k). To describe association, molecules are assigned bonding sites, and the interactions between these sites are modeled using the square-well potential. This introduces two additional parameters for pure associating fluids. The well depth of the interaction between sites A and B is described by the parameter AB/k and the well width characterized by the dimensionless parameter κAB. In the SAFT model, these association parameters are for the interaction between bonding sites on an isolated pair of molecules. The contribution to the residual Helmholtz free energy due to association is

1 ∆Hij ) (∆Hii + ∆Hjj) 2

(1)

1 ∆Sij ) (∆Sii + ∆Sjj) 2

where M is the number of bonding sites assigned to a molecule, XA is the mole fraction of the original molecules that are not bonded at site A, and the summation is over all the bonding sites of a molecule. The fractions of unbonded sites are determined from

(2)

and

aassoc RT

)

∑ A

[

ln XA -

XA ) [1 +

This implies that

1 ∆Gij ) (∆Gii + ∆Gjj) 2

(3)

and

Kij ) xKiiKjj

(4)

Using eq 4 one obtains estimates of the equilibrium constants for the cross-association reactions without performing additional quantum mechanical calculations. Application to Pure Components To describe the phase behavior of pure associating compounds, we initially used our calculated results in the statistical associating fluid theory (SAFT) developed by Chapman and co-workers (Chapman et al., 1988, 1989, 1990) from the first-order thermodynamic perturbation theory (TPT1) of Wertheim (1984a,b, 1986a,b). SAFT has been shown to result in accurate correlations

]

XA

1 + M 2 2

FXB∆AB]-1 ∑ B

(5)

(6)

The quantity ∆AB is a measure of the strength of association between sites A and B and is approximated (Chapman et al., 1988) as

[ ( ) ]

∆AB ) g(d)seg exp

AB - 1 (x2v∞κAB) kT

(7)

Here g(d)seg is the segmental radial distribution function, taken to be the hard sphere radial distribution function, evaluated at contact. Before using eq 5 to calculate the energy due to association, one must specify the allowed site-site interactions and the values of ∆AB needed to determine the unbonded fractions XA. The resulting expressions can be simplified by assuming the equivalence of specific bonds. Table 2 presents the association schemes used in this work, including bonding approximations and expressions for XA. The association scheme nomenclature used is that of Huang and Radosz (1990). Note that the quantities obtained from molecular orbital calculations and those in the SAFT equation are

Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 2919 Table 2. Bonding Schemes Used in This Work

3B

4C

*0

XA

example A

-1 + x1 + 4F∆AB 2F∆AB XA ) XB

∆AA ) ∆BB ) ∆CC ) 0 ∆AB ) 0 ∆AC ) ∆BC * 0

XA ) XB XC ) 2XA - 1

∆AA ) ∆BB ) ∆CC ) ∆DD ) 0 ∆AB ) ∆CD ) 0 ∆AC ) ∆AD ) ∆BC ) ∆BD * 0

XA ) XB ) XC ) XD

C

H

R

O

A

-1 + x1 + 4F∆AB 2F∆AB

•• O

H

B

R

-(1 - F∆AB) + x(1 + F∆AB)2 + 4F∆AB

A

4F∆AB

R

•• O

B

H

D

C

-1 + x1 + 8F∆AB 4F∆AB

C

••

∆AA ) ∆BB ) 0 ∆AB * 0

O

••

2B

XA approximations

••

1A

∆AA

••

∆ assumptions

association scheme

O

H A

H B

Table 3. Average Errors and Best-Fit Parameters for Pure Components method

association model

v∞ (cm3/mol)

m

u0/k (K)

AB/k (K)

κAB × 103

% error Pvap

% error Fliq

formic acid

original SAFT H-F B3LYP

1A 1A 1A

15.50 18.20 19.75

1.341 1.176 1.099

333.28 369.86 394.33

7522 4160 3771

1.625 582.3 464.5

0.62 0.76 0.62

0.49 0.27 0.10

acetic acid


1A 1A 1A

14.50 18.20 19.75

2.132 1.684 1.582

290.73 305.97 321.28

3941 5751 5490

39.26 3.465 4.157

1.90 0.86 0.94

0.65 1.58 0.63

propionic acid


1A 1A 1A

13.50 18.20 19.75

3.084 2.192 2.065

296.03 282.51 306.75

3400 6123 5395

8.193 0.855 1.715

0.35 0.87 0.75

0.09 0.67 0.87

butyric acid


1A 1A 1A

13.00 18.20 19.75

3.800 2.700 2.548

268.93 272.87 300.23

4155 6217 5604

3.700 0.580 0.728

0.20 0.76 0.69

0.85 1.19 0.72

1-pentanoic acid


1A 1A 1A

12.00 18.20 19.75

4.719 3.208 3.031

248.63 268.97 296.65

4322 6165 5593

5.103 0.631 0.649

0.32 0.85 0.68

0.98 1.25 0.40

water


3B 4C 4C

10.00 8.00 12.55

1.179 1.406 1.000

528.17 212.91 546.63

1810 1809 1237

15.93 91.09 22.03

1.22 1.74 1.94

3.15 2.89 3.69

methanol


2B 3B 3B

12.00 8.10 9.00

1.776 2.679 2.535

216.13 222.57 288.60

2714 1930 1470

48.56 44.72 16.41

0.83 1.23 6.88

0.88 0.64 0.66

compound

different. In an earlier study (Wolbach and Sandler, 1997c), we presented a procedure by which we related the results of our molecular orbital calculations to the association parameters in SAFT. We found that the relationship

∆AB ) C′RTK

densities obtained using both sets of quantum mechanics calculations, along with the errors from using the best-fit parameters of the original SAFT model, are presented in Table 3. Application to Mixtures

(8)

was true regardless of the association model assumed. In eq 8, K is the equilibrium constant for self-association, and the value of the constant C′ was dependent on the association model assumed. For the 1A model C′ was equal to 2, for the 2B model C′ was equal to 1, for the 3B model C′ was equal to 1/2, and for the 4C model C′ was equal to 1/4. The values for C′ were determined by matching the vapor-phase composition of monomers and dimers at low pressures predicted by the SAFT equation to the composition predicted by our molecular orbital calculations. Using eqs 7 and 8 we were able to relate the results of our molecular orbital calculations to the SAFT parameters for self-association and fit pure-component vapor pressure and liquid density data for associating compounds using fewer adjustable parameters than those used in the original SAFT model. The average errors in the calculated vapor pressures and liquid

Equation of State for Mixtures. As was the case for pure fluids, the SAFT equation of state for mixtures is defined in terms of the residual Helmholtz free energy:

ares ) ahs + adisp + achain + aassoc

(9)

Here, we will detail only two of the terms. The contribution to the residual mixture Helmholtz energy due to dispersion is calculated from

adisp 0 adisp )m RT RT

(10)

where m is an average segment number for the mixture and a0disp is an average dispersion energy per mole of segments. The average segment number for the mixture is determined from

2920 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998

m)

∑i ∑j

(

)

mi + mj

zizj

2

(11)

where zi is the superficial mole fraction of component i (that is, the mole fraction if there was no association), and the average dispersion energy is given by

adisp 0

)

RT

[ ][ ] u

s

∑s ∑t Dst kT

η

t

(12)

τ

where the mixing rule used for the average segment energy is

[] uij

u

)

∑i ∑j zizjmimj kT (v0)ij ∑i ∑j zizjmimj(v0)ij

kT

(13)

molecular orbital calculations for cross dimers. The use of this mixing rule allows these types of mixtures to be fit using only a single adjustable parameter (kij). Derivation of the Mixing Rule. Earlier (Wolbach and Sandler, 1997b), we concluded from the results of our molecular orbital calculations that the equilibrium constant for the cross-association reaction was the geometric average of the equilibrium constants for selfassociation (eq 4). We have also shown that the equilibrium constant for self-association can be related to the SAFT quantity ∆AiBj for any of the association schemes considered by the use of eq 8. To derive a mixing rule, we combined the results of eqs 4 and 8 and assumed that the equilibrium constant for cross association was related to the ∆ parameter for cross association in the following way:

Kij ) xKiKj )

x∆A B ∆A B xC′iC′jRT i i

j j

)

∆AiBj

(20)

xC′iC′jRT

We then approximated the ∆ parameter for cross association as

where

(v0)ij )

[21[(v )

0 1/3 i

+ (v0)1/3 j ]

3

]

( )

(14)

and

uij ) (1 - kij)(uiiujj)1/2

(15)

The binary interaction parameter kij in eq 15 is the only adjustable parameter we used to correlate mixture data. The Helmholtz energy due to association is a mole fraction-weighted average of that for pure components, as described by Chapman et al. (1990)

a

assoc

)

RT

[(

ln XA ∑i zi ∑ A i

) ]

Ai

X

2

i

1 + Mi 2

XAi ) [1 +

zjFXB ∆A B ]-1 ∑j ∑ B j

i j

(17)

[ ( ) ]

[21[(v )

∞ 1/3 i

AiBj )

(18)

+ (v∞j )1/3]

3

]

(

)

[AiBi + AjBj] (22) 2kT

AiBi + AjBj 2

(23)

However, the mixing rule for the well width is more complicated due to the mixture segment volume not being the geometric average of the pure segment

where

v∞ij )

j j

where the first equality is the definition of the ∆ parameter for cross association and the second results from our assumptions. Equating the exponential terms in eq 22 leads to the mixing rule for the well depth of cross association:

AiBj

-1 kT

AiBj ) kT

xκA B κA B x2xv∞i v∞j exp i i

The association strength ∆AiBj is calculated from

∆AiBj ) gij(dij)x2v∞ij κAiBj exp

( )

(16)

j

(21)

where we have assumed that the hard-sphere radial distribution function evaluated at contact is equal to 1 and have neglected the -1 relative to the exponential (AB . kT, see eq 7). Since we have determined the pure-component parameters using low-temperature vapor data, these assumptions should cause negligible errors. This allowed us to express ∆AiBj as

∆AiBj ) κAiBjx2v∞ij exp

where XAi is the mole fraction of molecules of type i not bonded at site A in the mixture and is given by

AiBi kT

∆AiBi = κAiBix2v∞i exp

(19)

volumes (v∞ij ) xv∞i v∞j ), but rather being given by eq 19. Consequently, the mixing rule for the well width is

Further details about the equation of state for mixtures can be found in Huang and Radosz (1991). Binary Mixtures of Self-associating Compounds We initially considered mixtures of two self-associating compounds. To correlate the mixture data, we needed to determine the values of the SAFT crossassociation parameters AiBj/k and κAiBj describing bonding between sites on molecules of two different species. We determined the values of these parameters from the values of the SAFT parameters for self-association by the use of a mixing rule based on the results of our

xκA B κA B xv∞i v∞j i i

κAiBj )

xκA B κA B xv∞i v∞j

j j

v∞ij

i i

)

j j

3 1 ∞ 1/3 [(vi ) + (v∞j )1/3] 2

(

)

(24)

which gives the cross association well width in terms of only the pure-component SAFT parameters. We used these mixing rules to determine the SAFT association parameters for all allowed site-site cross associations. Table 4 displays the allowed site-site bonding combinations for mixtures of compounds described by the association models used in this work.

Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 2921 Table 4. Mixture Association Schemes Used in This Work compound compound 1 2

mixture 1A + 1A

A1

A2

O C

H

1A + 2B

C

H

R

O A1

O

A2

•• O

••

O C

H

R

O

1A + 3B

•• O

•• H

R

O

1A + 4C

A1

C2

∆A1A1 * 0 ∆A2B2 ) ∆A2C2 * 0 ∆A1A2 ) ∆A1B2 ) ∆A1C2 * 0 others ) 0

H

R C2

D2

••

••

O

O

C

H

R

O

H

H

B2

B2

A2

B1

H

A2

•• O

••

••

•• O

R

H

C2

R

3B + 4C

C2

••

•• O

H

C1

O

H

R

D2

••

B1

••

A1

∆A1A1, ∆A2B2 * 0 ∆A1A2 ) ∆A1B2 * 0 others ) 0

H

B2

A2

C

B2

R

A1

O

A1

∆A1A1, ∆A2A2, ∆A1A2 * 0

O

R

2B + 3B

∆ assumptionsa

A2

H B2

∆A1A1 * 0 ∆A1C2 ) ∆A2D2 ) ∆B2C2 ) ∆B2D2 * 0 ∆A1A2 ) ∆A1B2 ) ∆A1C2 ) ∆A1D2 * 0 others ) 0 ∆A1B1 * 0 ∆A2C2 ) ∆B2C2 * 0 ∆A1C2 ) ∆B1A2 ) ∆B1B2 ) 0 others ) 0 ∆A1C1 ) ∆B1C1 * 0 ∆A2C2 ) ∆B2C2 ) ∆A2D2 ∆B2D2 * 0 ∆A1A2 ) ∆A1B2 ) ∆B1A2 ) ∆B1B2 ) ∆C1C2 ) ∆C1D2 * 0 others ) 0

∆ ) zero indicates that association does not occur between sites. The association scheme nomenclature is from Huang and Radosz (1990). a

Fitting Procedure. In the data correlations that follow, we performed either bubble-point pressure or bubble-point temperature calculations, depending on whether the data used had been obtained at constant temperature or pressure. The objective function minimized when fitting the data was, respectively,

z)

|Pcalc - Pexp| + Pexp

∑

∑|ycalc - yexp|

(25)

or

z)

∑|Tcalc - Texp| + 100∑|ycalc - yexp|

(26)

where y represents the superficial vapor-phase mole fractions. The only adjustable parameter was the binary interaction parameter of eq 15, as all other equation of state parameters for the mixture were determined from the values of the pure-component parameters. Mixtures of Acids The first set of cross-associating systems we examined were mixtures of two carboxylic acids. Each of the acids in the mixture was treated using the 1A association model. These compounds should form nearly ideal mixtures, as they are chemically very similar. Table 5 presents the results we obtained for bubble-pressure and bubble-temperature calculations for these mixtures. We see that smaller errors are obtained when using either set of quantum mechanically derived parameters than by using the originally published SAFT parameters. Also, the quantum mechanically derived param-

Figure 1. Propionic acid + formic acid at 760 mmHg. The u0/k for propionic acid was adjusted to match the pure component boiling temperature, and the data recorrelated. For the original SAFT, u0/k for propionic acid is 291.25, when using the H-F results u0/k is 277.25, and when using the B3LYP results u0/k is 300.75. Experimental data from Gmehling and Onken (1977).

eters resulted in smaller values for the binary interaction parameter. Figure 1 displays the results obtained for a mixture of formic and propionic acids at 760 mmHg. The parameters in Table 3 did not reproduce the experimental boiling point of propionic acid for this mixture, so we adjusted the value of u0/k for propionic acid to match the pure-component boiling point. The reoptimized binary interaction parameters are listed in the legend of Figure 1. The improvement obtained when using the quantum mechanically derived parameters is most noticeable for mixtures that include formic acid. This is because the carboxylic acids form nearly ideal mixtures, which requires that the two components have attractive forces that are nearly equal both in magnitude and character. This implies that the acids should have similar self-association equilibrium constants. We can invert eq 8 to compute the value of the equilibrium constant at 298 K from the values of the association parameters of the original SAFT model and from our two sets of quantum mechanically derived parameters. Using the originally published SAFT parameters, the values of Keq (298 K) for formic acid, acetic acid, and propionic acid were 65 902, 9.1, and 0.3, respectively. This indicates that the smallest carboxylic acids behave quite differently on hydrogen bonding as there is a difference of 5 orders of magnitude between the equilibrium constants for formic and propionic acids. In contrast, the equilibrium constants derived from either set of quantum mechanical parameters show a difference of less than a factor of 4. With this much smaller difference in the equilibrium constants, the SAFT parameters derived from our molecular orbital calculations are much better able to capture the near ideality of mixtures of these acids.

2922 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 Table 5. Correlations for Mixtures of Acidsa acid 1

a

acid 2

T (°C)

P (mmHg)

method

kij

∆T (K)

∆y

formic

acetic

100-119

760

SAFT H-F B3LYP

0.042 0.015 0.012

0.58 0.31 0.27

0.021 0.006 0.009

formic

propionic

101-136

760

SAFT H-F B3LYP

0.101 0.050 0.053

2.56 1.33 1.23

0.103 0.022 0.034

acetic

propionic

118-139

760

SAFT H-F B3LYP

0.015 0.018 0.021

0.41 0.45 0.46

0.034 0.019 0.021

acid 1

acid 2

T (°C)

P (mmHg)

method

kij

∆P (mmHg)

∆y

formic

1-pentanoic

100

27-665

SAFT H-F B3LYP

0.116 0.089 0.082

23.3 18.7 15.5

0.054 0.013 0.035

Data from Gmehling and Onken (1977).

Table 6. Bubble Temperature Correlations for Mixtures of Acids + a Second Associating Compounda T (°C)

P (mmHg)

method

kij

∆T (K)

∆y

formic

water

100-108

760

SAFT H-F B3LYP

-0.041 -0.034 -0.102

0.56 0.60 0.30

0.008 0.015 0.009

butyric

water

51-108

100

SAFT H-F B3LYP

-0.050 0.003 -0.132

3.73 1.65 3.36

0.053 0.041 0.055

butyric

water

99-164

760

SAFT H-F B3LYP

-0.043 0.016 -0.112

3.90 1.60 3.83

0.042 0.029 0.046

acetic

methanol

65-115

760

SAFT H-F

-0.070 0.069

1.11 1.81

0.016 0.031

propionic

methanol

65-139

760

SAFT H-F

-0.038 0.088

1.72 2.56

0.025 0.031

acid

a

compound 2


Mixtures of an Acid with Separate Water and Methanol We next considered mixtures of acids with either water or methanol. When using the quantum mechanically derived SAFT parameters, water was treated with the 4C association model, while methanol was treated with the 3B association model. These were the association models for these compounds found earlier (Wolbach and Sandler, 1997c,d) to yield the best correlation of pure-component VLE data when using the results of our quantum mechanical calculations. We used these same association models here when considering mixtures. Thus, water was treated with the 3B association model and methanol was treated with the 2B association model. A molecule described by the 4C association model possesses two hydrogen-bond donor and two hydrogen-bond acceptor sites. A molecule described by the 3B model possesses one donor and two acceptor (or two donor and one acceptor) sites, while a molecule described by the 2B association model possesses one donor and one acceptor site. The acids were treated with the 1A association model regardless of the parameter set used. The use of different association models for water and methanol complicates the comparison between the correlations obtained using the original SAFT parameters and those obtained from our molecular orbital calculations. However, as each set of parameters and association models performed equally well for the pure components, we expect this complication to be minimal. Table 6 presents the results of bubble-temperature calculations for these mixtures, and Table 7 displays

Figure 2. Acetic acid + water at 40 °C. Experimental data from Gmehling and Onken (1977).

the results of bubble-pressure calculations. When methanol was one of the components in a mixture, we only correlated the data using SAFT parameters derived from the results of our H-F calculations. Once again, the use of the quantum mechanically derived SAFT parameters results in the correlation of the experimental data of a quality equal to or better than that obtained using the originally published SAFT

Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 2923 Table 7. Bubble Pressure Correlations for Mixtures of Acids + a Second Associating Compounda acid

a

compound 2

T (°C)

P (mmHg)

method

kij

∆P (mmHg)

∆y

acetic

water

40

35-56

SAFT H-F B3LYP

-0.070 -0.016 -0.138

1.29 1.16 1.18

0.022 0.022 0.024

acetic

water

60

106-146

SAFT H-F B3LYP

-0.074 -0.030 -0.139

3.09 2.09 2.35

0.018 0.014 0.016

propionic

water

100

191-773

SAFT H-F B3LYP

-0.016 0.000 -0.105

30.73 17.53 26.49

0.036 0.020 0.027

propionic

water

60

29-151

SAFT H-F B3LYP

-0.025 0.007 -0.116

6.29 3.45 4.90

0.037 0.027 0.032

acetic

methanol

35

26-211

SAFT H-F

-0.021 0.122

2.24 2.38

0.021 0.006


Table 8. Results of Correlations for Methanol + Water Mixturesa

a

∆P (mmHg)

∆y

system

T (°C)

P (mmHg)

method

kij

MeOH + H2O

25

26-127

SAFT H-F

-0.138 0.027

MeOH + H2O

100

780-2529

SAFT H-F

-0.112 0.017

26.1 62.3

0.009 0.011

MeOH + H2O

150

3790-10 302

SAFT H-F

-0.089 0.014

59.6 114.0

0.014 0.017

MeOH + H2O

200

12 204-29 581

SAFT H-F

-0.075 0.005

195.5 262.5

0.013 0.012

0.63 0.93

0.004 0.009

system

T (°C)

P (mmHg)

method

kij

∆T (K)

∆y

MeOH + H2O

64.7-100

760

SAFT H-F

-0.100 0.036

0.42 0.43

0.014 0.009


Table 9. Parameters for the Diluents in These Cross-associating Mixturesa system

method

∆H (kcal/mol)


∆Cp0 (cal/(mol‚K))

acetone + MeOH diethyl ether + MeOH acetone + H2O

H-F H-F H-F B3LYP

-4.658 -4.004 -4.623 -4.033

-24.709 -25.108 -23.365 -24.300

1.27 1.44 1.32 1.19

compound

v∞ (cm3/mol)

m

u0/k (K)

acetone diethyl ether water (H-F: 3B)

7.765 10.220 8.80

4.504 4.430 1.278

210.92 192.50 385.12

system

a

method

κAB × 102

AB/k (K)

2286 AB/k

(K)

3.091 κAB

× 102

acetone + MeOH diethyl ether + MeOH

H-F H-F

1824 1662

4.661 2.003

acetone + H2O

H-F (4C) H-F (3B) B3LYP (4C)

1947 1913 1631

2.989 6.317 3.360


parameters, though here there is no consistent trend in the magnitude of the binary interaction parameters. Figure 2 displays the results for a mixture of acetic acid and water at 40 °C. While all three parameter sets yield small average deviations in the calculated bubble pressure and vaporphase composition, none of the resulting fits are able to capture the shape of the vapor-liquid-phase envelope. This may be due to the limitations of the association model used to represent acids, as all pure-acid dimers and cross-associated dimers are forced to be cyclic. We also correlated the data using the 2B association model for acetic acid, which would allow the formation of linear polymers, but found this to yield little improvement.

Mixtures of Methanol and Water The final mixtures of self-associating compounds we considered were of methanol and water. We used the SAFT parameters derived from the results of our H-F calculations and treated water with the 4C association model and methanol with the 3B association model. When using the originally published SAFT parameters, water was treated with the 3B association model and methanol with the 2B association model. Table 9 presents the results for these mixtures at a variety of conditions. Although both parameter sets fit the experimental data well, use of the original SAFT parameters generally yields smaller errors, though with larger values of the binary interaction parameter. The


constants, eq 27c is the balance on the true mole fractions, and eqs 27d and 27e are the component mass balances. The quantities Nbefore and Nafter are the number of molecules before and after association, respectively. The unknowns in this system of equations are the four true mole fractions and the ratio of Nbefore/ Nafter. Solving this system of equations, the mole fraction of the original acetone molecules that are unbonded is Figure 3. Bonding between methanol and acetone. The sites assigned to each molecule are labeled. The bond shown is the only allowed intermolecular bond.

results in Table 8 also indicate that the kij for these mixtures may be approximated as a linear function of temperature. Cross-associating Mixtures of an Associating and a Nonassociating Compound The final class of mixtures we considered were of a self-associating compound and a non-self-associating compound that can cross associate. Examples of this type of system are methanol or water + acetone and methanol + diethyl ether. This type of mixture is difficult to treat using the original SAFT model, as one cannot derive the cross-association parameters from a mixing rule, so that at least one additional adjustable parameter would have to be introduced. We are able to use the results of our molecular orbital calculations to determine the SAFT parameters for cross association a priori. The procedure is straightforward, although involved. We used the values of the quantum mechanically derived equilibrium constants to determine the mixture vapor-phase composition, assuming that only monomers and dimers are present. We then used this information to determine the mole fraction of unbonded molecules of the non-self-associating fluid at a number of temperatures, and then a linear leastsquares fit is used to determine the SAFT parameters for cross association. As an example, consider a mixture of acetone and methanol. Methanol is treated with the 3B association model, and only one of its association sites binds to the acetone. Acetone is assumed to have a lone pair as its one association site. This is depicted in Figure 3. To determine the true vapor-phase composition at a given temperature, the following system of equations was solved (in which methanol is represented by M and acetone by Ac):

KM-MP )

KM-AcP )

x2,M (x1,M)2

x2,M-Ac x1,Mx1,Ac

x1,M + x1,Ac + x2,M + x2,M-Ac ) 1

(27a)

(27b) (27c)

(x1,M + 2x2,M + x2,M-Ac)Nafter ) zMNbefore (27d) (x1,A + x2,M-Ac)Nafter ) zAcNbefore

(27e)

In these equations, xi represents the true mole fraction of a species and zi represents its superficial mole fraction which is known from the experimental data. Equations 27a and 27b are the definitions of the two equilibrium

XA2 )

x1,AcNafter zAcNbefore

(28)

The quantities of interest are the two SAFT parameters for cross association, A2C1/k and κA2C1. Both of these parameters are contained in ∆A2C1, and we solved a second system of equations:

XC1 )

XA1 )

1 1 + FX zM∆A1C1

(29a)

XB1 )

1 1 + FX zM∆A1C1

(29b)

C1

C1

1 1 + FXA1zM∆A1C1 + FXB1zM∆A1C1 + FXA2zAc∆A2C1 (29c) XA2 )

1 1 + FXC1zM∆A2C1

(29d)

to determine the value of the ∆A2C1. The unknowns in these equations are XA1, XB1, XC1, and ∆A2C1. After determining values of ∆A2C1 at a number of temperatures, we used eq 24 and a linear least-squares fit of ln(∆A2C1) versus 1/T to determine the values of the cross-association parameters. The equations necessary to determine the SAFT cross-association parameters for this type of mixture using other association models to describe the self-associating compound are given in the Appendix. We used the originally published SAFT parameters of Huang and Radosz (1990) for the nonassociating compounds; the pure-component parameters are given in Table 9. Also in this table are the results of our molecular orbital calculations for the cross dimers, along with the derived SAFT parameters for cross association for each system. Table 10 presents the average deviations from experimental data using the SAFT model with both our quantum mechanically derived parameters and the originally reported SAFT parameters. When correlating data using the original SAFT model, we assumed that the cross-association parameters were equal to zero in order to have only a single adjustable parameter in the fitting procedure and to test the importance of explicitly accounting for the cross association in these mixtures. We first considered a pair of systems in which methanol was the self-associating component and acetone or diethyl ether was the diluent. For these compounds, a single lone pair was assumed to be the only bonding site, and the association shown in Figure 3 was assumed to be the only allowed cross association. We used only parameters derived from the results of our H-F calculations when correlating the data, and

Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 2925 Table 10. Results for Mixtures of Associating Compounds and Diluents Where Cross Association Is Possiblea system MeOH + acetone MeOH + (C2H5)2O water + acetone

a

kij

∆T (K)

∆y

SAFT (2B) H-F (3B)

-0.050 -0.025

1.09 0.50

0.026 0.006

760

SAFT (2B) H-F (3B)

-0.026 -0.020

1.57 0.75

0.026 0.028

760

SAFT (3B) H-F (4C) H-F (3B) B3LYP (4C)

-0.145 -0.165 -0.050 -0.071

4.87 4.71 1.11 2.14

0.084 0.090 0.013 0.044

T (°C)

P (mmHg)

55-65

760

307-337 56-100

method (model)


Figure 4. Diethyl ether + methanol at 760 mmHg. Experimental data from Gmehling and Onken (1977).

Figure 5. Methanol + acetone at 760 mmHg. Experimental data from Gmehling and Onken (1977).

the resulting correlations of the experimental data are displayed in Figures 4 and 5. For the methanol systems we were able to explicitly account for the cross association and correlate azeotropic mixture data using a single binary interaction parameter with an accuracy greater than that of the SAFT

Figure 6. Water + acetone at 760 mmHg. Experimental data from Gmehling and Onken (1977).

model with original parameters. While either model is able to reproduce the proper azeotropic compositions for these two systems, the original SAFT underestimates the azeotropic temperature to a greater extent, which is equivalent to underestimating the extent of association in the liquid phase at the azeotropic composition. By explicitly accounting for the cross association, we are much better able to represent the azeotropic temperature. Figure 4 does show departure from the fitted data, even when explicitly accounting for the cross association. However, since the shape of the phase envelope is reproduced accurately, this may indicate that the SAFT model is correctly accounting for the associated species formed in the mixture. Using our quantum mechanically derived parameters, the results consistently underestimate the boiling temperature at a given liquid composition. This is equivalent to underestimating the extent of association in the liquid phase, which suggests that the quantum mechanical calculations may underestimate the strength of the association between methanol and diethyl ether. We next examined the water + acetone system, and the results we obtained are shown in Figure 6. While the correlation using the results of our B3LYP calculations is very good, we were not able to achieve the same degree of accuracy when using the results of our H-F calculations. The SAFT model with original parameters also performs poorly for this system. One possible explanation for the poor performance of


Figure A1. Bonding between a 4C molecule (water) and a generic carbonyl-containing compound. The sites assigned to each molecule are labeled. The bond shown is the only allowed type of intermolecular bond.

Figure 7. Water + acetone at 760 mmHg using the results of our H-F calculations. The legend refers to the association model used for water. Experimental data from Gmehling and Onken (1977).

the H-F parameters is the use of the 4C association model for water. In Figure 6, the experimental data shows a sharp decrease in bubble temperature with the addition of small amounts of acetone. This suggests a lowering of the extent of association in the liquid phase of the mixture, as the vapor phase of this mixture exhibits only small degrees of association at 760 mmHg. The 4C association model for water allows for large degrees of self-association and cross-association to occur simultaneously and may result in an overestimation of the extent of association in the liquid phase of this mixture. To correct this problem, we recorrelated the data using the 3B association model for water, assigning two lone pairs and one proton to be the bonding sites for water, as was the case for methanol. This reduction in the number of proton donor sites reduces the extent of association in the liquid phase of this mixture. The optimized SAFT parameters for water using this model are given in Table 9, the average errors in Table 10, and the correlation displayed in Figure 7. This association model does produce much better results when using the H-F derived parameters. Conclusions This work demonstrates the usefulness of quantum mechanically derived parameters of association for describing the phase behavior of mixtures containing hydrogen-bonding fluids in the SAFT equation of state. We have developed a mixing rule based on the results of our molecular orbital calculations that allows the SAFT parameters for cross association to be determined from the values of the self-association parameters for the species in the mixture. Using this mixing rule, we were able to correlate VLE data for mixtures of selfassociating compounds with an accuracy equal to or greater than that of the original SAFT model. In addition, we used the results of our molecular orbital calculations to determine a priori the SAFT parameters for cross association for mixtures containing a self-associating compound and a diluent that could

form cross dimers. We were able to correlate VLE data for mixtures of this type with a single adjustable parameter, and the correlations obtained were of accuracy greater than that of the original SAFT model also using a single adjustable parameter. We have also found that it is not necessary to use the more rigorous B3LYP calculations to obtain accurate correlations of the mixture data. Appendix In this appendix we will first derive the equations for a mixture of a species described by the 4C association model (two lone pairs + two protons) and a species assumed to have a lone pair as its one association site. We then repeat the derivation for a mixture of a species described by an alternate 3 site association model (2 protons + one lone pair) and a species assumed to have a lone pair as its one association site. Application to the 4C Association Model. The situation to be described is shown in Figure A1. To determine the true vapor-phase composition at a given temperature, the following system of equations was solved (in which water is represented by W and the carbonyl compound by Cb)

KW-WP )

KW-CbP )

x2,W (x1,W)2

x2,W-Cb x1,Wx1,Cb

x1,W + x1,Cb + x2,W + x2,W-Cb ) 1 (x1,W + 2x2,W + x2,W-Cb)Nafter ) zWNbefore (x1,Cb + x2,W-Cb)Nafter ) zCbNbefore

(A.1a)

(A.1b) (A.1c) (A.1d) (A.1e)

The origin of, and notation in, these equations are as in eq 27. The mole fraction of the original carbonyl-containing molecules that were unbonded is

XA2 )

x1,CbNafter zCbNbefore

(A.2)

The quantities in which we are interested are contained in the two SAFT parameters for cross association, A2C1/k and κA2C1. These parameters are assumed to be equal to the parameters A2D1/k and κA2D1. We then solved a system of equations to determine the value of the ∆A2C1

Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 2927

XC1 ) XA2 )

1 A1

A1C1

B1

A2C1

1 + FX zW∆

+ FXA2zCb∆A2C1

1 1 + FX zW∆

+ FXC1zW∆A2C1

(A.4c) (A.4d)

The unknowns in these equations are XA1, XB1, XC1, and ∆A2C1. Figure A2. Bonding between a 3B* molecule (water) and a generic carbonyl-containing compound. The sites assigned to each molecule are labeled. The bond shown is the only allowed type of intermolecular bond.

XA1 ) XB1 ) XC1 )

XD1 )

1 C1

A1C1

C1

A1C1

1 + FX zW∆

D1

A1C1

+ FX zW∆

1 1 + FX zW∆

A1

A1C1

1 + FX zW∆

+ FXD1zW∆A1C1

(A.3a)

(A.3b)

1 + FX zW∆A1C1 + FXA2zCb∆A2C1 (A.3c) B1

1 1 + FXA1zW∆A1C1 + FXB1zW∆A1C1 + FXA2zCb∆A2C1 (A.3d) XA2 )

1 C1

A2C1

1 + FX zW∆

+ FXD1zW∆A2C1

(A.3e)

The unknowns in these equations are XA1, XB1, XC1, XD1, and ∆A2C1. After determining values of ∆A2C1 at a number of temperatures, we used the approximation in eq 24 and performed a linear least-squares fit of ln(∆A2C1) versus 1/T to determine the values of the cross-association parameters. Application to an Alternate 3B Model. In the text we considered a three-site model with two lone pairs and a single hydrogen as the bonding sites. Here, we consider a three-site model with two hydrogens and a single lone pair as the bonding sites. The model is labeled 3B* and could be used to describe water. The specification of bonding sites has no effect on the purecomponent equation, since no explicit sites have been assigned in the SAFT association designations. The situation to be described is shown in Figure A2. Again, we label the water as W and the carbonyl compound as Cb. We may use eqs A.1a-e to determine the true vapor-phase compositions, and eq A.2 to calculate the fraction of unbonded carbonyl compound. We are interested in the two SAFT parameters for cross association, A2C1/k and κA2C1. Both of these parameters are contained in ∆A2C1, and we solve the system of equations below to determine the value of the ∆A2C1

XA1 ) XB1 )

1 B1

A1C1

1 + FX zW∆

+ FXC1zW∆A1C1

(A.4a)

1 (A.4b) 1 + FXA1zW∆A1C1 + FXA2zCb∆A2C1

Literature Cited Becke, A. D. Density Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98 (7), 5648. Chapman, W. G.; Jackson, G.; Gubbins, K. E. Phase Equilibria of Associating Fluids: Chain Molecules with Multiple Bonding Sites. Mol. Phys. 1988, 65, 1057. Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. SAFT: Equation-of-State Solution Model for Associating Fluids. Fluid Phase Equilib. 1989, 52, 31. Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. New Reference Equation of State for Associating Liquids. Ind. Eng. Chem. Res. 1990, 29, 1709. Economou, I. G.; Donohue, M. D. Chemical, Quasi-Chemical and Perturbation Theories for Associating Fluids. AIChE J. 1991, 37, 1875. Economou, I. G.; Donohue, M. D. Equation of State with Multiple Associating Sites for Water and Water-Hydrocarbon Mixtures. Ind. Eng. Chem. Res. 1992, 31, 2388. Gaussian 92/DFT, Revision G.2; Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Wong, M. W.; Foresman, J. B.; Robb, M. A.; Head-Gordon, M.; Replogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzalez, C.; Martin, R. L.; Fox, D. J.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J. A. Gaussian, Inc.: Pittsburgh, PA, 1993. Gaussian 94, Revision B.2; Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanayakkara, A.; Challacombe, M.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; Head-Gordon, M.; Gonzalez, C.; Pople, J. A. Gaussian, Inc.: Pittsburgh, PA, 1995. Gmehling, J.; Onken, U. Vapor-Liquid Equilibrium Data Collection; DECHEMA: Frankfurt, FRG, 1977. Huang, S. H.; Radosz, M. Equation of State for Small, Large, Polydisperse, and Associating Molecules. Ind. Eng. Chem. Res. 1990, 29, 2284. Huang, S. H.; Radosz, M. Equation of State for Small, Large, Polydisperse, and Associating Molecules: Extension to Fluid Mixtures. Ind. Eng. Chem. Res. 1991, 30, 1994. Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785. McQuarrie, D. A. Statistical Mechanics; Harper-Collins: New York, 1976. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. I. Statistical Thermodynamics. J. Stat. Phys. 1984a, 35, 19. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. II. Thermodynamic Perturbation Theory and Integral Equations. J. Stat. Phys. 1984b, 35, 35. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. III. Multiple Attraction Sites. J. Stat. Phys. 1986a, 42, 459. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. IV. Equilibrium Polymerization. J. Stat. Phys. 1986b, 42, 477.

2928 Ind. Eng. Chem. Res., Vol. 37, No. 8, 1998 Wolbach, J. P.; Sandler, S. I. Thermodynamics of Hydrogen Bonding from Molecular Orbital Theory. 1. Water. AIChE. J. 1997a, 43 (6), 1589. Wolbach, J. P.; Sandler, S. I. Thermodynamics of Hydrogen Bonding from Molecular Orbital Theory. 2. Organics. AIChE. J. 1997b, 43 (6), 1597. Wolbach, J. P.; Sandler, S. I. Using Molecular Orbital Calculations to Describe the Phase Behavior of Hydrogen-Bonding Fluids. Ind. Eng. Chem. Res. 1997c, 36, 4041.

Wolbach, J. P.; Sandler, S. I. The Use of Molecular Orbital Calculations to Describe the Phase Behavior of HydrogenBonding Mixtures. Int. J. Thermophys. 1997d, 18 (4), 1001.

Received for review November 5, 1997 Revised manuscript received April 14, 1998 Accepted April 14, 1998 IE970781L

Using Molecular Orbital Calculations To Describe the Phase Behavior

Recommend Documents