# A Multiscale Gibbs-Helmholtz Constrained Cubic  Equation of State

Authors: Angelo Lucia
Journal: Journal of Thermodynamics
 (2010)

---

## Abstract

This paper presents a radically new approach to cubic equations of state (EOS) in which the Gibbs-Helmholtz equation is used to constrain the attraction or energy parameter,a. The resulting expressions for a(T,p) for pure components and a(T,p,x) for mixtures contain internal energy departure functions and completely avoid the need to use empirical expressions like the Soave alpha function. Our approach also provides a novel and thermodynamically rigorous mixing rule for a(T,p,x). When the internal energy departure function is computed using Monte Carlo or molecular dynamics simulations as a function of current bulk phase conditions, the resulting EOS is a multiscale equation of state. The proposed new Gibbs-Helmholtz constrained (GHC) cubic equation of state is used to predict liquid densities at high pressure and validated using experimental data from literature. Numerical results clearly show that the GHC EOS provides fast and accurate computation of liquid densities at high pressure, which are needed in the determination of gas hydrate equilibria.


---

## Body

## 1. Introduction and Motivation

Carbon storage in deep ocean sediments is one of the many technologies being explored to mitigate carbon dioxide (CO2) emissions from the consumption of carbon-based fuels (e.g., oil, coal, and even biofuels) and the associated greenhouse gas (GHG) effects that have impacted global climate change. Carbon storage has been identified by the National Academy of Engineering as one of fourteen grand challenge problems of the 21st century [1]. The feasibility of permanent deep ocean storage rests on two factors—neutral buoyancy (Brewer et al. [2], Dornan et al. [3]) and hydrate formation—and recent impact estimates indicate that the associated displacement of seawater in porous sediments by CO2 emissions from 1000 coal-fired power plants over 100 years would cause only a 1-millimeter rise in sea levels (Combs [4]).Pumping carbon dioxide in deep ocean sediments is challenging. Ocean depth gives rise to a point of neutral buoyancy and geothermal heating in the sediments produces a second, deeper point of neutral buoyancy because the density of CO2 decreases more rapidly than pore seawater with depth due to heating. Thus there is a neutral buoyancy zone (NBZ). Additionally, pumping liquid CO2 into the sediments at temperatures near 275.15 K and pressures around 30 MPa produces locally supersaturated conditions of CO2 in seawater, creating an environment that favors the formation and sustainability of a liquid carbon dioxide and CO2 hydrates thermodynamically. Thus various combinations of a liquid CO2 phase, liquid seawater and one or more hydrate phases will coexist. There are also chemical reactions in seawater that produce carbonate and bicarbonate ions as well as the presence of sodium (Na) and chloride (Cl) ions in seawater; thus strong electrolyte solution behavior must be considered. To assess and predict the short and long-term impact of storing carbon dioxide in deep ocean sediments, accurate quantitative descriptions of the growth and possible dissolution of liquid CO2 and CO2 hydrate reservoirs as well as the interactions between these reservoirs and the ocean over multiple time and length scales are needed. This, in turn, requires a concerted effort of computer modeling, experimental work, and observational data to build understanding and to validate computer models.The focus of this paper is on one aspect of these complex simulations—the rapid and accurate computation of liquid properties using equations of state (EOS) for modeling the growth and dissolution of CO2 liquid and gas hydrate reservoirs. This class of reservoir simulations can require hundreds of thousands of property, phase split, and phase equilibrium computations, and it has long been recognized that cubic equations of state generally do a poor job predicting accurate liquid densities at high pressure. Moreover, more complicated equations of state like the statistical associating fluid theory (SAFT) equation and its many variants require nested iteration for hydrogen bonding compounds and thus will not provide rapid enough computations for reservoir simulations.The novel contributions of this paper are the following:(1)
the development of a new thermodynamically rigorous framework for pure component parameters in cubic equations of state (EOS) based on using the Gibbs-Helmholtz equation as a constraint;(2)
both closed-form and integral multiscale expressions for the energy parameters in cubic EOS that directly incorporate molecular level information (specifically internal energy departure functions) obtained from Monte Carlo or molecular dynamics simulations;(3)
the development of a new mixing rule based on the Gibbs-Helmholtz equation for the prediction of phase properties and phase equilibria for mixtures;(4)
numerical testing and validation with experimental data from the open literature.Thesignificance of the modeling effort contained in this article is the development and validation of a radically new, fast, reliable, and truly predictive cubic EOS framework using a novel multiscale modeling approach that makes combined use of molecular simulations and classical thermodynamics. Specifically, the Gibbs-Helmholtz equation is used to constrain the energy parameter, a(T, p) in the van der Waals family of cubic EOS through the internal energy departure function. The proposed framework is very general, uses Monte Carlo simulations to evaluate the internal energy departure function, and can readily accommodate the usual molecular interactions for nonelectrolyte species (e.g., van der Waals forces) as well as electrostatic interactions (i.e., charge-charge, charge-dipole, dipole-dipole, charge-quadrupole, etc.). Electrostatic interactions for weak and strong electrolyte systems (or other interactions) can be directly taken into account by using an appropriate potential energy function at the molecular level. Moreover, the proposed modeling framework, when coupled with excess Gibbs free energy (GE) models for solid and hydrate phases, has the potential to (1) rapidly, reliably, and accurately determine all types of phase equilibrium involving vapor (V), liquid (L), solid ice (I), and gas hydrate (H) phases and (2) be readily incorporated into a variety of (reservoir) simulation programs.


### 1.1. Preview

