A molecular equation of state for alcohols which includes steric hindrance in hydrogen bonding
11 A molecular equation of state for alcohols which includes steric hindrance in hydrogen bonding
Bennett D. Marshall
ExxonMobil Research and Engineering, 22777 Springwoods Village Parkway, Spring TX 77389 USA
Abstract
In this paper we develop the first equation of state for alcohol containing mixtures which includes the effect of steric hindrance between the two electron lone pair hydrogen bond acceptor sites on the alcohol’s hydroxyl oxygen. The theory is derived for multi-component mixtures within Wertheim’s multi-density statistical mechanics in a second order perturbation theory. The accuracy of the new approach is demonstrated by application to pure methanol and ethanol and binary ethanol / water mixtures. It is demonstrated that the new approach gives a substantial improvement in the prediction of hydrogen bonding structure of both pure alcohol and alcohol/water mixtures as compared to conventional approaches which do not include steric effects between the alcohol association sites. Finally, it is demonstrated that the inclusion of steric effects allows for more accurate binary phase equilibria and heats of mixing prediction with water. [email protected]
I: Introduction
Much of current thermodynamics research is focused on molecular simulations: new algorithms, force fields, applications, etc.… Indeed, molecular simulations allow for the “exact” solution of the thermodynamics / fluid structure for a given force field. However, even with the advances in molecular simulation, it is the general mixture equation of state (EoS) which forms the backbone of commercial thermodynamics packages used to predict fluid phase equilibria. EoS can be used to quickly estimate the phase equilibria of multi-component mixtures, which makes them amiable for implementation into process simulators. The predictive ability of an EoS often relies on a proper accounting of the physics of the relevant intermolecular interactions. Perturbation theory provides a relatively simple and accurate framework in which to develop EoS. Weeks, Chandler, Anderson and Barker, Hendersen where the first widely applied perturbation theories for spherically symmetric fluids. Extension of Perturbation theory to hydrogen bonding (associating) fluids is challenging due to the strength, anisotropy and limited valence of the hydrogen bond. Andersen , Dahl and Andersen and Chandler and Pratt where the first to develop statistical mechanical formalisms which could be used to develop perturbation theories for associating fluids. However, these approaches have not been widely applied because the final theoretical results are in the form of cluster expansions. In the 1980’s Wertheim developed a new multi-density form of statistical mechanics, where each hydrogen bonding state of a molecule is assigned a separate density. Taking this approach, Wertheim was able to develop a rigorous statistical mechanics formalism for associating fluids. The beauty of Wertheim’s approach is that general solutions to the theory can be obtained in perturbation theory (TPT). Chapman obtained a general solution to TPT in first order (TPT1) for a multi-component mixture where each molecule can have an arbitrary number and functionality of association sites. The simplicity and generality of this TPT1 solution has allowed for its wide industrial and academic application in the form of the statistical associating fluid theory (SAFT) class of EoS as well as the Cubic-Plus-Association EoS. While widely applied, TPT1 has many limitations. First order perturbation theory assumes that each association bond in a cluster is independent of the other association bonds. This allows for a simple general solution, but neglects effects such as steric hindrance , ring formation , double bonding of molecules and association sites which can receive multiple association bonds . To include these effects, one must go to higher order in perturbation. TPT2 allows for the interaction of two association bonds in a cluster, TPT3 allows for the simultaneous interaction of 3 association bonds etc.… Alcohols are a common class of associating molecules. Each alcohol has two hydrogen bond acceptor sites (oxygen lone pairs) and one hydrogen bond donor (hydroxyl hydrogen). This defines a 3C association scheme. There has been much debate in the literature on whether alcohols should be modeled with two (one donor and one acceptor in a 2B scheme) or three association sites.
Recently, Fouad et al. demonstrated that within the polar PC-SAFT equation of state, the 2B scheme gave better agreement with hydrogen bonding distributions in pure alcohols. However, it was simultaneously shown that the 3C model gave better agreement for hydrogen bonding distributions for alcohol-water mixtures, as well as well as improved phase behavior predictions with water. The results of Fouad et al. point to a deficiency in the application of TPT1 to alcohols. As TPT1 does not include steric effects, it will not account for the fact that bonding at one of the oxygen acceptor sites may partially block accessibility to the other. This would lead to an overprediction in the fraction of ethanol molecules which are bonded at both acceptor sites. Hence, a 2B model, which rejects second acceptor site completely, gave better results than the 3C model. However, when inserting an alcohol molecule in an aqueous phase, it may become bonded at all three sites. What is needed is a theoretical approach which accounts for the steric hindrance between the two oxygen sites in a 3C model. This will be the subject of this paper. In this paper we develop the first general multi-component solution of TPT2 which allows for the incorporation of steric effects. Each species can have an arbitrary number and functionality of association sites; however, for tractability we only allow one pair of sites per molecule to become sterically coupled. We then specialize this solution to mixtures which contain alcohols modelled with a 3C association scheme. We then incorporate the new perturbation theory into the perturbed chain statistical associating fluid theory (PC-SAFT) EoS and use it to study the hydrogen bonding and phase behavior of pure alcohol and alcohol-water binary mixtures. In this paper we focus on the alcohols methanol and ethanol. We demonstrate that the inclusion of steric hindrance through TPT2 gives a significantly improved EoS as compared to TPT1. II: Thermodynamic perturbation theory
In this section we extend thermodynamic perturbation theory (TPT) to account for steric hindrance between association sites in a multi-component fluid. We consider a mixture of N molecules consisting of n separate species of number density ρ ( k ) . Each species contains a set of Г (k) = {A, B, C,…,G} association sites, where the capitals letters represent distinct association sites. While each species can have any number and type of association sites, we restrict the theory such that only a single pair of association sites per species can sterically hinder each other. For molecules with a 3C association scheme, we call this approach the 3C-SH association model. Figure 1 gives a model representation of a 3C-SH alcohol model. There is one donor hydrogen labeled H and two acceptor oxygens lone pairs which we label O and O . However, only the two oxygen sites exhibit steric effects. Figure 1:
Example representation of ethanol using a 3-C association scheme. The two lone pairs of electrons are treated as acceptor sites (blue) and the hydrogen is a donor site (red) The potential of interaction between a molecule 1 of species k and a molecule 2 of species j is given by ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , , ,12
12 12 k j k j k j k jhs ABA B r = + (1) The distance between the centers of the molecules is r and the notation (1) represents the position and orientation of molecule 1. The term φ hs is the pair potential of the spherically symmetric hard sphere reference fluid and ( ) , k jAB is the potential of interaction between site A on species k and site B on species j . The theory is developed in Wertheim’s multi-density formalism. In this approach each bonding state of a molecule is treated as a distinct species and assigned a density 𝜌 𝛼(𝑘) , where α is the set of bonded sites. Hence, 𝜌 𝑜(𝑘) is the monomer density of species k . To aid in the topological reduction from fugacity to density graphs, Wertheim defines a set of density parameters ( ) ( ) k k = (2) where 𝜎 𝑜(𝑘) = 𝜌 𝑜(𝑘) and 𝜎 Г(𝑘) = 𝜌 (𝑘) . The total Helmholtz free energy is given by ( ) ( )( ) ( ) ( ) ( ) ln k ok k khs okkB A A cQVk T V − = + + − (3) where A hs is the free energy of the hard sphere reference fluid, V is the system volume and T is the absolute temperature. Equation (3) is mathematically exact. The fundamental graph sum Δ c (o) is an infinite series of integrals which encodes all association interactions between molecules. In thermodynamic perturbation theory (TPT) all contributions to Δ c (o) which contain more than a single associated cluster are neglected. This allows for the summation of Δ c (o) in terms of reference system correlation functions. TPT assumes that the structure of the fluid is unchanged due to association. Within perturbation theory, Δ c (o) is then ordered in terms of irreducible clusters which contain 1, 2, 3 etc… association bonds ( ) o jj c c = (4) where Δ c j contains irreducible contributions from cluster integrals with k association bonds. In first order perturbation theory (TPT1) Eq. (4) is truncated at j = 1, second order TPT2 at j = 2 etc… At the TPT1 level, Eq. (4) only includes contributions with a single irreducible bond. All reducible clusters can be created from the TPT1 result, however there is no steric hindrance between association sites. For a multi component mixture, where each species can have an arbitrary number and functionality of association sites, the TPT1 contribution is given by ( ) ( ) ( )( ) ( ) ( ) ( ) ,1 1 1 k jk j n n k j k jABA Bk j A B cV − −= = = (5) Where ( ) ( ) ( ) ( ) ( ) ( ) , , , 122 k j k j k jAB AB hs f g r d = (6) The association Mayer function is given by ( ) ( ) ( ) ( ) ,, k jk j ABAB b f k T = − − (7) and ( ) ( ) , 12 k jhs g r is the mixture pair correlation function of the hard sphere reference system. Including the TPT2 term includes irreducible graphs which contain clusters with two association bonds. This is the minimum level of theory to incorporate steric effects. Extending the second order approach of Wertheim to the case of a multi-component fluid with and arbitrary number of association sites ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( )( )( ) , ,2 1 1 1 k j ik j i i n n n k j i k i jACDBA B CDk j i A B C D cV − − −= = = = (8) with ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , , , , , ,12 234, , 12 23 13, ,12 23 k i j k i i j k i i jACDB AC DB hs hsk i jhsk i i jhs hs f f g r g rg r r r d dg r g r = − (9) Where ( ) ( ) , , 12 23 13 , , k i jhs g r r r is the triplet correlation function of the reference fluid. Equation (9) completes the definition of the graph sum Δc (o) . The last term to consider in Eq. (3) is Q ( k ) ( ) ( ) ( ) ( ) ( ) ( ) kk k k k k Q c − = − + (10) The functions 𝑐 𝛾(𝑘) are generated from the graph sum Δc (o) according to the relation ( ) ( ) ( ) ; ok k cc V − = (11) Now we assume that each molecule can have at most one pair of second order sites which sterically hinder each other. We label this set {C,D}. Evaluating Eq. (11) subject to Eqns. (4), (5) and (8) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( )( ) , , , , ,1 1 1 j j ij j i i n n nk j k j j i k i j j i kA AB ACDB BCDAB B CDj j iB B C D c − − −= = = = + + (12) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) , , , ,1 1 i ji j n n i j i k j i k jACDB ADCBA Bi j A BkEF for EF CDc otherwise − −= = + == (13) ( ) ( ) k c for n = (14) A key quantity in TPT1 is the fraction molecules not bonded at site A: 𝑋 𝐴(𝑘) = 𝜎 Г (𝑘) −𝐴(𝑘) /𝜌 . In this second order theory we will also require the fraction of molecules not bonded at both sites C and D : 𝑋 𝐶𝐷(𝑘) = 𝜎 Г (𝑘) −𝐶𝐷(𝑘) /𝜌 . With the current assumption of only a single pair of second order sites the theory will have similar structure to that of Marshall and Chapman (MC) . Generalizing the MC solution to a multi-component mixture we obtain ( ) ( )( ) ( ) ( ) kSkS k kL CD for S C or DcX c X otherwise += + (15) In Eq. (15) when S = C, L = D and when
S = D, L = C . The fractions 𝑋 𝐶𝐷(𝑘) are given by ( ) ( ) ( ) ( ) ( ) ( )
11 1 kCD k k kCD C D
X c c c = + + + (16) The monomer fraction is found to be ( ) ( )( ) ( ) ( ) ( ) k kk kCDo Ak k AC D XX XX X = (17) From these results Eq. (10) can be evaluated as ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) / 1 2 k k kk k k C DA kA CD X XQ X X = − + − (18) Combining these results, we simplify the free energy in Eq. (3) to ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) k kkn nk k khs CDAA k kk kAB C D A A XXXVk T X X = = − = − + + (19) 0 To maintain consistency with the PC-SAFT EoS, equation (6) is evaluated as ( ) ( ) ( ) ( , ), , ,3 exp 1 k jk j k j k jABAB kj AB hsb gk T = − (20) where σ kj is the cross-species diameter, ( ) , k jAB is the bond volume, ( , ) k jAB is the association energy and ( ) , k jhs g is the contact value of the hard sphere pair correlation function evaluated with the Carnahan and Starling EoS. The mixture quantities are evaluated with the following combining rules ( ) ( ) ( ) ( , ) ( , ), , ,3 3 3 ( , ) ; 2 k k j jk j k k j j k j AB ABkj AB kk AB jj AB AB += = (21) Turning our attention to the second order contribution in Eq. (9), we first approximate the triplet correlation function with the following simple superposition ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ,, , , , 1312 23 13 12 23 , , exp k jk i j k i i j hshs hs hs b rg r r r g r g r k T = − (22) The exponential term serves to provide the steric effects between the sites C and D on species i . Combining (9) and (22) we obtain ( ) ( ) ( ) ( ) ( ) , , , , , , k i j k i j k i i jACDB AC DB ACDB = − (23) Where ( ) , , k i j ACDB is the overlap integral defined as the fraction of associated states (where both C and D on i are bonded to site A on k and B on j respectively) which do not lead to overlap. When the sites CD do not sterically hinder each other ( ) , , k i j ACDB = 1. If there is complete steric hindrance, then the association at one site completely blocks the second site and ( ) , , k i j ACDB = 0. Mathematically the overlap integral is evaluated by 1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , , ,, , 13, ,
12 23 exp (2) (3)12 23 (2) (3) k i j k jk i i j hsAC DB bACDB k i i jAC DB rf f d dk Tf f d d − = (24) In practice we shall assume that the overlap integral of component i is independent of the species k and j ( ) ( ) , , k i j i ACDB = (25) The theory developed in this section is general for an arbitrary number of components, with an arbitrary number and functionality of association sites. The one restriction is that each molecule can have at most one pair of association sites which sterically hinder each other. In section III the theory is specialized to alcohols using the 3C-SH association model. 2
III: Application to alcohols
In this section we specialize the theory developed in section II to the case of mixtures which contain alcohols with the 3C-SH association model outlined in Fig. 1. Site H (red) is the hydrogen bond donor site, O is the first lone pair oxygen acceptor site (blue) and finally O is the second lone pair oxygen acceptor site (green). Sites O and O are assumed equivalent and are sterically coupled. Hence, association at O can block O and vice versa. There are no steric effects between either oxygen site and the donor hydrogen H . Specializing Eqns. (12)-(13) to this case and enforcing Eq. (25) ( ) ( ) ( ) ( ) ( ) ( ) ,1 j nk k j j k jO O B O Bj B c c X = = = (26) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , , ,1 1 1 j j n n nk j j k j i j j i k i i j iH B HB B O O HO O Bj j iB B c X X X = = = = + − (27) ( ) ( ) ( ) ( ) ( ) ( ) ( ) j nk j j j k kO O B BOj B c X = = − (28) ( ) ( ) k kO H O H c c = = (29) Specializing Eqns. (15)-(16) to the 3C-SH model ( ) ( ) ( ) ( ) ( ) ( ) ( ) ,1 j nk k k j j k jO O O O B O Bj B X X X X = = = + (30) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , , ,1 1 1
11 2 1 j j kH n n nj j k j i j j i k i i j iB HB A O O HO O Aj j iB A
X X X X = = = = + + − (31) 3 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
11 1 j j kO O n nj j k j k j j k jB O B B O Bj jB B
X X X = = = − + + (32)
Combining Eqns. (30) and (32) we obtain ( ) ( )( ) ( ) ( ) ( )( ) k kO Ok kO O k kO O O O
X XX X X − − + = (33) Solving Eq. (33) ( ) ( )( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( )
22 1 1 4 1 1 1 k kO Ok k k k k kO O O
XX X X = = − + + − − − + (34) Combining Eqns. (30) and (34) ( ) ( ) ( )( ) ( ) ( ) ( ) ,1 j kk kO O n j j k jB O Bj B X X X = = = + (35) Hence the bonding fractions of molecules with the 3C-SH association model are obtained by the simultaneous solution of Eqns. (31) and (35), from which ( ) kO O X can be calculated from Eq. (34). For hydrogen bonding species in a multi-component mixture which do not exhibit steric hindrance, the fraction of molecules which are not bonded at the acceptor site A is obtained from Eq. (35) by setting O = A and δ ( k ) = 1. Similarly, for hydrogen bonding molecules in a multi-component mixture which do not exhibit steric hindrance, the fraction of molecules not bonded at donor site D is obtained from Eq. (31) with the transformation H = D . 4 Figure 2:
Plot of δ versus the fraction of unbonded oxygen sites (either O or O ) for various values of Ψ Analyzing Eq. (34) we see that δ = 1 in the limit Ψ = 1, signifying a reduction to TPT1 in the absence of steric hindrance. Figure 2 plots δ versus the fraction of unbonded oxygen sites (either O or O ) for various values of Ψ. The Ψ → 0 + curve is particularly interesting because δ vanishes identically for 𝑋 𝑂 = 0.5. This demonstrates that the theory successfully reproduces complete blocking of the unbonded receptor site. When steric effects are included in models with a single acceptor and donor association sites, one must go to infinite order in perturbation to reproduce complete blockage. For the current 3C-SH model with steric hindrance between O and O , we obtain convergence at the TPT2 level. IV: Incorporation into simplified PC-SAFT
The association contribution alone is not sufficient to define the full equation of state. For this, we incorporate the association theory into the wider PC-SAFT approach. In PC-SAFT molecules are modelled as chains of length m of tangentially bonded hard spheres of diameter σ with segment-segment square well attractions of depth ε. The total excess (residual) Helmholtz free energy is defined as exex hs ch at asB Aa ma a a aNk T = = + + + (36) The association contribution to the free energy is described with the theory developed in this work Eq. (19). The contributions a hs and a ch are the excess free energy of the hard sphere reference fluid and the change in free energy associated with chain formation from a fluid of hard spheres respectively. The average chain length in a mixture is given by (where x i is the mole fraction of species i ) n i ii m x m = = (37) In this work we follow the simplified approach of evaluating these contributions using the pure component Carnahan and Starling forms of a hs and the contact value of the pair correlation function g hs ( ) ( ) ( )
22 3
14 3 2; 1 ln ;1 1 hs ch hs hs a a m g g −−= = − =− − (38) where the packing fraction is given by n i i ii m d = = (39) and d i is the temperature dependent hard sphere diameter 6 ii i B d k T = − − (40) It has been demonstrated that the use of Eqns. (38) instead of the standard mixture form results in an equally capable equation of state. The use of the simplified approach in Eq. (38) is particularly useful in the current work, as it allows for a comparatively simple form of the association contribution to the chemical potential. This quantity is derived in Appendix A. The contribution to the free energy due to isotropic attractions is given by a at in Eq. (36). Gross and Sadowski developed a at using a modified Barker-Hendersen second order perturbation theory (BH2) applied to chain molecules. This completes the description of our TPT2 modification of PC-SAFT for the description of steric effects in 3C-SH alcohol models. In section V we apply the new theory to study both pure component and mixture phase equilibria of methanol and ethanol. 7 V: Application to pure methanol and ethanol
In this section we apply the new 3C-SH theory to study the phase equilibria and hydrogen bonding structure of pure methanol and ethanol. Alcohols are parameterized by 6 physically meaningful parameters: chain length m , hard sphere diameter σ, isotropic square well attraction energy ε, hydrogen bond volume κ OH , hydrogen bond energy ε OH and the blockage integral Ψ. The standard parameterization approach for SAFT theories is to adjust the model parameters to vapor pressure ( P sat ) and saturated liquid density data ( ρ L ). In this work, due to the need to obtain the blockage integral Ψ, we also include data for the heat of vaporization ( h vap ) . Specifically we adjust model parameters to minimize the objective function J given as ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) , , , , , ,1 1 1, , , hP k k k k k kn nn sat data sat theory L data L theory vap data vap theoryk k kk k ksat data L data vap data P P h hJ P h = = = − − −= + + (41) where n p , n ρ and n h is the number of vapor pressure, saturated liquid density, and heat of vaporization points respectively. Note, we have weighted the heat of vaporization data by a factor of 0.2. We consider two cases in this study. For the first case we consider the 3C-SH model which treats alcohols with two oxygen lone pair acceptor and one hydrogen donor association sites where the two acceptor sites exhibit steric hindrance (TPT2). In the second case which we call TPT1, we treat the 3C association scheme in a first order perturbation theory which does not account for steric hindrance between the acceptor sites (Ψ = 1). 8 AAD Species level m σ ε / k B ε OH / k B κ OH Ψ P sat ρ L h vap Methanol TPT2 1.914 3.0361 206.565 2362.19 0.02898 0.2 0.70% 1.71% 3.87% Methanol TPT1 1.989 3.0123 224.729 2065.50 0.02699 1 0.28% 2.57% 6.57% Ethanol TPT2 2.746 3.0303 197.754 2294.91 0.02667 0.2 0.08% 0.40% 3.80% Ethanol TPT1 2.632 3.0978 214.965 2085.08 0.02065 1 0.06% 1.20% 5.60%
Table 1:
Model parameters for methanol and ethanol and average absolute deviations (AAD). The temperature ranges used in the data regression were P sat (240 K < T < 462 K), ρ L (160 K < T < 462 K) and h vap (300 K < T < 460 K). The parameters ε / k B and ε OH / k B are in units K and σ is in units of angstroms. The resulting parameters and average absolute deviations can be found in Table 1. The best fit value for the blockage integral was found to be near Ψ = 0.2 for both alcohols. To force consistency, we then set the Ψ = 0.2 for both TPT2 parameter sets and refit the remaining parameters. A value of Ψ = 0.2 means that 80% of the allowable associated states to where both oxygen sites are bonded would be rejected due to overlap. Otherwise, the two TPT1 and TPT2 parameter sets are similar for each molecule. TPT2 gives better agreement with the ab initio calculated hydrogen bonding energy of ethanol(zero point energy corrected MP2) ε OH = 2239.4 k B K. Both TPT2 and TPT1 accurately represent vapor pressure, liquid density and h vap , although TPT2 gives a significant improvement in the representation of both liquid density and h vap . However, TPT2 and TPT1 predict substantially different liquid phase hydrogen bond structure. Figure 3 compares model predictions to molecular dynamics simulations using the OPLS-AA force field for the fraction of ethanol molecules bonded k times in a pure saturated liquid 𝜒 𝑘 . The fractions 𝜒 𝑘 are derived in TPT2 in Appendix B. 9 Figure 3:
Comparison of MD simulations (circles) and theory predictions: TPT2-solid curve, TPT1-dashed curve of the fraction of ethanol molecules bonded k times in a saturated liquid We begin the discussion of Fig. 3 with a comparison of 𝜒 and 𝜒 . As can be seen, TPT2 predicts values of 𝜒 and 𝜒 which are consistent with the MD simulation results, while TPT1 under-predicts 𝜒 and over-predicts 𝜒 . This behavior is explicable by the fact that TPT1 does not account for steric effects between the two oxygen acceptor sites. Since there is no steric hindrance in TPT1, it overpredicts the fraction 𝜒 . For a molecule to be fully bonded it must be bonded at both acceptor sites. As the sum ∑ 𝜒 𝑘 = 1 must hold, the over-prediction 𝜒 must be debited from other fractions. This then results in the under-prediction of 𝜒 in TPT1. TPT2 is in better agreement with the simulated 𝜒 at lower temperatures while TPT1 appears to be in better agreement at high temperatures. Finally, both approaches under-predict the monomer fraction data (bottom left) as compared to both the MD simulations. 0 Please note, that the definition of a “hydrogen bond” will in general be different between simulation and theory. Hence, the comparison in Fig. 3 is qualitative in nature. Qualitatively, both TPT2 and MD predict a substantially lower 𝜒 than TPT1. However, for the T = 373 K data point, the comparison is not sufficiently quantitative to say whether TPT2 or TPT1 yields a more accurate prediction of 𝜒 . Table 2 compares model predictions to simulation data for the fractions 𝜒 𝑘 of liquid methanol at T = 300 K. The simulation results are taken from Ferrando et al. who employed the AUA4 force field, which has been demonstrated to accurately represent liquid methanol structure. The TPT2 predictions accurately represent this simulation data, while TPT1 overpredicts the fraction 𝜒 and underpredicts 𝜒 for the same reasons discussed above for ethanol. Fraction MC TPT2 TPT1 𝜒 𝜒 𝜒 𝜒 Table 2:
Comparison of model predictions to Monte Carlo simulation predictions for the fraction of methanol molecules bonded k times in a saturated liquid at T = 300 K. VI: Water-ethanol mixtures
In section V it was demonstrated that TPT2 gave an improved representation of the hydrogen bonding structure of saturated liquid methanol and ethanol. In this section we focus on water-ethanol binary mixtures. However, before jumping into binary calculations, we first digress briefly to develop a benchmark water model. 1 a. Development of a benchmark water model
PC-SAFT has found wide application as a general-purpose equation of state for multi-component phase equilibria; however, it has been demonstrated by several authors that PC-SAFT does not give an entirely satisfactory representation the thermodynamics of pure water. The parameter sets for water which yield the best agreement with vapor pressures and saturated liquid densities, give a poor representation of the hydrogen bonded structure in the fluid. Recently, Marshall demonstrated that if one accounts for the transition of water to tetrahedral symmetry within the second order Barker-Hendersen perturbation theory for free energy a at , an accurate equation of state for water could be obtained which accurately represents both the thermodynamics and hydrogen bonding structure of pure water. In this work we develop a standard 4C (2 acceptor and 2 donor association sites) for water within TPT1, but we only develop the model for use over a limited temperature range for which it can act as an accurate “reference” when performing binary hydrogen bond structure and thermodynamics calculations with ethanol. When developing the water model, we include vapor pressure and liquid density data in the temperature range 293 K ≤ T ≤ 473 K. In addition, to optimize the choice of hydrogen bonding parameters, we include hydrogen bond structure data on the fraction of free OH groups ( X A ) as measured by Luck in the data regression over this same temperature range. We adjust model parameters to minimize the following objective function ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) , , , , , ,1 1 1, , , P X k k k k k knn nsat data sat theory L data L theory A data A theoryk k kk k ksat data L data A data
P P X XJ P X = = = − − −= + + (42) where n X is the number of X A data points. The results can be found in Table 3. Over this limited temperature range, PC-SAFT can accurately represent vapor pressure, saturated liquid density and the fraction of free OH groups. 2 AAD m σ (Å) ε /k b (K) ε OH /k b (K) κ OH P sat ρ L X A
1 3.063 350.624 1502.33 0.031196 0.28% 1.50% 2.10%
Table 3:
Regression results for 4C water model in the temperature range 293 K ≤ T ≤ 473 K b. Predictions of Ethanol-water hydrogen bonding
For the ethanol-water pair, the cross-association bond volume and association energy are given by the combining rules in Eq. (21). There is no binary parameter in the association contribution of the theory to tune to experimental data. Hence any calculation of binary hydrogen bonding structure is necessarily a prediction. In Fig. 4 we compare TPT1 and TPT2 predictions to MD simulations (same reference and methodology as Fig. 3 for ethanol, with water modeled with the i -AMOEBA force field) for the average number of hydrogen bonds per ethanol molecule 𝑛 𝐻𝐵 (𝑥) , where x is the mole fraction of ethanol, in an ethanol-water binary liquid mixture at T = 298 and 373 K. As can be seen, the 𝑛 𝐻𝐵 (𝑥) composition dependence of TPT2 is in much better agreement with simulation than TPT1. Most notably, TPT2 predicts much less hydrogen bonding in the ethanol dilute region, due to the inclusion of steric effects. 3 Figure 4:
Comparison of model predictions (dashed curve-TPT1, solid curve-TPT2) to MD simulation results (circles) for the average number of bonds per ethanol molecule in a water/ethanol binary mixture at T = 298 K (left) and T = 373 K (right) Figure 5:
Same as Fig. 4, except for all results scaled by the average number of hydrogen bonds in pure ethanol 𝑛 𝐻𝐵 (1) The MD calculation of 𝑛 𝐻𝐵 (𝑥) in Fig. 4 is ambiguous in the sense that one must impose on the simulation the definition of what it means to be hydrogen bonded. Hence, the absolute value of the simulated 𝑛 𝐻𝐵 (𝑥) is less telling than the composition dependence. In Fig. 5 we present the same results as Fig. 4, but we scale all results by the average number of hydrogen bonds in pure 4 ethanol 𝑛 𝐻𝐵 (1) . As can be seen, the composition dependence of TPT2 is in good agreement with the simulation results, while TPT1 clearly shows a stronger composition dependence. c. Ethanol-water phase equilibria and heat of mixing
To describe the binary phase equilibria and heat of mixing for the ethanol-water binary system we must first adjust the binary interaction parameter k ij ; which is used in the calculation of the cross square well depth ε ij in the free energy term due to isotropic attractions a at . We use standard Lorentz-Berthelot mixing rules ( ) ij i j ij k = − (43) The binary parameter k ij is adjusted to minimize the error in the description of binary vapor-liquid equilibria (VLE). Method k ij TPT1 -0.02424 TPT2 -0.06413
Table 4:
Binary parameters for the ethanol-water pair Table 4 gives regressed values of k ij using both TPT1 and TPT2, and Fig. 6 compares TPT1 and TPT2 model results to experimental data for the phase diagram at atmospheric pressure. As can be seen, both TPT1 and TPT2 are able to accurately correlate the binary VLE data, although TPT2 does give an improved result. 5 Figure 6:
Binary phase equilibria model results (TPT1-dashed cure, TPT2-solid curve) compared to experimental data (circles) for ethanol-water binary at atmospheric pressure While both TPT1 and TPT2 give a similar correlation of binary VLE, TPT2 is substantially more accurate in the prediction of the heat of mixing h mix . This can be seen in Fig. 7 which compares model predictions of h mix to experimental data at several temperatures. While neither approach is accurate at the T = 323.15 K, TPT2 gives much more accurate predictions of h mix at the higher temperatures T = 373.15 K and T = 398.15 K. For low ethanol concentrations, both TPT1 and TPT2 give similar predictions for h mix . However, for ethanol rich mixtures, TPT1 substantially under-predicts h mix , while TPT2 is in good agreement with the experimental data. This is due to the inclusion of steric effects in the TPT2 theory. The disagreement between model and experiment at T = 323.15 K occurs at a temperature which coincides with a change from h mix > 0 to h mix < 0 as temperature is decreased. It is also in this temperature regime that water is undergoing a structural transition towards tetrahedral symmetry. This transition is not included in a standard PC-SAFT model for water. Hence, the disagreement at T = 323.15 K for h mix is likely the result of the PC-SAFT water model and not of the theoretical treatment of ethanol itself. Recently, Marshall included water’s transition to 6 tetrahedral symmetry in the PC-SAFT equation of state by means of an associated reference fluid. A point of future research could be to combine the approach of Marshall, which accounts for the transition to tetrahedral symmetry of fully hydrogen bonded water, with the TPT2 alcohol model developed in this work. Figure 7:
Heat of mixing (change of enthalpy of mixing) as a function of ethanol mole fraction in a water ethanol mixture. Solid curve give TPT2 predictions, dashed curve TPT1 predictions and symbols give experimental data
7 When inserting a water molecule into an ethanol rich liquid phase, TPT1 will predict overly accessible oxygen acceptor sites on the ethanol molecule, while TPT2 correctly accounts for the fact that when one ethanol acceptor site is occupied, it partially blocks the second. This results in an increase in h mix as compared to TPT1. Figure 9 demonstrates this hydrogen bonding effect by showing the theory predicted average number of hydrogen bonds per water molecule in an ethanol / water binary mixture at a temperature of 373.15 K. As the mole fraction of ethanol increases, so does the deviation between TPT1 and TPT2. When water is dilute, TPT2 predicts substantially less water hydrogen bonding than TPT1. Figure 8:
Theory predictions of the average number of liquid phase hydrogen bonds per water molecule in an ethanol / water binary mixture at T = 373.15 K VI: Conclusions
In this paper we have extended TPT2 to account for steric effects in associating molecules. The theory is derived for a multi-component fluid where each molecule can have an arbitrary number of acceptor and donor association sites. However, for mathematical tractability, we have restricted the theory such that only a single pair of association sites are allowed to sterically hinder each other. The theory was then applied to a 3C-SH alcohol model where steric hindrance between the two oxygen acceptor sites was accounted for in the model. It was demonstrated in Fig. 2 that complete blockage could be achieved at the TPT2 level. We then paired this new association theory with the PC-SAFT equation of state to study hydrogen bonding and phase equilibria of alcohols. It was demonstrated that accounting for steric hindrance in TPT2 allows for a more capable equation of state as compared to a traditional first order approach.
Appendix A: Derivation of the association contribution to the chemical potential
In this appendix we derive the association contribution to the chemical potential in Wertheim’s multi-density formalism. In general, the association contribution to the chemical potential is given by the relation ( ) ( ) ( )( ) /ln k okas o kb c VXk T = − (44) From Eqns. (5)-(8) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ln/ 12 I o n i i i khs A Ak k i A gc V X c F = = + (45) Specializing to the 3C-SH association model 9 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ln 1ln 1 1 j n nk j j j i i i ihs A AO O Ok i j B in Oi i ihs O Ok ii O O gF X XXg X X = = = = − = − − (46) Combining Eqns. (15), (45) and (46) (47) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ln/ 1 1 1 1 12 i i io n n O Oi i i i ihs A O Ok k i ii ia O O O O X Xgc V X XX X = = = − + − + − −
Where the second sum on the right-hand side of Eq. (A4) is over the species of 3C-SH molecules only.
Appendix B: Derivation of the fraction of molecules bonded k times 𝝌 𝒌 In Wertheim’s multi-density statistical mechanics, the density of molecules bonded at the set of sites γ is given by ( )( ) ( )( ) k kk Po c = = (48) where P ( γ )is the partition of the set γ into non-empty subsets. For the 3C-SH association model, the density of molecules bonded once at any association site is given by ( )( ) ( ) ( ) k k kO Hko c c = + (49) The density bonded twice 0 ( )( ) ( ) ( ) ( ) ( ) ( ) k k k k kO H O O Oko c c c c = + + (50) And the density bonded three times ( )( ) ( ) ( ) ( ) ( ) ( ) k k k k kO H O O Hko c c c c = + (51) The fraction of molecules bonded j times are calculated through the relation ( ) ( )( ) kjkj k = (52) References: J.D. Weeks, D. Chandler, and H.C. Andersen, J. Chem. Phys. , 5422 (1971). J.A. Barker and D. Henderson, J. Chem. Phys. , 2856 (1967). H.C. Andersen, J. Chem. Phys. , 4714 (1973). H.C. Andersen, J. Chem. Phys. , 4985 (1974). L.W. Dahl and H.C. Andersen, J. Chem. Phys. , 1962 (1983). L.R. Pratt and D. Chandler, J. Chem. Phys. , 147 (1977). M.S. Wertheim, J. Stat. Phys. , 19 (1984). M.S. Wertheim, J. Stat. Phys. , 35 (1984). M.S. Wertheim, J. Stat. Phys. , 459 (1986). M.S. Wertheim, J. Stat. Phys. , 477 (1986). W.G. Chapman,
PhD Disertation (Cornell University, Ithaca, NY, 1988). W.G. Chapman, K.E. Gubbins, G. Jackson, and M. Radosz, Ind. Eng. Chem. Res. , 1709 (1990). 1 J. Gross and G. Sadowski, Ind. Eng. Chem. Res. , 5510 (2002). S. Dufal, T. Lafitte, A.J. Haslam, A. Galindo, G.N.I. Clark, C. Vega, and G. Jackson, Mol. Phys. , 948 (2015). L.F. Vega and F. Llovell, Fluid Phase Equilib. , 150 (2016). N. von Solms, M.L. Michelsen, and G.M. Kontogeorgis, Ind. Eng. Chem. Res. , 1098 (2003). G.M. Kontogeorgis, M.L. Michelsen, G.K. Folas, S. Derawi, N. von Solms, and E.H. Stenby, Ind. Eng. Chem. Res. , 4855 (2006). B.D. Marshall and W.G. Chapman, Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys. , 52307 (2013). A. Haghmoradi, L. Wang, and W.G. Chapman, Mol. Phys. , 2548 (2016). R.P. Sear and G. Jackson, Phys. Rev. E , 386 (1994). B.D. Marshall and W.G. Chapman, J. Chem. Phys. , 54902 (2013). R.P. Sear and G. Jackson, Mol. Phys. , 1033 (1994). J. Janecek and P. Paricaud, J. Phys. Chem. B , 7874 (2012). Y.V. Kalyuzhnyi, B.D. Marshall, W.G. Chapman, and P.T. Cummings, J. Chem. Phys. , (2013). B.D. Marshall and W.G. Chapman, Soft Matter , (2014). S.H. Huang and M. Radosz, Ind. Eng. Chem. Res. , 2284 (1990). N. von Solms, L. Jensen, J.L. Kofod, M.L. Michelsen, and G.M. Kontogeorgis, Fluid Phase Equilib. , 272 (2007). W.A. Fouad, L. Wang, A. Haghmoradi, D. Asthagiri, and W.G. Chapman, J. Phys. Chem. B , 3388 (2016). 2 M.S. Wertheim, J. Chem. Phys. , 7323 (1987). N.F. Carnahan and K.E. Starling, J. Chem. Phys. , 635 (1969). J. Gross and G. Sadowski, Ind. Eng. Chem. Res. , 1244 (2001). N.F. Carnahan, J. Chem. Phys. , 635 (1969). J.A. Barker and D. Henderson, J. Chem. Phys. , 4714 (1967). E.E. Fileti, P. Chaudhuri, and S. Canuto, Chem. Phys. Lett. , 494 (2004). W. Jorgesen, D. Maxwell, and J. Tirado-Rivas, J. Am. Chem. Soc. , 11225 (1996). N. Ferrando, J.-C. de Hemptinne, P. Mougin, and J.-P. Passarello, J. Phys. Chem. B , 367 (2012). X. Liang, B. Maribo-Mogensen, I. Tsivintzelis, and G.M. Kontogeorgis, Fluid Phase Equilib. , 2 (2016). G.M. Kontogeorgis, I. Tsivintzelis, N. von Solms, A. Grenner, D. Bøgh, M. Frost, A. Knage-Rasmussen, and I.G. Economou, Fluid Phase Equilib. , 219 (2010). B.D. Marshall, J. Chem. Phys. , 174104 (2017). B.D. Marshall, Ind. Eng. Chem. Res. , 4070 (2018). B.D. Marshall, Phys. Rev. E , 52602 (2017). W.A.P. Luck, Angew. Chemie Int. Ed. English , 28 (1980). L.-P. Wang, T. Head-Gordon, J.W. Ponder, P. Ren, J.D. Chodera, P.K. Eastman, T.J. Martinez, and V.S. Pande, J. Phys. Chem. B , 9956 (2013). H.-S. Lai, Y.-F. Lin, and C.-H. Tu, J. Chem. Thermodyn. , 13 (2014). J.R. Errington and P.G. Debenedetti, Nature , 318 (2001). J.. Ott, C.. Stouffer, G.. Cornett, B.. Woodfield, C. Guanquan, and J.. Christensen, J. Chem. Thermodyn. , 337 (1987). 3 J.B. Ott, G.V. Cornett, C.E. Stouffer, B.F. Woodfield, C. Guanquan, and J.J. Christensen, J. Chem. Thermodyn.18