Figure1 is reproduction of Figure 4 in Brewer et al. [2] who report experimental mass densities for liquid CO2 at seven distinct pressures and two temperatures (273.15 K and 283.15 K). The o’s and +’s in Figure 1 represent experimental data points while the curves shown by …. and - - - in Figure 1 are extrapolated from least squares fits of the experimental data. Also shown in Figure 1 is an estimate of the upper point of neutral buoyancy (i.e., where the mass density of CO2 and seawater are equal). Note that this occurs at 2600 m where T=275.15K and p=26.5 MPa and the mass density of seawater is 1.045 g/cm3. For comparison, mass densities of liquid CO2 at experimental conditions were calculated in two different ways for a temperature of 275.15K : (1) using the classical Soave-Redlich-Kwong (SRK) equation and (2) by the novel multiscale Gibbs-Helmholtz constrained (GHC) equation of state for pure components proposed in this work. The classical SRK equation (shown in green) does extremely poorly, predicting densities that are way too low (e.g., at a depth of 2600 m, the density of liquid CO2 predicted by SRK is 0.9627 g/cm3 (just outside the figure). On the other hand, the novel Gibbs-Helmholtz constrained equation for pure components that will be described in this paper predicts liquid CO2 densities (shown as the red curve) much closer to extrapolated experimental values (shown by - - -). More specifically, at 275.15K and 26.5 MPa, the proposed Gibbs-Helmholtz constrained equation of state predicts a liquid CO2 density of 1.0378 g/cm3, which compares favorably with the extrapolated experimental value of 1.045 g/cm3. Not surprisingly, the proposed approach also predicts an upper point of neutral buoyancy at approximately 2845 m (29.0 MPa), which is in good agreement with the experimental value of 2600 m reported by Brewer et al. [2]. Details of these computations are given in Section 4.Figure 1
Comparison of experimental and calculated liquid CO2 densities.The key conclusions that can be drawn from the numerical comparisons shown in Figure1 are the following.(1)
The classical SRK equation does a very poor job of matching liquid carbon dioxide densities at high pressure.(2)
The novel multiscale Gibbs-Helmholtz constrained (GHC) cubic equation of state for pure components proposed in this paper predicts a good approximate upper point of neutral buoyancy between pure liquid CO2 and seawater.


### 1.2. Current Needs

To be successful in modeling deep ocean sedimentary storage of CO2, any equation of state must:(1)
accurately model multiple liquid phases involving both nonelectrolyte (liquid CO2) and strong electrolyte (seawater with dissolved CO2) phase behavior,(2)
accurately predict themass density inversions that occur between seawater and liquid CO2 in the neutral buoyancy zone (i.e., both the upper and lower points of neutral buoyancy),(3)
be computationally tractable and provide a high level of computational reliability since phase equilibrium computations need to be performed over many finite elements of any reservoir simulation.(4)
be able to be coupled to models for ice and/or gas hydrate phases in order to successfully model liquid-liquid-hydrate (LLH) or other equilibrium.


## 1.1. Preview

Figure1 is reproduction of Figure 4 in Brewer et al. [2] who report experimental mass densities for liquid CO2 at seven distinct pressures and two temperatures (273.15 K and 283.15 K). The o’s and +’s in Figure 1 represent experimental data points while the curves shown by …. and - - - in Figure 1 are extrapolated from least squares fits of the experimental data. Also shown in Figure 1 is an estimate of the upper point of neutral buoyancy (i.e., where the mass density of CO2 and seawater are equal). Note that this occurs at 2600 m where T=275.15K and p=26.5 MPa and the mass density of seawater is 1.045 g/cm3. For comparison, mass densities of liquid CO2 at experimental conditions were calculated in two different ways for a temperature of 275.15K : (1) using the classical Soave-Redlich-Kwong (SRK) equation and (2) by the novel multiscale Gibbs-Helmholtz constrained (GHC) equation of state for pure components proposed in this work. The classical SRK equation (shown in green) does extremely poorly, predicting densities that are way too low (e.g., at a depth of 2600 m, the density of liquid CO2 predicted by SRK is 0.9627 g/cm3 (just outside the figure). On the other hand, the novel Gibbs-Helmholtz constrained equation for pure components that will be described in this paper predicts liquid CO2 densities (shown as the red curve) much closer to extrapolated experimental values (shown by - - -). More specifically, at 275.15K and 26.5 MPa, the proposed Gibbs-Helmholtz constrained equation of state predicts a liquid CO2 density of 1.0378 g/cm3, which compares favorably with the extrapolated experimental value of 1.045 g/cm3. Not surprisingly, the proposed approach also predicts an upper point of neutral buoyancy at approximately 2845 m (29.0 MPa), which is in good agreement with the experimental value of 2600 m reported by Brewer et al. [2]. Details of these computations are given in Section 4.Figure 1
Comparison of experimental and calculated liquid CO2 densities.The key conclusions that can be drawn from the numerical comparisons shown in Figure1 are the following.(1)
The classical SRK equation does a very poor job of matching liquid carbon dioxide densities at high pressure.(2)
The novel multiscale Gibbs-Helmholtz constrained (GHC) cubic equation of state for pure components proposed in this paper predicts a good approximate upper point of neutral buoyancy between pure liquid CO2 and seawater.


## 1.2. Current Needs

To be successful in modeling deep ocean sedimentary storage of CO2, any equation of state must:(1)
accurately model multiple liquid phases involving both nonelectrolyte (liquid CO2) and strong electrolyte (seawater with dissolved CO2) phase behavior,(2)
accurately predict themass density inversions that occur between seawater and liquid CO2 in the neutral buoyancy zone (i.e., both the upper and lower points of neutral buoyancy),(3)
be computationally tractable and provide a high level of computational reliability since phase equilibrium computations need to be performed over many finite elements of any reservoir simulation.(4)
be able to be coupled to models for ice and/or gas hydrate phases in order to successfully model liquid-liquid-hydrate (LLH) or other equilibrium.


## 2. Literature Survey

A brief survey of literature relevant to phase equilibrium in electrolyte and gas hydrate systems is presented.


### 2.1. Existing Models

Current models for multiphase equilibrium capable of modeling strong electrolyte systems include the segment-based NRTL activity coefficient (SAC) methods developed by Chen and Song [5, 6], the predictive Soave-Redlich-Kwong (PSRK) equation of Kiepe et al. [7], the electrolyte SAFT equation (Galindo et al. [8], Bezhadi et al. [9]), and a variety of very specific models (e.g., Chen and Millero [10]).   Each of these approaches has advantages and disadvantages.


### 2.2. EOS-Based Models

The PSRK equation of  Kiepe et al. [7] is a cubic equation of state that uses GE-based mixing rule for the attraction parameter in the SRK equation coupled to the Liquid Functional Activity Coefficient (LIFAC) group contribution method of Yan et al. [11]. LIFAC accounts for the influence of short (van der Waals), middle (charge-dipole), and long (charge-charge) effects of strong electrolytes on the activity coefficients of the nonelectrolyte species. That is, phase behavior for nonelectrolyte components is directly modeled using the SRK equation with a Huron-Vidal type or GE-based mixing rule (Huron and Vidal [12]), where activity coefficients are first predicted by LIFAC and then modified by iso-activity constraints to obtain salt-free estimates of activity coefficients for use in GE mixing rules. Also because of the presence of ions, a large number of parameters are needed to capture the interactions between gases and ionic groups in the PSRK electrolyte approach. Many results for gas solubility, Henry’s law constants, and pressure-composition diagrams are presented and compared to experimental data, and while the work of Kiepe et al. [7]. does address strong electrolyte systems, no studies of liquid-liquid equilibrium at conditions relevant to CO2 storage are reported. Additionally, recent work by Yoon et al. [13] shows that the PSRK model (not the PSRK electrolyte model) can be used in conjunction with van der Waals-Platteeuw theory to predict various equilibria between vapor, liquid, ice, and hydrates. However, while Yoon and Yamamoto report liquid-liquid-hydrate equilibrium for cyclopropane hydrates, they only consider nonelectrolyte systems with water.The Statistical Associating Fluid Theory (SAFT) equation has also been extended to electrolyte systems. All approaches thus far use a mean spherical approximation (MSA) for long-range Coulombic interactions and calculate solvent-ion interactions using a potential functions (e.g., square-well, Yukawa potentials). See , for example, Galindo et al. [8]. and Bezhadi et al. [9]. who use the SAFT-Variable Range (SAFT-VR) equation to study vapor-liquid equilibrium in strong electrolyte mixtures. While SAFT has proven to be a very popular approach, the iterative computation of densities has embedded in it another level of iteration for finding mole fractions of unbonded sites in mixtures that exhibit hydrogen bonding (i.e., strong electrolytes). This coupled with the fact that SAFT tends to yield high-order transcendental density functionality makes it difficult to guarantee that all density roots can be calculated reliably.


### 2.3. Activity Coefficient Models

The recently developed segment-based NRTL activity coefficient (SAC) models of Chen and Song [5, 6] have been tested on a wide array of organic electrolyte mixtures common to the pharmaceutical industry and very good results have been obtained. However, no hydrate forming gases are included in the SAC solute data, and because this approach is an activity coefficient approach, it is more suitable for low-pressure applications. Thus, in our opinion, the SAC approach is not applicable to deep ocean sedimentary CO2 storage.


### 2.4. Perspective

The complex functionality and nested iteration structure for hydrogen bonding systems makes SAFT unsuitable for use in reservoir simulation where fast and reliable phase equilibrium computations are required. On the other hand, the high-pressure environment of carbon storage and the lack of model parameters for light gases eliminate the use of the segment-based activity coefficient (SAC) approach. This leaves cubic equations of state from the van der Waals family, which we believe represent a reasonable starting point for the development of an EOS for use in reservoir simulation. Cubic EOS are a reasonable approach because it is straightforward to guarantee that all density roots at specified conditions of temperature, pressure, and composition can be found using a reliable equation solver coupled with polynomial deflation. However, while this makes them computationally tractable and provides a high level of computational reliability, the main difficulty in using cubic EOS has always been their inability to accurately predict liquid density and their limited development for strong electrolyte solutions. Inaccurate density calculations have a domino effect since liquid densities are used to compute fugacity coefficients, which in turn are used to compute other thermodynamic properties(e.g.,GE=RTlnϕM-RT∑xilnϕi). Additionally, the PSRK approach requires a very large number of parameters and has not been tested on high-pressure liquid-liquid equilibrium. Therefore one of the main focuses of the computer modeling aspects of this paper is to improve the ability of cubic EOS to predict accurate liquid densities for nonelectrolyte and electrolyte solutions. In this paper, we focus specifically on nonelectrolyte systems but do consider systems with electrostatic effects.


## 2.1. Existing Models

Current models for multiphase equilibrium capable of modeling strong electrolyte systems include the segment-based NRTL activity coefficient (SAC) methods developed by Chen and Song [5, 6], the predictive Soave-Redlich-Kwong (PSRK) equation of Kiepe et al. [7], the electrolyte SAFT equation (Galindo et al. [8], Bezhadi et al. [9]), and a variety of very specific models (e.g., Chen and Millero [10]).   Each of these approaches has advantages and disadvantages.


## 2.2. EOS-Based Models

The PSRK equation of  Kiepe et al. [7] is a cubic equation of state that uses GE-based mixing rule for the attraction parameter in the SRK equation coupled to the Liquid Functional Activity Coefficient (LIFAC) group contribution method of Yan et al. [11]. LIFAC accounts for the influence of short (van der Waals), middle (charge-dipole), and long (charge-charge) effects of strong electrolytes on the activity coefficients of the nonelectrolyte species. That is, phase behavior for nonelectrolyte components is directly modeled using the SRK equation with a Huron-Vidal type or GE-based mixing rule (Huron and Vidal [12]), where activity coefficients are first predicted by LIFAC and then modified by iso-activity constraints to obtain salt-free estimates of activity coefficients for use in GE mixing rules. Also because of the presence of ions, a large number of parameters are needed to capture the interactions between gases and ionic groups in the PSRK electrolyte approach. Many results for gas solubility, Henry’s law constants, and pressure-composition diagrams are presented and compared to experimental data, and while the work of Kiepe et al. [7]. does address strong electrolyte systems, no studies of liquid-liquid equilibrium at conditions relevant to CO2 storage are reported. Additionally, recent work by Yoon et al. [13] shows that the PSRK model (not the PSRK electrolyte model) can be used in conjunction with van der Waals-Platteeuw theory to predict various equilibria between vapor, liquid, ice, and hydrates. However, while Yoon and Yamamoto report liquid-liquid-hydrate equilibrium for cyclopropane hydrates, they only consider nonelectrolyte systems with water.The Statistical Associating Fluid Theory (SAFT) equation has also been extended to electrolyte systems. All approaches thus far use a mean spherical approximation (MSA) for long-range Coulombic interactions and calculate solvent-ion interactions using a potential functions (e.g., square-well, Yukawa potentials). See , for example, Galindo et al. [8]. and Bezhadi et al. [9]. who use the SAFT-Variable Range (SAFT-VR) equation to study vapor-liquid equilibrium in strong electrolyte mixtures. While SAFT has proven to be a very popular approach, the iterative computation of densities has embedded in it another level of iteration for finding mole fractions of unbonded sites in mixtures that exhibit hydrogen bonding (i.e., strong electrolytes). This coupled with the fact that SAFT tends to yield high-order transcendental density functionality makes it difficult to guarantee that all density roots can be calculated reliably.


## 2.3. Activity Coefficient Models

The recently developed segment-based NRTL activity coefficient (SAC) models of Chen and Song [5, 6] have been tested on a wide array of organic electrolyte mixtures common to the pharmaceutical industry and very good results have been obtained. However, no hydrate forming gases are included in the SAC solute data, and because this approach is an activity coefficient approach, it is more suitable for low-pressure applications. Thus, in our opinion, the SAC approach is not applicable to deep ocean sedimentary CO2 storage.


## 2.4. Perspective

The complex functionality and nested iteration structure for hydrogen bonding systems makes SAFT unsuitable for use in reservoir simulation where fast and reliable phase equilibrium computations are required. On the other hand, the high-pressure environment of carbon storage and the lack of model parameters for light gases eliminate the use of the segment-based activity coefficient (SAC) approach. This leaves cubic equations of state from the van der Waals family, which we believe represent a reasonable starting point for the development of an EOS for use in reservoir simulation. Cubic EOS are a reasonable approach because it is straightforward to guarantee that all density roots at specified conditions of temperature, pressure, and composition can be found using a reliable equation solver coupled with polynomial deflation. However, while this makes them computationally tractable and provides a high level of computational reliability, the main difficulty in using cubic EOS has always been their inability to accurately predict liquid density and their limited development for strong electrolyte solutions. Inaccurate density calculations have a domino effect since liquid densities are used to compute fugacity coefficients, which in turn are used to compute other thermodynamic properties(e.g.,GE=RTlnϕM-RT∑xilnϕi). Additionally, the PSRK approach requires a very large number of parameters and has not been tested on high-pressure liquid-liquid equilibrium. Therefore one of the main focuses of the computer modeling aspects of this paper is to improve the ability of cubic EOS to predict accurate liquid densities for nonelectrolyte and electrolyte solutions. In this paper, we focus specifically on nonelectrolyte systems but do consider systems with electrostatic effects.


## 3. A New Multiscale Modeling Approach to PVT Behavior of Fluids

To improve the ability of cubic EOS to predict accurate liquid densities for nonelectrolyte, we have developed a new multiscale methodology for cubic equations of state by using the Gibbs-Helmholtz equations to constrain the attraction parameter,a(T,p). The resulting multiscale modeling framework(1)
is thermodynamically rigorous,(2)
easily allows molecular level information to be included,(3)
has been developed for pure components and mixtures,


### 3.1. Pure Components

Any thermodynamic departure function is defined as the difference between a real property and an ideal gas property at the same temperature and pressure. In the material that follows, we denote departure functions using the superscriptD. From the differential form of the Gibbs free energy and the definition of fugacity, f=ϕp, it follows that G/RT=Gig/RT+GD/RT→lnf=lnp+lnϕ and that
(1)GDRT=lnϕ,
where GD and Gig denote the Gibbs free energy of departure  and ideal gas Gibbs free energy respectively, evaluated at the same temperature (T) and pressure (p), ϕ is a fugacity coefficient, and R is the gas constant. From the Gibbs-Helmholtz equation, we have that
(2)[∂(G/RT)∂T]p=[∂(Gig/RT)∂T]p+[∂(GD/RT)∂T]p=-HigRT2-HDRT2.
Therefore [∂(GD/RT)/∂T]p=-HD/RT2=[∂lnϕ/∂T]p. Using the Soave form of the Redlich-Kwong (RK) equation for a pure component given by p=RT/(V-b)-a(T)/[V(V+b)] and setting b=Vs=1/ρs, where Vs is the molar volume of the solid phase, we have that
(3)lnϕ=ln[z(1-bV)]+z-1-[a(T,p)bRT]ln[(V+b)V].
At high pressure, (∂V/∂T)p≈0 is a reasonable assumption. To see this consider any pair of isotherms in the compressed liquid region of a generic pV diagram like the one shown in Figure 2. Note that the isotherms are packed closely together and therefore (∂V/∂T)p is small. Thus we let (∂V/∂T)p=0, which gives ∂ln[z(1-b/V)]/∂T=-1/Tand∂z/∂T=-pV/RT2. Subsequent differentiation of lnϕ with respect to temperature at constant pressure yields
(4)[∂lnϕ∂T]p=-HDRT2=-1T-pVRT2-[a′bRT-abRT2]ln[V+bV].
Use of the expressions: (1) pV/RT2=pVig/RT2+pVD/RT2, (2) pVig/RT2=1/T, and (3) HD/RT2=UD/RT2+pVD/RT2 in the expression for [∂lnϕ/∂T]p and some algebraic rearrangement gives
(5)Ta′-a=[bUD-2bRT]ln[(V+b)/V]=f(T).
Equation (5) is a first-order, inhomogeneous differential equation that defines a(T) at fixed pressure and is easily solved. The general solution to the homogeneous form of (5) is given by ah(T)=cT. The particular solution to (5) can be calculated using variation of parameters or simple one-sided Green’s functions. To do this, we assume that there is a solution of the form ap(T)=w(T)T, where w(T) is a weighting parameter that must be determined. Straightforward application of variation of parameters applied to (5) gives the condition w′(T)=f(T)/T2, which yields
(6)w(T)=∫T[f(T)T2]dT.
Thus ap(T)={∫[f(T)/T2]dT}T. Differentiation of ap(T) and substitution of ap(T) and ap′(T) easily show that this form of the particular solution satisfies (5). Substituting the expression for f(T) from (5) gives
(7)ap(T)={∫T[bUD-2bRT](T2ln[(V+b)/V])dT}T.
Therefore the complete solution to (5) is
(8)a(T)={c+∫T[bUD-2bRT](T2ln[(V+b)/V])dT}T,
where the constant c can be evaluated by applying the boundary condition a(Tc)=0.42748RTc2/pc at T=Tc. That is,(9)c=a(Tc)Tc-{∫T[bUD-2bRT](T2ln[(V+b)/V])dT}=0.42748RTcpc-{∫T[bUD-2bRT](T2ln[V+b/V])dT}|T=Tc.
The foregoing analysis yields the following thermodynamically constrained integral expression for a(T):(10)a(T)={0.42748R2Tcpc-{∫T[bUD-2bRT](T2ln[(V+b)/V])dT}|T=Tc}T+{∫T[bUD-2bRT](T2ln[V+b/V])dT}T.Figure 2
pV Diagram for a typical light gas fluid.Observations
(1)Approximating b, the excluded volume, by b=Vs=1/ρs, where Vs is the solid molar volume is reasonable since the solid phase represents a “close packed” phase. The standard expression b=0.08664RTc/pc often overestimates the excluded volume. Moreover, solid density data, ρs, for many pure components are readily available.
(2)Equation (10) involves UD(T,p). The effect of pressure on a(T) is incorporated implicitly through the evaluation of the internal energy departure function and thus (10) is really an expression for a(T,p).
(3)Because a(T,p) was developed using the Gibbs-Helmholtz equation as a constraint, it is thermodynamically rigorous. We do not require the acentric factor or empirical expressions for a(T) such as the alpha function given by Soave [14]
(11)a(T)=a(Tc)α(T)=a(Tc){1+m[1-(TTc)1/2]}2,wherem=0.480-1.574ω-0.176ω2.
Other similar empirical expressions (e.g., Mathias [15], Mathias and Copeman [16], Melham [17], Twu et al. [18]) are not required either. Also, the boundary condition used to evaluate the constant c in (8) is the critical point and often represents an experimentally validated condition for most pure compounds.
(4)UD(T,p) can be readily determined using isothermal-isobaric (NTP) Monte Carlo (MC) or molecular dynamics (MD) simulations. While this requires an empirical potential energy model for the component(s) of interest, it has the flexibility to readily accommodate a variety of interactions that take place between particles (van der Waals forces, Coulombic forces, etc.). Moreover, any temperature and pressure dependence of UD(T,p) can be automatically included in the molecular simulations. We use NTP Monte Carlo simulations to measure UD(T,p).
(5)Note that the expression for a(T,p) involves V and therefore (10) is coupled to the EOS. To decouple the equation of state from (10) we invoke the high-pressure limit that V=b. Under this assumption, ln[(V+b)/V]=ln2, in which case (10) becomes
(12)a(T,p)={0.42748R2Tcpc-{∫T[bUD-2bRT](T2ln2)dT}|T=Tc}T+{∫T[bUD-2bRT](T2ln2)dT}T.
The assumption ln[(V+b)/V]=ln2 is only approximate at finite pressure but, as we show, only introduces a small error in calculated liquid density.
(6)If one further assumes that temperature functionality of UD(T,p) is a weak enough to not affect the integrals in (12), then the following closed-form solution for a(T,p) results
(13)a(T,p)={0.42748R2Tcpc+bUD(Tcln2)+2bRlnTcln2}T-bUDln2-(2bRln2)TlnT.
Note that the functionality of a(T,p) involves a linear term in temperature and a TlnT term. Also note that (13) still allows UD(T,p) to be determined from molecular simulations and thus still includes implicit pressure effects on a(T,p). Equation (13) represents a radically new expression for a(T,p) that relates energies measured at the molecular scale, UD, to the energy parameter in cubic equations of state. There is nothing even close to (13) in the open literature!.


### 3.2. Mixtures

Modeling and computations involving mixtures arealways handled using mixing rules. Our analysis for mixtures begins with the fundamental thermodynamic relationship:(14)GMDRT=lnϕM-∑xilnϕi,
where GMD is the Gibbs free energy departure function of the mixture, ϕM is the fugacity coefficient for the mixture, xi denotes the mole fraction of the ith component in the mixture, and ϕi is the fugacity coefficient for pure component i. Differentiation and application of the Gibbs-Helmholtz equation for mixtures give
(15)-HMDRT2=[∂(GMD/RT)∂T]p=(∂lnϕM∂T)p,x-∑xi(∂lnϕi∂T)p,
where HMD is the enthalpy departure function of the mixture.To develop aclosed form mixing rule for aM(T,p,x) we do the following.(1)
Differentiatelnϕi=ln[zi(1-bi/Vi)]+zi-1-[ai(T)/biRT]ln[(Vi+bi)/Vi] with respect to T at constant p.(2)
LetbM=∑xibi, where each pure component excluded volume bi=1/ρis, where ρis is the density of pure solid component i.(3)
DifferentiatelnϕM=ln[zM(1-bM/VM)]+zM-1-[aM(T)/bMRT]ln[(VM+bM)/VM] with respect to T at constant p and x.(4)
Use the approximationsVi=bi and VM=bM, which imply ln[(Vi+bi)/Vi]=ln2 and ln[(VM+bM)/VM]=ln2, respectively.These assumptions lead to the differential equation(16)TaM′-aM=f(T,p,x)=bM{UMD/ln2-∑xi(pViD)+∑(xi/bi)[Tai′-ai]}ln2,
where UMD(T,p,x) is the internal energy departure function of the mixture, and where the functions ai(T,p) and their first derivatives are known and given by (13) applied to pure component i. Moreover, assumptions (4) and (5) show that Tai′-ai=[biUiD-2biRT]/ln2, and since ViD must be computed to calculate lnϕi, given UMD and UiD, the entire right hand side of (16) is well defined.Equation (16) can be solved in exactly the same way that (5) was solved. This yields the following solution or mixing rule for aM(17)aM(T,p,x)={c+∫T[f(T,x)T2]dT}T=cT-bM{[UMD-∑xiUiD]+∑xi[pViD+2RlnT]}ln2.
Note that there are three terms on the right-hand side of (17). The last two terms in the mixing rule are quadratic in composition (since they involve the product of bM(∑xi)), the last term explicitly involves the pressure, and the second term involves the quantity UMD(T,p,x), which accounts for temperature, pressure and composition dependence in an implicit manner. Note also that there is no need to use empirical relationships such as the geometric combining rule to define the coefficients of the cross composition terms (i.e., the xixj terms). Rather all terms including the cross composition terms in the proposed mixing rule are defined in terms of physically meaning quantities (i.e., bi’s, Ui’s, pViD’s, etc.). Moreover, the cross composition terms are not assumed to be symmetric.All that remains is to establish a boundary condition foraM(T,p,x) to evaluate the constant in (17). To do this we use Kay’s rules [19] to calculate mixture critical properties TcM and pcM and assume that aM(TcM,pcM)=0.42748R2TcM2/pcM for any mixture of composition x. This gives
(18)c=aM(TcM,pcM,x)Tc+bM{[UMD-∑xiUiD]+∑xi[pcMViD+2RlnTcM]}(Tcln2)
and the following expression for a(T,p,x)(19)aM(T,p,x)={aM(TcM,pcM,x)TcM+[bMTcMln2]×{[UMD-∑xiUiD]+∑xi[pcMViD+2RTcMlnTcM]}aM(TcM,pcM,x)TcM+[bMTcMln2]}T-[bMln2]{[UMD-∑xiUiD]+∑xi[pViD+2RTlnT]}.Observations
(1)The proposed mixing rule for aM(T,p,x) is thermodynamically rigorous since it has been derived using the Gibbs-Helmholtz equation as a constraint for both the mixture and for each pure component.
(2)aM(T,p,x) has the same temperature dependence for mixtures as for pure components. It involves a linear term in T as well as TlnT dependence.
(3)There is explicit pressure dependence in the mixing rule through the term pVi.
(4)This new mixing rule contains a number of terms that are quadratic in composition as well as internal energy departure functions for the mixture and pure components. The departure functions account for temperature, pressure, and composition dependence of aM in an implicit manner.
(5)The proposed mixing rule is a multiscale mixing rule since it can easily incorporate UMD and UiD from molecular simulations.


## 3.1. Pure Components

Any thermodynamic departure function is defined as the difference between a real property and an ideal gas property at the same temperature and pressure. In the material that follows, we denote departure functions using the superscriptD. From the differential form of the Gibbs free energy and the definition of fugacity, f=ϕp, it follows that G/RT=Gig/RT+GD/RT→lnf=lnp+lnϕ and that
(1)GDRT=lnϕ,
where GD and Gig denote the Gibbs free energy of departure  and ideal gas Gibbs free energy respectively, evaluated at the same temperature (T) and pressure (p), ϕ is a fugacity coefficient, and R is the gas constant. From the Gibbs-Helmholtz equation, we have that
(2)[∂(G/RT)∂T]p=[∂(Gig/RT)∂T]p+[∂(GD/RT)∂T]p=-HigRT2-HDRT2.
Therefore [∂(GD/RT)/∂T]p=-HD/RT2=[∂lnϕ/∂T]p. Using the Soave form of the Redlich-Kwong (RK) equation for a pure component given by p=RT/(V-b)-a(T)/[V(V+b)] and setting b=Vs=1/ρs, where Vs is the molar volume of the solid phase, we have that
(3)lnϕ=ln[z(1-bV)]+z-1-[a(T,p)bRT]ln[(V+b)V].
At high pressure, (∂V/∂T)p≈0 is a reasonable assumption. To see this consider any pair of isotherms in the compressed liquid region of a generic pV diagram like the one shown in Figure 2. Note that the isotherms are packed closely together and therefore (∂V/∂T)p is small. Thus we let (∂V/∂T)p=0, which gives ∂ln[z(1-b/V)]/∂T=-1/Tand∂z/∂T=-pV/RT2. Subsequent differentiation of lnϕ with respect to temperature at constant pressure yields
(4)[∂lnϕ∂T]p=-HDRT2=-1T-pVRT2-[a′bRT-abRT2]ln[V+bV].
Use of the expressions: (1) pV/RT2=pVig/RT2+pVD/RT2, (2) pVig/RT2=1/T, and (3) HD/RT2=UD/RT2+pVD/RT2 in the expression for [∂lnϕ/∂T]p and some algebraic rearrangement gives
(5)Ta′-a=[bUD-2bRT]ln[(V+b)/V]=f(T).
Equation (5) is a first-order, inhomogeneous differential equation that defines a(T) at fixed pressure and is easily solved. The general solution to the homogeneous form of (5) is given by ah(T)=cT. The particular solution to (5) can be calculated using variation of parameters or simple one-sided Green’s functions. To do this, we assume that there is a solution of the form ap(T)=w(T)T, where w(T) is a weighting parameter that must be determined. Straightforward application of variation of parameters applied to (5) gives the condition w′(T)=f(T)/T2, which yields
(6)w(T)=∫T[f(T)T2]dT.
Thus ap(T)={∫[f(T)/T2]dT}T. Differentiation of ap(T) and substitution of ap(T) and ap′(T) easily show that this form of the particular solution satisfies (5). Substituting the expression for f(T) from (5) gives
(7)ap(T)={∫T[bUD-2bRT](T2ln[(V+b)/V])dT}T.
Therefore the complete solution to (5) is
(8)a(T)={c+∫T[bUD-2bRT](T2ln[(V+b)/V])dT}T,
where the constant c can be evaluated by applying the boundary condition a(Tc)=0.42748RTc2/pc at T=Tc. That is,(9)c=a(Tc)Tc-{∫T[bUD-2bRT](T2ln[(V+b)/V])dT}=0.42748RTcpc-{∫T[bUD-2bRT](T2ln[V+b/V])dT}|T=Tc.
The foregoing analysis yields the following thermodynamically constrained integral expression for a(T):(10)a(T)={0.42748R2Tcpc-{∫T[bUD-2bRT](T2ln[(V+b)/V])dT}|T=Tc}T+{∫T[bUD-2bRT](T2ln[V+b/V])dT}T.Figure 2
pV Diagram for a typical light gas fluid.Observations
(1)Approximating b, the excluded volume, by b=Vs=1/ρs, where Vs is the solid molar volume is reasonable since the solid phase represents a “close packed” phase. The standard expression b=0.08664RTc/pc often overestimates the excluded volume. Moreover, solid density data, ρs, for many pure components are readily available.
(2)Equation (10) involves UD(T,p). The effect of pressure on a(T) is incorporated implicitly through the evaluation of the internal energy departure function and thus (10) is really an expression for a(T,p).
(3)Because a(T,p) was developed using the Gibbs-Helmholtz equation as a constraint, it is thermodynamically rigorous. We do not require the acentric factor or empirical expressions for a(T) such as the alpha function given by Soave [14]
(11)a(T)=a(Tc)α(T)=a(Tc){1+m[1-(TTc)1/2]}2,wherem=0.480-1.574ω-0.176ω2.
Other similar empirical expressions (e.g., Mathias [15], Mathias and Copeman [16], Melham [17], Twu et al. [18]) are not required either. Also, the boundary condition used to evaluate the constant c in (8) is the critical point and often represents an experimentally validated condition for most pure compounds.
(4)UD(T,p) can be readily determined using isothermal-isobaric (NTP) Monte Carlo (MC) or molecular dynamics (MD) simulations. While this requires an empirical potential energy model for the component(s) of interest, it has the flexibility to readily accommodate a variety of interactions that take place between particles (van der Waals forces, Coulombic forces, etc.). Moreover, any temperature and pressure dependence of UD(T,p) can be automatically included in the molecular simulations. We use NTP Monte Carlo simulations to measure UD(T,p).
(5)Note that the expression for a(T,p) involves V and therefore (10) is coupled to the EOS. To decouple the equation of state from (10) we invoke the high-pressure limit that V=b. Under this assumption, ln[(V+b)/V]=ln2, in which case (10) becomes
(12)a(T,p)={0.42748R2Tcpc-{∫T[bUD-2bRT](T2ln2)dT}|T=Tc}T+{∫T[bUD-2bRT](T2ln2)dT}T.
The assumption ln[(V+b)/V]=ln2 is only approximate at finite pressure but, as we show, only introduces a small error in calculated liquid density.
(6)If one further assumes that temperature functionality of UD(T,p) is a weak enough to not affect the integrals in (12), then the following closed-form solution for a(T,p) results
(13)a(T,p)={0.42748R2Tcpc+bUD(Tcln2)+2bRlnTcln2}T-bUDln2-(2bRln2)TlnT.
Note that the functionality of a(T,p) involves a linear term in temperature and a TlnT term. Also note that (13) still allows UD(T,p) to be determined from molecular simulations and thus still includes implicit pressure effects on a(T,p). Equation (13) represents a radically new expression for a(T,p) that relates energies measured at the molecular scale, UD, to the energy parameter in cubic equations of state. There is nothing even close to (13) in the open literature!.


## 3.2. Mixtures

Modeling and computations involving mixtures arealways handled using mixing rules. Our analysis for mixtures begins with the fundamental thermodynamic relationship:(14)GMDRT=lnϕM-∑xilnϕi,
where GMD is the Gibbs free energy departure function of the mixture, ϕM is the fugacity coefficient for the mixture, xi denotes the mole fraction of the ith component in the mixture, and ϕi is the fugacity coefficient for pure component i. Differentiation and application of the Gibbs-Helmholtz equation for mixtures give
(15)-HMDRT2=[∂(GMD/RT)∂T]p=(∂lnϕM∂T)p,x-∑xi(∂lnϕi∂T)p,
where HMD is the enthalpy departure function of the mixture.To develop aclosed form mixing rule for aM(T,p,x) we do the following.(1)
Differentiatelnϕi=ln[zi(1-bi/Vi)]+zi-1-[ai(T)/biRT]ln[(Vi+bi)/Vi] with respect to T at constant p.(2)
LetbM=∑xibi, where each pure component excluded volume bi=1/ρis, where ρis is the density of pure solid component i.(3)
DifferentiatelnϕM=ln[zM(1-bM/VM)]+zM-1-[aM(T)/bMRT]ln[(VM+bM)/VM] with respect to T at constant p and x.(4)
Use the approximationsVi=bi and VM=bM, which imply ln[(Vi+bi)/Vi]=ln2 and ln[(VM+bM)/VM]=ln2, respectively.These assumptions lead to the differential equation(16)TaM′-aM=f(T,p,x)=bM{UMD/ln2-∑xi(pViD)+∑(xi/bi)[Tai′-ai]}ln2,
where UMD(T,p,x) is the internal energy departure function of the mixture, and where the functions ai(T,p) and their first derivatives are known and given by (13) applied to pure component i. Moreover, assumptions (4) and (5) show that Tai′-ai=[biUiD-2biRT]/ln2, and since ViD must be computed to calculate lnϕi, given UMD and UiD, the entire right hand side of (16) is well defined.Equation (16) can be solved in exactly the same way that (5) was solved. This yields the following solution or mixing rule for aM(17)aM(T,p,x)={c+∫T[f(T,x)T2]dT}T=cT-bM{[UMD-∑xiUiD]+∑xi[pViD+2RlnT]}ln2.
Note that there are three terms on the right-hand side of (17). The last two terms in the mixing rule are quadratic in composition (since they involve the product of bM(∑xi)), the last term explicitly involves the pressure, and the second term involves the quantity UMD(T,p,x), which accounts for temperature, pressure and composition dependence in an implicit manner. Note also that there is no need to use empirical relationships such as the geometric combining rule to define the coefficients of the cross composition terms (i.e., the xixj terms). Rather all terms including the cross composition terms in the proposed mixing rule are defined in terms of physically meaning quantities (i.e., bi’s, Ui’s, pViD’s, etc.). Moreover, the cross composition terms are not assumed to be symmetric.All that remains is to establish a boundary condition foraM(T,p,x) to evaluate the constant in (17). To do this we use Kay’s rules [19] to calculate mixture critical properties TcM and pcM and assume that aM(TcM,pcM)=0.42748R2TcM2/pcM for any mixture of composition x. This gives
(18)c=aM(TcM,pcM,x)Tc+bM{[UMD-∑xiUiD]+∑xi[pcMViD+2RlnTcM]}(Tcln2)
and the following expression for a(T,p,x)(19)aM(T,p,x)={aM(TcM,pcM,x)TcM+[bMTcMln2]×{[UMD-∑xiUiD]+∑xi[pcMViD+2RTcMlnTcM]}aM(TcM,pcM,x)TcM+[bMTcMln2]}T-[bMln2]{[UMD-∑xiUiD]+∑xi[pViD+2RTlnT]}.Observations
(1)The proposed mixing rule for aM(T,p,x) is thermodynamically rigorous since it has been derived using the Gibbs-Helmholtz equation as a constraint for both the mixture and for each pure component.
(2)aM(T,p,x) has the same temperature dependence for mixtures as for pure components. It involves a linear term in T as well as TlnT dependence.
(3)There is explicit pressure dependence in the mixing rule through the term pVi.
(4)This new mixing rule contains a number of terms that are quadratic in composition as well as internal energy departure functions for the mixture and pure components. The departure functions account for temperature, pressure, and composition dependence of aM in an implicit manner.
(5)The proposed mixing rule is a multiscale mixing rule since it can easily incorporate UMD and UiD from molecular simulations.


## 4. Numerical Results and Comparisons of Modeling and Experiments

In this section, numerical results for liquid densities using the new proposed framework are presented and compared with experimental data. Specifically, liquid densities for pure carbon dioxide, pure water, and CO2-water mixtures are calculated using (13) and (19) and compared to high pressure experimental density data. All bulk phase property computations and molecular simulations were performed in double precision arithmetic on a Dell High Precision 670 workstation using the Lahey-Fijitsu (LF95) compiler.


### 4.1. Numerical Results in Figure1

The results for the classical SRK equation shown in Figure1 were computed using Tc=304.20K, pc=73.80 bar, and ω=0.224, which yield b=0.08664RTc/pc=29.693 cm3/mol and a(T)=a(Tc)α(T)=0.42748R2Tc2/pc, where α(T) is the Soave alpha function. The numerical results for the Gibbs-Helmholtz constrained cubic EOS shown in Figure 1 were determined using b=Vs=28.169 cm3/mol and (13) with UD=-1.2×10-5 cm3bar/mol, which was determined by performing a set of 25 NTP Monte Carlo simulations at low pressure (N=500,T=275.15K,p=10bar). The potential energy model for CO2 consisted of Lennard-Jones forces and electrostatic forces. Table 1 shows the calculated densities for both SRK and the G-H constrained (GHC) EOS.Table 1
Calculated liquid CO2 densities at 275.15 K.
p (MPa)ρSRK(g/cm3)ρGHC(g/cm3)20.20.93231.018122.80.94561.026525.30.95741.034328.00.96911.042130.50.97921.049035.50.99751.061840.61.01421.0737The numerical values shown in Table1 are those that have been plotted in Figure 1.


### 4.2. Other Validated Liquid CO2 Density Predictions

Table2 shows the liquid CO2 experimental data of Magee and Ely [20].Table 2
Experimental molar density data for pure carbon dioxide*.
T(K)p(bar)ρ(mol/cm3)T(K)p(bar)ρ(mol/cm3)233.383112.9260.025953254.81191.3330.023936235.585142.8550.025959257.119115.7250.023924237.777172.4620.025924259.454140.3190.023912239.962201.8110.025910261.761164.5170.023900242.142230.9280.025896264.058188.5290.023889244.308259.7110.025882266.348212.3840.023877246.474288.3470.025868268.627236.0400.023866248.634316.6180.025854270.910159.6500.023854273.188283.1320.023843271.34983.8990.022139285.807204.0780.022081275.311116.9360.022123290.551243.2350.022062277.704136.8560.022113292.916262.7020.022052280.089156.6830.022103295.278282.0940.022043281.032164.5120.022100297.644301.4730.022033283.421184.3230.022090*Reproduced from Magee and Ely [20, runs 900 through 1100 in Table V, page 1178–1179].It is important to note that runs 900 through 1100 are theonly runs that measure sample mass using a gravimetric method for fixed experimental apparatus volume; all other runs determine molar density using a 32-term extended Benedict-Rubin-Webb equation. See Magee and Ely [20, page 1167].Figure3 shows a comparison of the liquid molar densities calculated using the SRK equation and the proposed multiscale Gibbs-Helmholtz constrained equation of state with the experimental data of Magee and Ely [20]. Note that the proposed multiscale method using the expression for a(T,p) given by (13) does an exceptionally good job of matching experimental liquid CO2 density data over a wide range of temperatures and pressures while the Soave expression of a(T)=a(Tc)α(T) consistently underestimates the liquid molar density.Figure 3
Comparison of calculated and experimental liquid molar densities forCO2.Figure 4
Comparison of specific volumes for water at 273.15 K. experimental data are reproduced from Kell and Whalley [21].


### 4.3. Validated Numerical Results for Densities of Compressed Water

We have calculated liquid densities for water at pressures between 5.30 and 1026.39 bars at various temperatures using the SRK and the proposed GHC equation and compared those numerical results to the experimental data of Kell and Whalley [21]. Since water is a polar compound, we have included the Peneloux [22] volume translation given by the expression
(20)V′SRK=VSRK-0.40768(RTcpc)[0.29441-pcVcRTc]
in the computations using the SRK equation, where VSRK is the volume predicted by the traditional SRK equation and V′SRK is the SRK predicted volume modified by the Peneloux volume translation. The critical properties for water used in these calculations were Tc=647.37K, pc=221.20bar, Vc=56cm3/mol, and ω=0.345. Note that the value of b calculated from critical properties is 21.082 cm3/mol and represents a very poor estimate of the excluded volume.For the GHC equation, thesolid volume, Vs, of very high density amorphous ice was used to estimate the excluded volume. This value, Vs=14.286 cm3/mol, gives a much better estimate of the excluded volume. The calculated values of UD ranged from -0.03 to -0.035 cm3bar/mol. Comparisons of specific volumes at 273.15 K are shown in Table 3 and Figure 4 in cm3/g, which is the form in which Kell and Whalley [21] reported their experimental results.Table 3
Experimental and calculated specific volumes of water at 273.15 K*.
p(bar)Vexp(cm3/g)VSRK(cm3/g)VGHC(cm3/g)p(bar)Vexp(cm3/g)VSRK(cm3/g)VGHC(cm3/g)5.300.9998820.9454221.001958268.570.9870330.9415550.98965812.600.9995120.9453721.001615326.950.9843290.9407390.98632220.000.9991380.9452501.001236385.220.9816810.9399220.98394627.400.9984620.9451271.000875443.490.9790740.9391110.98165142.110.9980250.9448891.000136501.850.9765260.9383820.98025856.810.9972820.9446710.999416560.110.9740140.9375500.97649271.510.9965500.9444610.998698618.470.9715600.9367660.97440286.310.9958130.9442170.998967676.710.9691420.9360160.972390101.010.9950880.9440050.998267735.060.9667700.9352710969725115.710.9943620.9437940.997550793.290.9644420.9345220.967833130.520.9936370.9435830.996852851.630.9621500.9338040.965999151.890.9925940.9432500.994888909.850.9599040.9330880.963553210.280.98997820.9424000.992223968.080.9576990.9323720.9618181026.390.9555290.9316880.959509*Experimental data reproduced from Kell and Whalley [21, Table 4, page 586]. V′SRK includes Peneloux volume translation (V′SRK=VSRK-0.40768((RTc)/pc)[0.29441-(pcVc)/RTc]).From Table3, it is easily seen that the SRK equation with the Peneloux [22]  volume translation does poorly at matching the experimental data. However, this is not surprising since the Peneloux volume translation term, 0.40768(RTc/pc)[0.29441-pcVc/RTc], is an empirical correction. For water, the volume translation has a value of 0.354111 cm3/g or 6.374 cm3/mol, overcompensates for polarity, and produces specific volumes that are roughly 1 cm3/g smaller than they should be. On the other hand, without the Peneloux volume translation the specific volumes calculated by the SRK equation are extremely poor (V=~1.29 cm3/g). In contrast, the GHC equation does an excellent job of matching experimental specific volumes. Moreover, we believe that the excellent match provided by the GHC equation is due to the fact that it makes combined use of molecular scale information and a rigorous classical relationship to determine a(T,p).


### 4.4. Validated Density Calculations for CO2-H2O Mixtures

To validate the mixing rule given by (19) we compared numerical density calculations for CO2-H2O mixtures with the high-pressure experimental data of Teng et al. [23] at 278 K. For the SRK calculations we included the Peneloux volume translation for water. For the GHC equation, we used one fixed set of values of pure component internal energy departure functions: UMD=-0.085×106 cm3bar/mol, UCO2D=-0.24×106 cm3bar/mol, and UH2OD=-0.035×106 cm3bar/mol in order to illustrate that repeated and costly molecular simulations are not required for good phase property predictions. The comparison is shown in Table 4.Table 4
Comparison of densities of water-rich liquid for various CO2-H2O mixtures at 278 K*.
p(bar)xCO2ρexp(kg/m3)ρ′SRK(kg/m3)ρGHC(kg/m3)64.40.02931018.101075.051016.9898.70.03081019.771076.971018.61147.70.03201020.631078.581020.92196.80.03111022.011080.291023.22245.80.03411023.331081.921025.49294.90.03491025.331083.391027.71*Experimental data reproduced from Teng et al. [23, Table 2, page 1306].ρ′SRK includes the Peneloux correction for water.Again, the GHC equation clearly outperforms the SRK equation with the Peneloux translation.


### 4.5. Reliability and Computational Speed

At the equation of state level, the liquid density calculations presented in Section4
are 100% reliable when the global terrain methodology of Lucia and Feng [24] and polynomial deflation are used. Moreover, phase density computations using the terrain approach for either pure components or mixtures require roughly 0.01 s to find all three roots. At the molecular scale, NTP Monte Carlo calculations require substantial computational resources. If MC simulations are to be performed each time the temperature, pressure, and/or composition change, then significant computational costs may result. On the other hand, if these molecular scale calculations are performed at nominal conditions of temperature, pressure and composition and interpolated or extrapolated, then considerable reductions in computational overhead will result with little error at the bulk phase length scale.


## 4.1. Numerical Results in Figure1

The results for the classical SRK equation shown in Figure1 were computed using Tc=304.20K, pc=73.80 bar, and ω=0.224, which yield b=0.08664RTc/pc=29.693 cm3/mol and a(T)=a(Tc)α(T)=0.42748R2Tc2/pc, where α(T) is the Soave alpha function. The numerical results for the Gibbs-Helmholtz constrained cubic EOS shown in Figure 1 were determined using b=Vs=28.169 cm3/mol and (13) with UD=-1.2×10-5 cm3bar/mol, which was determined by performing a set of 25 NTP Monte Carlo simulations at low pressure (N=500,T=275.15K,p=10bar). The potential energy model for CO2 consisted of Lennard-Jones forces and electrostatic forces. Table 1 shows the calculated densities for both SRK and the G-H constrained (GHC) EOS.Table 1
Calculated liquid CO2 densities at 275.15 K.
p (MPa)ρSRK(g/cm3)ρGHC(g/cm3)20.20.93231.018122.80.94561.026525.30.95741.034328.00.96911.042130.50.97921.049035.50.99751.061840.61.01421.0737The numerical values shown in Table1 are those that have been plotted in Figure 1.


## 4.2. Other Validated Liquid CO2 Density Predictions

Table2 shows the liquid CO2 experimental data of Magee and Ely [20].Table 2
Experimental molar density data for pure carbon dioxide*.
T(K)p(bar)ρ(mol/cm3)T(K)p(bar)ρ(mol/cm3)233.383112.9260.025953254.81191.3330.023936235.585142.8550.025959257.119115.7250.023924237.777172.4620.025924259.454140.3190.023912239.962201.8110.025910261.761164.5170.023900242.142230.9280.025896264.058188.5290.023889244.308259.7110.025882266.348212.3840.023877246.474288.3470.025868268.627236.0400.023866248.634316.6180.025854270.910159.6500.023854273.188283.1320.023843271.34983.8990.022139285.807204.0780.022081275.311116.9360.022123290.551243.2350.022062277.704136.8560.022113292.916262.7020.022052280.089156.6830.022103295.278282.0940.022043281.032164.5120.022100297.644301.4730.022033283.421184.3230.022090*Reproduced from Magee and Ely [20, runs 900 through 1100 in Table V, page 1178–1179].It is important to note that runs 900 through 1100 are theonly runs that measure sample mass using a gravimetric method for fixed experimental apparatus volume; all other runs determine molar density using a 32-term extended Benedict-Rubin-Webb equation. See Magee and Ely [20, page 1167].Figure3 shows a comparison of the liquid molar densities calculated using the SRK equation and the proposed multiscale Gibbs-Helmholtz constrained equation of state with the experimental data of Magee and Ely [20]. Note that the proposed multiscale method using the expression for a(T,p) given by (13) does an exceptionally good job of matching experimental liquid CO2 density data over a wide range of temperatures and pressures while the Soave expression of a(T)=a(Tc)α(T) consistently underestimates the liquid molar density.Figure 3
Comparison of calculated and experimental liquid molar densities forCO2.Figure 4
Comparison of specific volumes for water at 273.15 K. experimental data are reproduced from Kell and Whalley [21].


## 4.3. Validated Numerical Results for Densities of Compressed Water

We have calculated liquid densities for water at pressures between 5.30 and 1026.39 bars at various temperatures using the SRK and the proposed GHC equation and compared those numerical results to the experimental data of Kell and Whalley [21]. Since water is a polar compound, we have included the Peneloux [22] volume translation given by the expression
(20)V′SRK=VSRK-0.40768(RTcpc)[0.29441-pcVcRTc]
in the computations using the SRK equation, where VSRK is the volume predicted by the traditional SRK equation and V′SRK is the SRK predicted volume modified by the Peneloux volume translation. The critical properties for water used in these calculations were Tc=647.37K, pc=221.20bar, Vc=56cm3/mol, and ω=0.345. Note that the value of b calculated from critical properties is 21.082 cm3/mol and represents a very poor estimate of the excluded volume.For the GHC equation, thesolid volume, Vs, of very high density amorphous ice was used to estimate the excluded volume. This value, Vs=14.286 cm3/mol, gives a much better estimate of the excluded volume. The calculated values of UD ranged from -0.03 to -0.035 cm3bar/mol. Comparisons of specific volumes at 273.15 K are shown in Table 3 and Figure 4 in cm3/g, which is the form in which Kell and Whalley [21] reported their experimental results.Table 3
Experimental and calculated specific volumes of water at 273.15 K*.
p(bar)Vexp(cm3/g)VSRK(cm3/g)VGHC(cm3/g)p(bar)Vexp(cm3/g)VSRK(cm3/g)VGHC(cm3/g)5.300.9998820.9454221.001958268.570.9870330.9415550.98965812.600.9995120.9453721.001615326.950.9843290.9407390.98632220.000.9991380.9452501.001236385.220.9816810.9399220.98394627.400.9984620.9451271.000875443.490.9790740.9391110.98165142.110.9980250.9448891.000136501.850.9765260.9383820.98025856.810.9972820.9446710.999416560.110.9740140.9375500.97649271.510.9965500.9444610.998698618.470.9715600.9367660.97440286.310.9958130.9442170.998967676.710.9691420.9360160.972390101.010.9950880.9440050.998267735.060.9667700.9352710969725115.710.9943620.9437940.997550793.290.9644420.9345220.967833130.520.9936370.9435830.996852851.630.9621500.9338040.965999151.890.9925940.9432500.994888909.850.9599040.9330880.963553210.280.98997820.9424000.992223968.080.9576990.9323720.9618181026.390.9555290.9316880.959509*Experimental data reproduced from Kell and Whalley [21, Table 4, page 586]. V′SRK includes Peneloux volume translation (V′SRK=VSRK-0.40768((RTc)/pc)[0.29441-(pcVc)/RTc]).From Table3, it is easily seen that the SRK equation with the Peneloux [22]  volume translation does poorly at matching the experimental data. However, this is not surprising since the Peneloux volume translation term, 0.40768(RTc/pc)[0.29441-pcVc/RTc], is an empirical correction. For water, the volume translation has a value of 0.354111 cm3/g or 6.374 cm3/mol, overcompensates for polarity, and produces specific volumes that are roughly 1 cm3/g smaller than they should be. On the other hand, without the Peneloux volume translation the specific volumes calculated by the SRK equation are extremely poor (V=~1.29 cm3/g). In contrast, the GHC equation does an excellent job of matching experimental specific volumes. Moreover, we believe that the excellent match provided by the GHC equation is due to the fact that it makes combined use of molecular scale information and a rigorous classical relationship to determine a(T,p).


## 4.4. Validated Density Calculations for CO2-H2O Mixtures

To validate the mixing rule given by (19) we compared numerical density calculations for CO2-H2O mixtures with the high-pressure experimental data of Teng et al. [23] at 278 K. For the SRK calculations we included the Peneloux volume translation for water. For the GHC equation, we used one fixed set of values of pure component internal energy departure functions: UMD=-0.085×106 cm3bar/mol, UCO2D=-0.24×106 cm3bar/mol, and UH2OD=-0.035×106 cm3bar/mol in order to illustrate that repeated and costly molecular simulations are not required for good phase property predictions. The comparison is shown in Table 4.Table 4
Comparison of densities of water-rich liquid for various CO2-H2O mixtures at 278 K*.
p(bar)xCO2ρexp(kg/m3)ρ′SRK(kg/m3)ρGHC(kg/m3)64.40.02931018.101075.051016.9898.70.03081019.771076.971018.61147.70.03201020.631078.581020.92196.80.03111022.011080.291023.22245.80.03411023.331081.921025.49294.90.03491025.331083.391027.71*Experimental data reproduced from Teng et al. [23, Table 2, page 1306].ρ′SRK includes the Peneloux correction for water.Again, the GHC equation clearly outperforms the SRK equation with the Peneloux translation.


## 4.5. Reliability and Computational Speed

At the equation of state level, the liquid density calculations presented in Section4
are 100% reliable when the global terrain methodology of Lucia and Feng [24] and polynomial deflation are used. Moreover, phase density computations using the terrain approach for either pure components or mixtures require roughly 0.01 s to find all three roots. At the molecular scale, NTP Monte Carlo calculations require substantial computational resources. If MC simulations are to be performed each time the temperature, pressure, and/or composition change, then significant computational costs may result. On the other hand, if these molecular scale calculations are performed at nominal conditions of temperature, pressure and composition and interpolated or extrapolated, then considerable reductions in computational overhead will result with little error at the bulk phase length scale.


## 5. Conclusions

A radically new cubic equation of state approach was presented in which the excluded volume parameter,b, was approximated using solid molar volumes and new expressions for the energy parameter, a, were derived by using the Gibbs-Helmholtz equation to constrain the value of a. The resulting expressions for a(T,p) and a(T,p,x) necessarily (1)
are thermodynamically rigorous,(2)
avoid the need for empirical correlations like the alpha function of Soave,(3)
have temperature functionality with a linear term, and aTlnT term.(4)
Involve internal energy departure functions, and thus make the new approach a multiscale cubic equation of state approach.Moreover, the mixing rule or expression fora(T,p,x) has a pressure explicit term clearly showing that the energy parameter is pressure dependent, albeit weak. The resulting new cubic equation approach is called the Gibbs-Helmholtz constrained (GHC) equation and is truly predictive.The GHC equation was compared to the classical SRK equation and validated using experimental high-pressure liquid density data for CO2, water, and CO2-H2O mixtures. Numerical testing clearly shows that the GHC equation clearly uniformly outperforms the SRK equation and matches experimental liquid density data exceptionally well.Future work with the GHC equation will involve further numerical testing and experimental validation of liquid densities for a wide variety of mixtures as well as predictions of phase equilibrium for systems involving multiple liquid phases and gas hydrates.



---

*Source: 238365-2009-11-17.xml*