A highly accurate equation of state (EOS) for chain molecules formed from spherical segments interacting through Mie potentials (i.e., a generalized Lennard-Jones form with variable repulsive and attractive exponents) is presented. The quality of the theoretical description of the vapour-liquid equilibria (coexistence densities and vapour pressures) and the second-derivative thermophysical properties (heat capacities, isobaric thermal expansivities, and speed of sound) are critically assessed by comparison with molecular simulation and with experimental data of representative real substances. Our new EOS represents a notable improvement with respect to previous versions of the statistical associating fluid theory for variable range interactions (SAFT-VR) of the generic Mie form. The approach makes rigorous use of the Barker and Henderson high-temperature perturbation expansion up to third order in the free energy of the monomer Mie system. The radial distribution function of the reference monomer fluid, which is a prerequisite for the representation of the properties of the fluid of Mie chains within a Wertheim first-order thermodynamic perturbation theory (TPT1), is calculated from a second-order expansion. The resulting SAFT-VR Mie EOS can now be applied to molecular fluids characterized by a broad range of interactions spanning from soft to very repulsive and short-ranged Mie potentials. A good representation of the corresponding molecular-simulation data is achieved for model monomer and chain fluids. When applied to the particular case of the ubiquitous Lennard-Jones potential, our rigorous description of the thermodynamic properties is of equivalent quality to that obtained with the empirical EOSs for LJ monomer (EOS of Johnson et al.) and LJ chain (soft-SAFT) fluids. A key feature of our reformulated SAFT-VR approach is the greatly enhanced accuracy in the near-critical region for chain molecules. This attribute, combined with the accurate modeling of second-derivative properties, allows for a much improved global representation of the thermodynamic properties and fluid-phase equilibria of pure fluids and their mixtures.

The goal of this contribution is the development of an accurate equation of state (EOS) for fluids based on a generic and versatile force field that captures the key features of the interactions for a broad range of real substances. The gradual unmasking of the role played by intermolecular forces in determining the thermophysical properties of fluids and the marked sensitivity of the observed behaviour on the nature and range of the interaction (the softness/hardness of the repulsions or the extent and range of the electrostatic attractions) has a long and illustrious history.1 In order to place our generic approach in the appropriate context we first summarize the key findings, as the subtleties are not always fully appreciated or have even been forgotten altogether. This is particularly true in the development of EOSs where the precise form of the underlying force field is often deemed as secondary to the theoretical formulation. A detailed knowledge of the various contributions to the intermolecular potential is central to the accurate description of molecular systems; a robust link between the force field and the thermodynamic properties is provided in our paper.

Much of the studies on “corpuscular” interactions in the late seventeenth and eighteenth centuries centered around the ideas of Newton, who had already considered attractive forces between particles that varied as an (undetermined) inverse power of their separation in his Principia Mathematica,2 but these were largely speculative with minimal empirical validation. The existence and nature of repulsive forces was even more disputed and poorly understood, though Newton and his supporters Rowning and Bošković recognized that repulsions must also be active at short distances to account for the incompressibility of dense fluids and materials; some of the sketches of the interactions between “Boscovichian” particles3 look strikingly similar to the potential functions that we are familiar with today. Perhaps the exception in terms of firm experimental confirmation was the discovery by Priestley,4 Cavendish,5 and Coulomb6 that charged bodies (no matter what their physical dimension) interact through an inverse-square force law, i.e., a potential function of the gravitational form which involves the reciprocal of the separation, 1/r. Young7 and Laplace8 were convinced at the turn of the nineteenth century that one had to take into account both repulsive and attractive interactions between particles to explain the observed cohesion and capillarity in fluids, but it took much longer to establish the precise forms of these interactions.

In trying to interpret the attractive potential between particles in terms of electrostatic forces, Mosotti9 was one of the first to develop an explicit potential for the interactions in dense electrically neutral fluids as an exponentially damped version of the Coulombic interaction (i.e., ua(r) ∝ exp (−λ/r)/r with λ as a system-dependent potential-range parameter) based on Poisson-Laplace electrostatics; this form of force field has since made an appearance in a number of physical scenarios and is now commonly referred to as the Yukawa10 potential.11 By 1857 Clausius12 fully recognized and explicitly stated that molecules repel each other at short relative separations and attract each other over a moderate range, starting a passionate quest amongst his contemporaries for the specific mathematical form of the overall interaction potential that best described the transport properties of dilute gases. In his seminal paper on the dynamical theory of gases, Maxwell13 extended his earlier work on perfectly elastic hard spheres (which borrowed heavily from the kinetic theory of Clausius) to particles interacting through a pairwise repulsive force taken to be inversely proportional to a power of the separation between their centres of mass (i.e., the repulsive analogue of the Newtonian attractive form). Maxwell13 noted that the treatment of molecules as spherical hard cores leads to an inadequate description of the diffusion of gases. By examining his experimental data for the viscosity of air at different temperatures (where he observed a direct proportionality with temperature, which is now known to be erroneous), Maxwell13 was able to deduce that, according to his theory, the exponent characterizing the repulsive force had to take the value of five to enable one to reproduce the observations, i.e., fr(r) ∝ 1/r5; the repulsive pair interaction potential between these Maxwellian particles would thus be of the form ur(r) ∝ 1/r4. Though in modern terms this may appear to be an overly soft interaction, Maxwell pointed out that he is representing the gas as “small bodies or groups of molecules repelling one another”; at the end of our paper we will return to the interesting consequences of this type of coarse graining of the interactions between collections of molecules in terms of single spherical cores. In work along the same vein, Boltzmann14 came to the contrasting conclusion that a completely different intermolecular force law (including a generic attractive contribution) could be used to reproduce the experimental gas viscosity.

The emergence of these ideas enabled van der Waals15 to develop his equation of state for fluids which, for the first time, allowed one to represent vapour-liquid equilibria (VLE). Within the van der Waalsian model, molecules are represented as perfectly hard-spherical cores (which maintain the low compressibility of the dense fluid) with additional long-range attractive interactions (which stabilize the condensed phase); however, the specific form of the attractive forces remained unspecified, the contribution of the interactions to the thermodynamic properties of the fluid being described instead with an integrated energy which conceals the precise algebraic dependence of the interaction on the molecular separation. Van der Waals toyed with the idea of a potential of the Yukawa form but did not pursue it further. In a series of contributions to the Philosophical Magazine spanning the last two decades of the nineteenth century, Sutherland16–18 followed van der Waals’ lead and coupled a hard-sphere model with pair attractive forces which he postulated would vary as the inverse fourth power of the separation (i.e., a pair potential of the form ua(r) ∝ 1/r3, obtained rather mysteriously as the second term of an expansion of the gravitational force which also depended on the particle's mass). With this inverse-fourth-power law for the attractions, Sutherland16–18 was able in large part to reconcile the existing experimental data for the thermophysical properties such as the equation of state,16 the vapour-liquid interfacial tension,17 and the viscosity of gases.18 

A significant leap in unveiling the empirical form of the repulsive interactions between particles was made by Mie19 in 1903, almost half a century after the pioneering works of Clausius12 and Maxwell.13 Like Maxwell, Mie represented the repulsive pair potential between molecules in terms of a generic power law with a term of the form

$u_{r} (r) \propto 1/r^{\lambda _r}$
ur(r)1/rλr⁠, where λr is the exponent characterizing the nature of the repulsions, but he now also included the attractive interactions implicitly with an integrated van der Waals energy. Mie's treatment of the potential is generic in the sense that it can be used to describe explicitly repulsive interactions of varying softness/hardness, and attractive interactions of different range as described through the magnitude of the van der Waals constant. We will make use of this particular versatility of the Mie potential form to represent the interactions between molecular segments later in our paper. By examining the equation of state and thermal properties of metals (namely, mercury, copper, nickel, and zinc), Mie showed that in order to provide a consistent description of the experimental values of the compressibility, a steeper repulsive exponent was required (in his case λr ∼ 5 instead of the Maxwellian value of λr = 4), something that Boltzmann had also become increasingly convinced about in his later years. The key finding that thermodynamic derivative properties such as the compressibility are particularly sensitive to the steepness of the repulsive pair potential is a principal theme running through our current work. It is not widely recognized, however, that Grüneisen20,21 appears to have been the first to write down an explicit expression for the overall interaction potential between pairs of particles (again for metallic systems) as the sum of a generic Maxwellian repulsive contribution
$u_{r} (r) \propto 1/r^{\lambda _r}$
ur(r)1/rλr
and a generic Newtonian attractive contribution
$u_{a} (r) \propto 1/r^{\lambda _a}$
ua(r)1/rλa
, where λa characterizes the variable range of the attractive interaction; this form of interaction is often ascribed to Mie or, even more-commonly referred to as the generalized Lennard-Jones potential. In his work Keesom22,23 focussed on the form of the attractive potential by examining the second virial coefficient of molecules modeled as rigid hard spheres with Newtonian attractive interactions; this type of force field is now generally referred to as the model of Sutherland in recognition of the Australian's seminal contribution. Keesom22,23 found that the description of the temperature dependence of the second virial coefficient (another thermodynamic derivative property) of argon was relatively insensitive to the value of the attractive exponent (at least in the range λa = 3 to 5), with a marginal preference for λa = 4 as the optimal value.

Lennard-Jones24–26 is now widely credited as having explicitly expressed the full pair potential in the form

$u(r)\break =C_r/r^{\lambda _r} -C_a/r^{\lambda _a}$
u(r)=Cr/rλrCa/rλa (Cr and Ca representing the corresponding repulsive and attractive constants) which, as has already been mentioned, is the generic expression that had been reported by Grüneisen twelve years earlier. In his paper of 1924 on the temperature dependence of the viscosity of gases, Lennard-Jones (Jones at the time)24 assumed an interaction potential with a variable λr and λa = 2, making this choice for the attractive exponent to render the triple integration for the collision dynamics tractable; he was able to describe the experimental viscosity of argon gas over a broad temperature range for repulsive exponents spanning
$\lambda _r=6\frac{2}{3}$
λr=623
to 20. The representation of the second virial coefficient of argon was the focus of a subsequent well-known paper,25 in which Lennard-Jones used the more-generic form of the potential; he concluded that values of
$\lambda _r=13\frac{1}{3}$
λr=1313
and λa = 4 provided the best description of the equation of state at low density, though again no unique set of exponents could really be favoured over another. The particle diameter that he obtained from the viscosity and virial coefficients of argon were found to be equivalent to within 10%. By the time that Lennard-Jones26 delivered his lecture to the Physical Society in 1931, he had settled on the ubiquitous LJ potential form which bears his name with exponents of λr = 12 and λa = 6, the latter value being consistent with the findings of the budding discipline of quantum mechanics.

With the advent of quantum mechanics in the early part of the twentieth century, a more rigorous and detailed description of the electrostatic interactions between particles became possible. The repulsive interactions could now be seen as the result of the Pauli exclusion principle between electronic clouds (or the exchange interaction discovered by Heisenberg and Dirac), though the precise dependence on the interparticle separation was still difficult to describe algebraically; a compact exponential form (exp (−r0), with σ0 representing a size parameter related to the Bohr radius) can be used for the repulsive potential resulting from the overlap of the electron density of s orbitals but only for the simplest atoms (cf. the papers of Slater et al.27,28). On the other hand, Wang29 had shown in 1927 that the attractive interactions (dispersion forces) between atoms are the result of fluctuating electric dipoles, which when treated at the quantum-mechanical level give rise to a contribution proportional to 1/r6; this form of attractive potential was later confirmed by London and Eisenschitz,30 and is now often referred to as the London dispersion interaction in recognition of London's accessible theoretical description.31 The fluctuations of the electrons in larger molecules also leads to quadrupoles and higher multipoles which give rise to a sum of contributions in powers of λa = 6, 8, 10, etc.32,33 It is therefore apparent that in order to retain a simple yet rigorous quantum-mechanical form of the pair interaction between molecules one could combine an exponential repulsive dependence with a sum of attractive terms in

$1/r^{\lambda _a}$
1/rλa⁠; in the particular case of an attractive contribution of the London form this corresponds to the Buckingham exp-6 potential.28,34,35 The reader is directed to the monograph by Stone36 for an up-to-date account of the different types of intermolecular multipolar interactions. Before we conclude our brief overview of the important quantum-mechanical insights into the development of classical force fields we should also mention that the square-well (SW) pair potential (or more precisely its quantum-mechanical analogue the “particle-in-the-box”) is also rooted in the birth of the wave mechanics of Schrödinger, Fermi, Dirac, and others. The popularity of the square-well potential in the description of the interactions between classical particles stems from its very simple mathematical form, which allows for exact and tractable statistical mechanical calculations (e.g., the algebraic evaluation of the second and third virial coefficients37,38).

By the 1950s the nature of the intermolecular pair potential was well understood, as was the sensitivity of the thermophysical properties of fluids to the precise form of the force field.39 Molecular-beam experiments allowed for a direct determination of the pair interaction, and both the exponential form of the repulsions at short distances and the London form of the attractions at moderate distances were confirmed.40 It was also well recognised that a universal form of pair potential (e.g., with fixed exponents characterizing the repulsive and attractive interactions) would not be generally appropriate to represent the intermolecular forces of different substances.39,40 As Hirshfelder et al.39 pointed out in 1954: “The only conclusion that can be drawn [from the body of existing work] is that a potential function more realistic than the Lennard-Jones (12-6) potential must be chosen if consistent information about the intermolecular forces is to be derived from the various physical properties.” The fact that the description of the intermolecular interactions in fluids and materials with the Lennard-Jones (12-6) potential is still so prevalent in modern theoretical and computer-simulation studies is rather surprising. The reason can in part be ascribed to the larger number of intermolecular parameters that need to be determined for potentials of a more-generic form. This is hampered by the lack of an accurate algebraic theory for the thermodynamic properties of fluids described with more-realistic models, which would facilitate the efficient development of the force fields for broad classes of systems.

In our paper we describe the fluid-phase equilibria and thermodynamic derivative properties of fluids of varying complexity by retaining a generic representation of the softness/hardness of the repulsive interactions and the range of the attractive interactions. The Mie (generalized Lennard-Jones) pair potential affords the appropriate versatility in this regard making it an ideal choice of force field; a power law of variable range can be used to provide a good representation of the more-realistic exponential form of the repulsive interactions, while still retaining mathematical tractability. A prerequisite to assessing the adequacy of the Mie potential in representing the interactions in real substances is an accurate description of the thermodynamics and structure of the fluid.

The development of an accurate molecular equation of state for the thermodynamic properties of fluids of associating chain molecules is an important goal as it provides a framework for the representation of the behaviour of a variety of real fluids and their mixtures with minimal computational cost. Moreover, a rigorous theory of fluid systems is invaluable in enabling an in-depth molecular-level understanding of the complex phenomena often encountered in engineering applications. The statistical associating fluid theory (SAFT),41,42 which is based on the first-order thermodynamic perturbation theory (TPT1) of Wertheim,43–48 is among the most successful of the modern algebraic descriptions of associating chain fluids. The literature on the development and use of the SAFT EOS is too extensive to be enumerated here; the reader is directed to excellent reviews which provide a comprehensive description of the theory and its successful application to the phase equilibria of complex fluids ranging from small associating compounds to multi-component systems including polymers, surfactants, micelles, liquid crystals, asphaltenes, and electrolyte solutions.49–55 The main feature of the SAFT approach is that the thermodynamic and structural properties of a fluid of associating chain molecules are described in terms of those of a corresponding monomer (reference) system comprising the segments that form the molecule. One must first define explicitly the intermolecular potential of the reference fluid, typically an isotropic potential, and then evaluate the system's free energy and structure via the radial distribution function (RDF). There are now many different embodiments of the generic SAFT approach, differing mainly in the interaction potential used to model the reference fluid. The earlier versions of the SAFT EOS were based on a reference system of hard spheres. In more-recent versions of the theory, reference fluids represented with the square-well,56–60 Lennard-Jones,61–67 Yukawa,56,68,69 and Sutherland56 potentials have been developed; examples of the more-popular versions of SAFT in current use include the variable-range SAFT-VR,56,70 the LJ soft-SAFT,66,67 and the perturbed-chain PC-SAFT71 EOSs. A detailed comparison of PC-SAFT71 and SAFT-VR (the latter implemented for both the SW56 and LJ65 potentials) has been reported,72–74 and it was demonstrated that EOSs based on a hard-core or LJ pair potential do not generally provide one with an adequate representation of thermodynamic second-derivative properties such as the isothermal compressibility or the speed of sound in the compressed liquid phase; the description of these properties is found to be no better than that obtained with EOSs of the conventional van der Waals form such as the Peng-Robinson75 (PR) or Soave-Redlich-Kwong76 (SRK) cubic relations. In order to overcome these limitations a generic description of the interactions can be introduced,72 using the Mie potential to describe the monomer-monomer interactions between the segments making up the molecules. This potential can be expressed as19–21,25,26

\begin{equation}u^{Mie}(r)=\mathcal {C} \epsilon \left(\left(\frac{\sigma }{r} \right)^{\lambda _{r}}-\left(\frac{\sigma }{r} \right)^{\lambda _{a}} \right)\end{equation}
uMie(r)=Cεσrλrσrλa
(1)

with

\begin{equation}\mathcal {C}=\frac{\lambda _{r}}{\lambda _{r}-\lambda _{a}}\left(\frac{\lambda _{r}}{\lambda _{a}} \right)^{\frac{\lambda _{a}}{\lambda _{r}-\lambda _{a}}},\end{equation}
C=λrλrλaλrλaλaλrλa,
(2)

where r is the distance between the spherical segments, σ is the position at which the potential is zero (segment diameter), ε is the potential depth, and λa and λr are the attractive and repulsive exponents, respectively; the pre-factor

$\mathcal {C}$
C ensures that the minimum of the potential is at −ε regardless of the values of the exponents. In early developments within the SAFT-VR framework56,65 expressions were derived for the EOS and RDF in the particular case of Lennard-Jones chains (SAFT-VR LJC). To extend further the predictive capabilities of the SAFT-VR treatment for the more-general case of Mie potentials, a description of the RDF at contact of the monomer Mie reference fluid was then proposed;72 here we will refer to the resulting EOS as SAFT-VR Mie 2006 to avoid any confusion with our current work. The SAFT-VR Mie 2006 EOS has been shown to provide a significantly improved description of the dense liquid phase, especially for the second-derivative properties, while still retaining an accurate representation of the vapour-liquid equilibria.72,73 As a particular example of the marked improvement that can be achieved in the description of thermodynamic second-derivative properties by using a Mie interaction with a variable repulsive exponent to control the softness/hardness of the potential, the speed of sound (a stereotypical derivative property which is related to the compressibility and heat capacity of the fluid) of selected n-alkanes calculated with SAFT-VR Mie 200672 is compared in Table I to the corresponding description with SAFT-VR SW, SAFT-VR LJC, and PC-SAFT, as well as with the PR and SRK cubic EOSs. The previously reported values of the molecular parameters for each approach are adopted in this assessment; a more-complete comparison of the EOSs (also including soft-SAFT) with parameters estimated from thermodynamic derivative properties as well as VLE will be made in Sec. V B. The added flexibility of varying the range of the repulsive interactions is clearly the key to a reliable description of the derivative properties of the dense fluid, in addition to providing an accurate and versatile platform for the representation of the fluid-phase equilibria. The adequacy of the description of the density and thermodynamic derivative properties of model fluids with the SAFT-VR Mie 2006 EOS has been assessed by comparison with molecular-dynamics (MD) simulations of spherical molecules interacting via λr-6 Mie potentials with λr = 8, 10, and 12; the theory is found to provide a good description of the simulated thermodynamic properties of Mie fluids, with the greatest deviations seen for the heat capacity.77 An assessment of the thermodynamic properties of model chain molecules formed from Mie segments with the SAFT-VR Mie EOS developed in 200672 has not been made until now, and a detailed analysis will be presented later in our current paper. It is important to realise that the quality of the theoretical representation will depend not only on the accuracy of the approximations inherent in the Wertheim TPT1 description of the chain contribution but also on that of the underlying theory for the EOS and contact value of the RDF of the reference monomeric Mie fluid. Such a comparison will therefore represent a stringent test of the approximations used at each level of the development of our new SAFT-VR Mie EOS.

Table I.

Speed of sound u of selected n-alkanes in the condensed-liquid phase. Comparison of experimental data with the description of equations of state including the Soave-Redlich-Kwong76 (SRK), Peng-Robinson75 (PR), PC-SAFT,71 SAFT-VR SW,56 SAFT-VR LJC,65 and the SAFT-VR Mie approach developed in 2006.72 The average absolute deviation AAD(%) is computed as

${\rm AAD}(\%)=\frac{1}{N_p}\sum _{i=1}^{N_p} | \frac{X_i^{{{\rm exp}}}-X_i^{{{\rm calc}}}}{X_i^{{{\rm exp}}}}| \times 100$
AAD (%)=1Npi=1Np|Xi exp Xi calc Xi exp |×100⁠, where
$X_i^{{{\rm exp}}}$
Xi exp
and
$X_i^{{{\rm calc}}}$
Xi calc
are the experimental and computed properties, respectively, and Np is the number of experimental data points used in the calculation.

 AAD(%), u
SubstanceSRKPRPC-SAFTSAFT-VR SWSAFT-VR LJCSAFT-VR Mie 2006
n-hexane 14.1 13.0 10.3 4.7 15.9 2.2 
n-heptane 16.4 15.7 11.4 6.5 17.7 2.1 
n-tridecane 38.2 38.7 15.4 11.4 22.7 2.4 
n-tetradecane 38.6 38.5 16.0 18.2 23.0 2.5 
n-pentadecane 38.1 38.3 16.1 12.6 23.3 2.7 
n-heptadecane 46.7 47.5 16.6 12.2 23.6 2.6 
 AAD(%), u
SubstanceSRKPRPC-SAFTSAFT-VR SWSAFT-VR LJCSAFT-VR Mie 2006
n-hexane 14.1 13.0 10.3 4.7 15.9 2.2 
n-heptane 16.4 15.7 11.4 6.5 17.7 2.1 
n-tridecane 38.2 38.7 15.4 11.4 22.7 2.4 
n-tetradecane 38.6 38.5 16.0 18.2 23.0 2.5 
n-pentadecane 38.1 38.3 16.1 12.6 23.3 2.7 
n-heptadecane 46.7 47.5 16.6 12.2 23.6 2.6 

An accurate algebraic expression for the EOS of simple monomeric fluids can be developed by using an empirical or semi-empirical scheme based on exact computer-simulation data. This has proven to be a very successful approach in the case of the LJ fluid as exemplified by the work of Johnson et al.78 or Kolafa and Nezbeda79 who made use of multiparametric relations to correlate the simulation data over a wide range of densities and temperatures. As far as chains formed from LJ segments are concerned, a marked sensitivity of the Wertheim TPT1 description to the precise representation used for the RDF of the reference momoner LJ fluid has been demonstrated.80 The EOS of Johnson et al.78 has been coupled with an accurate empirical description of simulation data for the contact value of the RDF of LJ fluids to develop a high-fidelity EOS for LJ chains;63 this EOS for the LJC system is the precursor to the subsequent SAFT approaches for associating LJ and LJC fluids81,82 including the soft-SAFT approach.66,67 Despite the accuracy of the resulting description with these approaches, it should be pointed out that semi-empirical expressions for the EOS and RDF of this type are restricted to a specific intermolecular potential model (the LJ (12-6) potential in this particular case), and the applicability of the methodology is, of course, limited to the temperature and density ranges of the computer-simulation data used in the parameterization.

A more-promising route to the description of the equation of state and structural properties of a generic reference fluid is through a rigorous statistical mechanical theory. Among these methodologies, thermodynamic perturbation theories are more suitable than integral-equation approaches for the development of an algebraic formulation. The perturbation theories of Weeks, Chandler, and Andersen (WCA),83 and of Barker and Henderson (BH)84,85 are the most-ubiquitous approaches in this context. It has been shown,80,86 however, that the WCA recipe does not provide the most-accurate description of the properties of chain molecules within a SAFT-like treatment. In our current paper we will show that the perturbation approach used to develop the SAFT-VR Mie EOS in 200672 is also somewhat inaccurate in its description of the vapour-liquid equilibrium of long LJ chains. Here, we propose a new formulation for both the EOS and RDF of the monomeric Mie fluid. The novel approach is based on a rigorous application of the Barker and Henderson84,85 high-temperature perturbation expansion of the Helmholtz free energy of the Mie fluid up to third order. In addition, a significantly improved expression for the RDF of Mie monomers at contact is presented, following a self-consistent description of the pressure from the Clausius virial theorem and the density derivative of the free energy. In the special case of LJ fluids, our new generic SAFT-VR Mie EOS will be seen to provide a description of similar fidelity to that of the empirical EOS of Johnson et al.78 (cf. Sec. V A 1).

In order to assess the adequacy of our new SAFT-VR Mie EOS, the calculated thermodynamic properties are first compared with the corresponding computer-simulation data for fluids of spherical Mie particles, for chains comprising up to 100 LJ segments, and for chains of Mie λr-6 segments. The association contribution of the SAFT-VR Mie EOS is then also revisited in order to capture the effect of placing the bonding sites at different distances from the centre of the reference core. This turns out to be a useful additional feature of the new theory in the development of more-realistic models of the hydrogen-bonding interaction. A generic analytical expression of the association contribution to the free energy for an arbitrary geometry of the bonding sites is presented. Its accuracy is evaluated by comparison with the available molecular-simulation data for selected thermodynamic properties. As a final demonstration of the enhanced capabilities that our new SAFT-VR Mie EOS has to offer, the theory is used to describe the vapour-liquid equilibria and second-derivative properties for representative real substances where we demonstrate the importance and versatility of using a variable repulsive range with Mie potentials of the λr-6 form. The accurate theoretical description allows for the size and energy intermolecular parameters and the repulsive exponents to be estimated with confidence from the bulk thermodynamic data of the different systems.

In this section we derive the Helmholtz free energy and radial distribution function at contact of Mie fluids using the BH perturbation theory for soft-core potentials. Different approximations for the EOS and RDF are available and the corresponding description is compared with that obtained from our new methodology. A comparison is also made with computer-simulation data for a wide range of values of both the repulsive and attractive exponents λr and λa of the Mie potential. Based on this analysis, we propose simple and compact expressions for the EOS and RDF by making use of the mean-value theorem (MVT). The detailed implementation of these expressions in the development of the SAFT-VR Mie EOS for associating chains of Mie segments is described in Sec. IV.

1. New SAFT-VR Mie formulation

The starting point of the BH thermodynamic perturbation theory is to decompose the Mie potential given in Eq. (1) into a sum of a reference u0(r) and a perturbation u1(r) term:84 

\begin{equation}u^{Mie}(r)=u_{0}(r)+u_{1}(r),\end{equation}
uMie(r)=u0(r)+u1(r),
(3)

with

\begin{equation}u_{0}(r)=\left\lbrace \begin{array}{c @{} c}u^{Mie}(r) &\hspace{28.45274pt} r \le \sigma \\[4pt]0 &\hspace{28.45274pt} r > \sigma, \\[4pt]\end{array} \right.\end{equation}
u0(r)=uMie(r)rσ0r>σ,
(4)
\begin{equation}u_{1}(r)=\left\lbrace \begin{array}{c @{} c}0 &\hspace{28.45274pt} r \le \sigma \\[4pt]u^{Mie}(r) &\hspace{28.45274pt} r > \sigma . \\[4pt]\end{array} \right.\end{equation}
u1(r)=0rσuMie(r)r>σ.
(5)

Barker and Henderson84 showed that it is then possible to write the residual Helmholtz free energy ares = A/(NskT) as

\begin{equation}a^{res} = \sum _{n=0}^{\infty }(\beta )^{n}a_{n} ,\end{equation}
ares=n=0(β)nan,
(6)

where a0 is the Helmholtz free energy of the reference fluid, a1, a2, …, an are the first-order, second-order, and nth-order perturbation terms, respectively, β = 1/(kT) (with k representing the Boltzmann constant, and T the absolute temperature), and Ns is the number of particles (spherical segments).

The reference (repulsive) system defined by Eq. (4) is not convenient in practice because its Helmholtz free energy a0 and RDF g0(r) are not generally known. However, this reference system can be mapped onto a system of effective hard spheres, whose diameter d remains to be determined.85 In order to evaluate this diameter, an expansion of the Helmholtz free energy has to be carried out with respect to two coupling parameters (α and γ) which vary the steepness and depth of the potential in such a way that α = γ = 1 corresponds to the original Mie potential uMie(r) and α = γ = 0 to the effective hard-sphere potential uHS(r); the reference potential u0(r) is obtained in the limit of γ = 0 and α = 1. Barker and Henderson84,85 demonstrated that one can nullify the first-order term of the Helmholtz free energy in the α-expansion about a hard-sphere system of diameter d when

\begin{equation}d=\int _{0}^{\sigma } \big(1-{\rm exp} (-\beta u^{Mie}(r) )\big){\rm d}r.\end{equation}
d=0σ1 exp (βuMie(r))dr.
(7)

All the thermodynamic functions of the reference system as well as its RDF can therefore be approximated by those of this effective hard-sphere system. For instance, a0 and g0(r) may be obtained from

\begin{equation}a_0 \approx a^{HS},\end{equation}
a0aHS,
(8)
\begin{equation}g_0(r) \approx g^{HS}_d(r),\end{equation}
g0(r)gdHS(r),
(9)

where aHS and

$g^{HS}_d(r)$
gdHS(r) are the free energy and radial distribution function of the hard-sphere fluid of effective diameter d. The integral involved in the evaluation of d can be computed using an appropriate numerical integration method, e.g., one based on Gaussian quadrature.87 It follows that the Helmholtz free energy of a soft-core potential can be written as a series of the form

\begin{equation}a^{res}=a^{HS}+\beta a_1+(\beta )^{2}a_2+(\beta )^{3}a_3+ \cdots .\end{equation}
ares=aHS+βa1+(β)2a2+(β)3a3+.
(10)

The reference hard-sphere free energy aHS is calculated by integrating the well-known Carnahan and Starling88 EOS, aHS = ∫(PHS/(ρskT) − 1)/ρss,

\begin{equation}a^{HS}=\frac{4\eta -3\eta ^{2}}{(1-\eta )^{2}} ,\end{equation}
aHS=4η3η2(1η)2,
(11)

where η = ρsπd3/6 is the packing fraction of the reference hard-sphere system, and ρs = Ns/V is the segment density of the Mie fluid (with V the volume of the system).

In our current work, we propose a truncation of the series at third order instead of the more-widely used second-order expansion. This is motivated by the fact that for large values of the repulsive and/or attractive exponents, the vapour-liquid critical point can be located below T* = kT/ε = 1; in such cases truncating the high-temperature expansion at second order generally results in a significant overshoot of the critical temperature.89 It will be shown in Sec. V A that the incorporation of a third-order term considerably improves the description of the VLE in the vicinity of the critical region.

The first-order perturbation term a1 (mean-attractive energy) is calculated from

\begin{equation}a_{1}=2\pi \rho _s\int _{\sigma }^{\infty }g^{HS}_d(r)u_{1}(r)r^{2}{\rm d}r .\end{equation}
a1=2πρsσgdHS(r)u1(r)r2dr.
(12)

The key quantity is the reference hard-sphere RDF

$g^{HS}_d(r)$
gdHS(r)⁠, which we obtain by solving the Ornstein-Zernike integral equation using the reference hypernetted chain (RHNC) closure and the Malijevsky and Labik90 formula for the hard-sphere bridge function. A comparison of the mean-attractive energy obtained in this way with the corresponding molecular-simulation data is shown in Fig. 1. At this stage, we emphasize that Eq. (12) is slightly different to the approximations made in the previous work on the SAFT-VR Mie EOS of 200672 where

\begin{equation}u_{1}(r)\approx u_{1}^{approx}(r)=\left\lbrace \begin{array}{c @{} c}0 &\hspace{28.45274pt} r \le d \\[4pt]\mathcal {C}\epsilon \left(\left(\frac{d}{r} \right)^{\lambda _{r}} -\left(\frac{d}{r} \right)^{\lambda _{a}}\right) &\hspace{28.45274pt} r > d, \\[4pt]\end{array} \right.\end{equation}
u1(r)u1approx(r)=0rdCεdrλrdrλar>d,
(13)

and

\begin{equation}a_{1}\approx 2\pi \rho _s\int _{d}^{\infty }g^{HS}_d(r)u_{1}^{approx}(r)r^{2}{\rm d}r.\end{equation}
a12πρsdgdHS(r)u1approx(r)r2dr.
(14)

Note that though the latter approach is not a formally correct application of the BH theory, it is a useful approximation to derive a simple and compact algebraic expression for a1 when the SAFT-VR mapping56 is used. The corresponding description for a1 using Eq. (14) with the Malijevsky and Labik90 RHNC recipe for

$g_d^{HS}(r)$
gdHS(r) is also shown in Fig. 1(a). It can be seen that the rigorous application of the BH theory leads to a more-accurate representation of a1. In order to develop an EOS for use in engineering applications an algebraic expression for a1 is highly advantageous; such an expression will be presented in Secs. III A 2 and III A 3. The approximations inherent in Eqs. (13) and (14) have also been employed recently by Betancourt-Cardenas et al.91 to derive the perturbation terms in the Helmholtz energy expansion of the LJ fluid; this explains in part why these authors find discrepancies between their simulation data for the perturbation terms and the original values of Smith et al.92 

FIG. 1.

The density

$\rho_s^\ast =N_s\sigma^3/V$
ρs*=Nsσ3/V dependence of the first- and second-order Barker and Henderson perturbation terms for the Lennard-Jones (12-6) fluid at the temperature T* = kT/ε = 1. (a) First perturbation term
$a_1^*=a_1/\epsilon$
a1*=a1/ε
. The circles represent the exact results obtained by Monte Carlo simulation, and the continuous and dashed curves are obtained from Eq. (12) and the approximate relation, Eq. (14), respectively. (b) Second perturbation term
$a_2^*=a_2/\epsilon ^2$
a2*=a2/ε2
. The circles represent the exact results obtained by Monte Carlo simulation, and the continuous, dashed, and dotted curves are obtained from the modified MCA approximation, Eq. (15), the MCA approximation, and the formulation of Gil-Villegas et al.56 using the LCA approximation,85 respectively.

FIG. 1.

The density

$\rho_s^\ast =N_s\sigma^3/V$
ρs*=Nsσ3/V dependence of the first- and second-order Barker and Henderson perturbation terms for the Lennard-Jones (12-6) fluid at the temperature T* = kT/ε = 1. (a) First perturbation term
$a_1^*=a_1/\epsilon$
a1*=a1/ε
. The circles represent the exact results obtained by Monte Carlo simulation, and the continuous and dashed curves are obtained from Eq. (12) and the approximate relation, Eq. (14), respectively. (b) Second perturbation term
$a_2^*=a_2/\epsilon ^2$
a2*=a2/ε2
. The circles represent the exact results obtained by Monte Carlo simulation, and the continuous, dashed, and dotted curves are obtained from the modified MCA approximation, Eq. (15), the MCA approximation, and the formulation of Gil-Villegas et al.56 using the LCA approximation,85 respectively.

Close modal

An exact evaluation of the second-order term a2 (fluctuation contribution) is much more difficult as it requires a knowledge of two-, three-, and four-body correlation functions of the reference system.92 When the superposition approximation is used to evaluate the three- and four-body functions, the integrals may be evaluated numerically. Such calculations were performed by Smith et al.92 showing significant improvement compared to the local compressibility approximation (LCA)85 and macroscopic compressibility approximation (MCA)85 that are commonly used to represent a2. Despite the good description obtained with this more-rigorous approach compared to the simpler LCA and MCA recipes, Smith et al.92 showed that the agreement with molecular simulation remains only fair at higher densities. Moreover, the integrals involved in the calculation of a2 using the superposition approximation cannot be calculated analytically. With the aim of deriving an algebraic equation of state, our objective is therefore to find a simpler, but accurate, expression for the fluctuation term a2. The approach that we advocate here is based on the work of Zhang93 who proposed an improved MCA expression by assuming that the numbers of molecules in neighbouring coordination shells are correlated. Within the improved MCA approach, the term a2 can be expressed as93 

\begin{equation}a_{2}=-\pi \rho _s K^{HS}(1+\chi )\int _{\sigma }^{\infty }g^{HS}_d(r) \big(u_{1}(r)\big)^{2}r^{2} {\rm d}r,\end{equation}
a2=πρsKHS(1+χ)σgdHS(r)u1(r)2r2dr,
(15)

where χ = 8.23η2 is a simple correction pre-factor, and KHS = kT(∂ρs/∂PHS)T is the isothermal compressibility of the hard-sphere reference fluid, which is obtained from the density derivative of the Carnahan and Starling88 expression for the compressibility factor Z = PHS/(ρskT):

\begin{equation}K^{HS}=\frac{(1-\eta )^{4}}{1+4\eta +4\eta ^{2}-4\eta ^{3}+\eta ^{4}}.\end{equation}
KHS=(1η)41+4η+4η24η3+η4.
(16)

By introducing this simple correction, Zhang93 showed that a significant improvement in the representation of a2 is achieved in the case of the square-well potential of range λ = 1.5, especially for systems at intermediate to high densities. This correction fails however in the case of soft-core potentials. A recent modification was proposed by Paricaud87 to address this deficiency by including an additional correction contribution proportional to η5.

A similar approach is used in our work, but employing a more-generic correction factor:

\begin{equation}\chi =f_1(\alpha )\eta x_0^3+f_2(\alpha )\big (\eta x_0^3\big )^5+f_3(\alpha )\big (\eta x_0^3\big )^8 ,\end{equation}
χ=f1(α)ηx03+f2(α)ηx035+f3(α)ηx038,
(17)

where x0 = σ/d. The functions fi (i = 1, 2, 3) depend on both the attractive and repulsive exponents of the Mie potential through a single dimensionless van der Waals-like attractive constant α defined as

\begin{equation}\alpha = \frac{1}{\epsilon \sigma ^3}\int _{\sigma }^{\infty }u^{Mie}(r)r^{2} {\rm d}r=\mathcal {C} \left(\frac{1}{\lambda _a-3}-\frac{1}{\lambda _r-3} \right) .\end{equation}
α=1εσ3σuMie(r)r2dr=C1λa31λr3.
(18)

The specific functional dependence of fi on α is given in Eq. (20).

It is important to emphasize that while the fluctuation term displays a complex non-monotonic behaviour with density, the third-order term is known to exhibit a relatively simple dependence, at least for the square-well potential. This was shown by Espindola-Heredia et al.94 who evaluated the perturbation terms up to fourth order by means of Monte Carlo simulation for fluids characterized by square-well potentials with ranges in the interval 1.1 ⩽ λ ⩽ 3; they then proposed an accurate equation of state using a set of empirical functions that simultaneously reproduced the exact values of the perturbation expansion up to fourth order and the virial coefficients up to third order. In view of the success of their approach, a similar strategy is chosen here for the Mie fluid with the use of the following empirical expression:

\begin{equation}a_{3}=-\epsilon ^3f_4(\alpha )\eta x_0^3\, {\rm exp}\big(f_5(\alpha ) \eta x_0^3+f_6(\alpha )\eta ^2 x_0^6\big) ,\end{equation}
a3=ε3f4(α)ηx03 exp f5(α)ηx03+f6(α)η2x06,
(19)

where the expression is cast in a functional form so that it is independent of temperature. This approximation is chosen in order to reduce the number of adjustable empirical coefficients and is a reasonable assumption as the third-order term is found to be essentially insensitive to the temperature.

There is no unique way of determining the functions fi. A rigorous approach would be to first determine a2 and a3 as functions of density and temperature for different combinations of λr and λa by means of molecular simulation. The functions fi could then be determined by correlating the simulation data as in the work of Espindola-Heredia et al.94 The drawback of this methodology is the relatively large uncertainty in calculating the third-order term a3. Moreover, very small differences in the behaviour of a3 with density are found to lead to significant differences in the location of the critical point. A different methodology is therefore employed using a combination of a knowledge of the accurate simulated values of the fluctuation term a2 and selected simulated VLE and critical-point data for different Mie potentials. The main difficulty of such an approach lies in finding a closed functional form for fi which can be used to represent the selected data using the unique constant α. We have found a satisfactory description with the following functional form:

\begin{equation}f_{i}(\alpha )=\left(\sum _{n=0}^{n=3} \phi _{i,n} \alpha ^n \right)\bigg / \left(1+\sum _{n=4}^{n=6} \phi _{i,n} \alpha ^{n-3}\right) \ i=1,\ldots,6.\end{equation}
fi(α)=n=0n=3ϕi,nαn/1+n=4n=6ϕi,nαn3i=1,...,6.
(20)

Two steps are made in order to determine the coefficients ϕi, n. First, we consider only the coefficients ϕi,n  (i = 1, …, 3) that are involved in the representation of a2, obtained by simultaneously correlating the exact values of a2 and the VLE data, as determined by molecular simulation. For the exact values of a2, we have performed traditional Metropolis Monte Carlo simulations for reference Mie (λra) fluids with selected values of the repulsive exponent, keeping the attractive exponent fixed at its London value of λa = 6 (8-6, 12-6, 14-6, 20-6, 30-6) for a fixed reduced temperature of T* = kT/ε = 1; the simulation procedure is described in the review by Barker and Henderson.85 The literature values of the VLE data95,96 for the saturated liquid densities and vapour pressures for Mie fluids characterized with exponents (8-6), (10-6), (12-6), (15-6), (20-6), (32-6), (36-6), (14-7), (18-9), and (24-12) are considered. Once the coefficients ϕi, n  (i = 1, …, 3) are determined, the remaining terms ϕi, n (i = 4, …, 6) of Eq. (20), corresponding to the third-order perturbation contribution, are estimated from the VLE data, including critical temperatures and critical densities. We should stress that the resulting a3 term therefore corresponds to an effective perturbation contribution which accounts for higher-order terms that have been neglected in the expansion. The resulting values of the coefficients ϕi, n are summarized in Table II.

Table II.

Coefficients ϕi, n for Eq. (20).

nϕ1,nϕ2,nϕ3,nϕ4,nϕ5,nϕ6,nϕ7,n
7.5365557 −359.44 1550.9 −1.19932 −1911.28 9236.9 10 
−37.60463 1825.6 −5070.1 9.063632 21390.175 −129430 10 
71.745953 −3168.0 6534.6 −17.9482 −51320.7 357230 0.57 
−46.83552 1884.2 −3288.7 11.34027 37064.54 −315530 −6.7 
−2.467982 −0.82376 −2.7171 20.52142 1103.742 1390.2 −8 
−0.50272 −3.1935 2.0883 −56.6377 −3264.61 −4518.2 … 
8.0956883 3.7090 40.53683 2556.181 4241.6 … 
nϕ1,nϕ2,nϕ3,nϕ4,nϕ5,nϕ6,nϕ7,n
7.5365557 −359.44 1550.9 −1.19932 −1911.28 9236.9 10 
−37.60463 1825.6 −5070.1 9.063632 21390.175 −129430 10 
71.745953 −3168.0 6534.6 −17.9482 −51320.7 357230 0.57 
−46.83552 1884.2 −3288.7 11.34027 37064.54 −315530 −6.7 
−2.467982 −0.82376 −2.7171 20.52142 1103.742 1390.2 −8 
−0.50272 −3.1935 2.0883 −56.6377 −3264.61 −4518.2 … 
8.0956883 3.7090 40.53683 2556.181 4241.6 … 

The reliability of the semi-empirical expressions for a2 can be assessed by comparing the fluctuation term proposed in Eq. (15), which makes use of the correction function χ defined in Eq. (17), to the traditional MCA approach for the Lennard-Jones fluid (cf. Fig. 1(b)). As can be seen, the corrected MCA expression leads to a significant improvement of the description of a2, especially at the higher densities where the correct trend is captured. Note that the traditional MCA is appropriate at low and intermediate densities but unreliable at high densities. In our comparison we also include calculations with the previous form of a2 employed in the SAFT-VR Mie 2006 EOS,72 which was based on the approximate perturbation potential described by Eq. (13) and the LCA approach.

As can be seen from Fig. 1(b) the LCA is too crude an approximation to accurately represent the second-order perturbation terms of Mie fluids. In Sec. V A 1 we will show that despite this the LCA expression together with the approximate formulation for a1 given in Eq. (14) provides a surprisingly good description of the VLE in the specific case of the Lennard-Jones fluid. This is due to the cancellation of errors between the first- and second-order perturbation contributions.

The behaviour of the perturbation terms of the Helmholtz free energy of Mie fluids is further examined in Fig. 2 where the a1, a2, and a3 terms calculated with our new methodology are compared with exact simulation data for different values of the repulsive exponent (λr = 8, 12, 14, 20, and 30) keeping the attractive exponent fixed at λa = 6. Metropolis Monte Carlo simulations of the corresponding BH reference fluids are carried out here to provide exact values for the three perturbation contributions (see Ref. 85 for details of the methodology). Excellent agreement with the simulation data is apparent for the first-order term for systems characterized by this broad range of repulsive exponents. Small systematic discrepancies can, however, be observed at high densities due to the slight inaccuracies resulting from the calculation of the RDF of the reference hard-sphere fluid with the RHNC integral equation. It is important to recall that the representation of the first-order mean-attractive term a1 follows the rigorous application of the BH theory and is not optimized to simulation data. This negligible difference between the theory and the simulation data for a1 does, however, have some repercussions on the evaluation of the fluctuation term. The second-order contribution a2 is simultaneously correlated to exact values of the fluctuation term and to VLE data. Our methodology therefore results in some discrepancies between the computed values of a2 and the simulation data at high densities (particularly in the case of the more-repulsive potentials) due to a compensation of errors with the first-order term in this region. Nevertheless, the agreement at low densities is very satisfactory, and the overall behaviour for a2 is qualitatively correct for systems with varying repulsive exponents. We should also note that the qualitative trend observed for the fluctuation term of the Mie fluids is very similar to the one reported by Espindola-Heredia et al.94 for the square-well fluids. The novelty of our approach lies in the evaluation of the third-order term using VLE data, including the critical region, for different Mie fluids. One could argue that this constitutes a mere correlation of empirical functions to the coexistence curves of the various systems, but it is reassuring to observe that the resulting density dependence of a3 (cf. Fig. 2(c)) follows the same general behaviour as the simulation data for the third-order term obtained by Espindola-Heredia et al.94 for square-well fluids of variable range. Indeed, it appears that the magnitude of the a3 term is largest near the critical density. The trends observed in Fig. 2(c) are consistent with the slow convergence of the high-temperature expansion in the near-critical region. Interestingly, the magnitude of a3 decreases for increasing values of the repulsive exponent (corresponding to a decrease in the range of the interaction) which is somewhat counterintuitive. This behaviour is again in agreement with the findings of Espindola-Heredia et al.94 We will show in Secs. V A 1 and V A 2 that the inclusion of a third-order term in the high-temperature perturbation expansion allows for a very good description of the critical region without using a specific cross-over treatment.

FIG. 2.

The density

$\rho _s^\ast=N_s \sigma^3/V$
ρs*=Nsσ3/V dependence of the (a) first
$a_1^*=a_1/\epsilon$
a1*=a1/ε
, (b) second
$a_2^*=a_2/\epsilon ^2$
a2*=a2/ε2
, and (c) third
$a_3^*=a_3/\epsilon ^3$
a3*=a3/ε3
Barker and Henderson perturbation terms for selected Mie (8-6), (12-6), (14-6), (20-6), and (30-6) potentials at a temperature T* = kT/ε = 1. The continuous curves represent the calculation using our theoretical expressions, Eqs. (12), (15), and (19). The circles represent the exact values we have determined by Monte Carlo simulation.

FIG. 2.

The density

$\rho _s^\ast=N_s \sigma^3/V$
ρs*=Nsσ3/V dependence of the (a) first
$a_1^*=a_1/\epsilon$
a1*=a1/ε
, (b) second
$a_2^*=a_2/\epsilon ^2$
a2*=a2/ε2
, and (c) third
$a_3^*=a_3/\epsilon ^3$
a3*=a3/ε3
Barker and Henderson perturbation terms for selected Mie (8-6), (12-6), (14-6), (20-6), and (30-6) potentials at a temperature T* = kT/ε = 1. The continuous curves represent the calculation using our theoretical expressions, Eqs. (12), (15), and (19). The circles represent the exact values we have determined by Monte Carlo simulation.

Close modal

2. Algebraic expressions for the first- and second-order perturbation terms

In order to avoid using numerical quadrature techniques to evaluate the integrals in Eqs. (12) and (15) with the reference RDF calculated with the RHNC integral equation, we propose a formal mapping of the integrals by making use of the mean-value theorem. The procedure is at the heart of the SAFT-VR methodology for general spherically symmetrical potentials with hard-repulsive cores.56 The approach is somewhat more complicated to implement when a soft-core reference fluid is considered as will now be explained in detail.

The first-order perturbation term a1 in the BH approach for a Mie potential is obtained by re-expressing Eq. (12) as

\begin{eqnarray}a_{1} &=& 2\pi \rho _s\int _{\sigma }^{\infty }g^{HS}_d(r)u^{Mie}(r)r^{2} {\rm d}r \nonumber \\[4pt]&=& 2\pi \rho _s\int _{d}^{\infty }g^{HS}_d(r)u^{Mie}(r)r^{2} {\rm d}r \nonumber \\[4pt]&& - 2\pi \rho _s\int _{d}^{\sigma }g^{HS}_d(r)u^{Mie}(r)r^{2} {\rm d}r \nonumber \\[4pt]&=& I_{1A} + I_{1B} .\end{eqnarray}
a1=2πρsσgdHS(r)uMie(r)r2dr=2πρsdgdHS(r)uMie(r)r2dr2πρsdσgdHS(r)uMie(r)r2dr=I1A+I1B.
(21)

The computation of the first integral I1A is facilitated by using the algebraic expressions of the first-order term of the hard-core Sutherland potential as follows:

\begin{eqnarray}I_{1A}&=&2\pi \rho _s\int _{d}^{\infty }g^{HS}_d(r)\mathcal {C}\epsilon \left(\left(\frac{\sigma }{r} \right)^{\lambda _{r}}-\left(\frac{\sigma }{r} \right)^{\lambda _{a}} \right)r^{2} {\rm d}r \nonumber \\[4pt]&=&2\pi \rho _s\epsilon \mathcal {C}\bigg [x_0^{\lambda _{r}} d^{3} \int _{1}^{\infty } \frac{1}{x^{\lambda _{r}}}g^{HS}_d(xd)x^{2} {\rm d}x \nonumber \\[4pt]&& -\, x_0^{\lambda _{a}} d^{3} \int _{1}^{\infty } \frac{1}{x^{\lambda _{a}}}g^{HS}_d(xd)x^{2} {\rm d}x \bigg ] ,\end{eqnarray}
I1A=2πρsdgdHS(r)Cεσrλrσrλar2dr=2πρsεC[x0λrd311xλrgdHS(xd)x2dxx0λad311xλagdHS(xd)x2dx],
(22)

where x = r/d is the reduced centre-centre distance between two hard spheres of diameter d.

The expression

\begin{equation}a_1^S=2\pi \rho _s\epsilon d^{3} \int _{1}^{\infty } \left(-\frac{1}{x^{\lambda }} \right) g^{HS}_d(xd) x^{2} {\rm d}x\end{equation}
a1S=2πρsεd311xλgdHS(xd)x2dx
(23)

represents the first-order term

$a_{1}^{S}(\eta ;\lambda )$
a1S(η;λ) of the Helmholtz free energy of a system of hard-core Sutherland particles of diameter d characterized by the range parameter λ. This contribution can be evaluated within the SAFT-VR MVT mapping,56 which allows one to express
$a_{1}^{S}(\eta ;\lambda )$
a1S(η;λ)
as a simple analytical function without loss of accuracy compared to using numerical quadrature. Details of the procedure are discussed in Sec. III A 3. With the availability of an algebraic form for the Sutherland first-order term
$a_{1}^{S}(\eta ;\lambda )$
a1S(η;λ)
, the integral I1A can be evaluated analytically as

\begin{equation}I_{1A}=\mathcal {C}\big(x_0^{\lambda _{a}}a_{1}^{S}(\eta ;\lambda _{a}) - x_0^{\lambda _{r}}a_{1}^{S}(\eta ;\lambda _{r}) \big).\end{equation}
I1A=Cx0λaa1S(η;λa)x0λra1S(η;λr).
(24)

In order to obtain an analytical expression for the second integral in Eq. (21), we approximate

$g^{HS}_d(r)$
gdHS(r) with a linear relation obtained from a first-order Taylor expansion of the contact value of
$g^{HS}_d(r)$
gdHS(r)
and its first derivative with respect to the separation r. Nezbeda and Iglesias-Silva97 have pointed out that such a description is expected to be reasonable for short ranges of integration. The hard-sphere RDF
$g^{HS}_d(r)$
gdHS(r)
can therefore be evaluated with the following expression:

\begin{equation}g^{HS}_d(xd) \approx g^{HS}_d(d)+(x-1)\left(\frac{{\rm d} g^{HS}_d(xd)}{{\rm d}x}\right)_{x=1} .\end{equation}
gdHS(xd)gdHS(d)+(x1)dgdHS(xd)dxx=1.
(25)

The second integral I1B can then be written as

\begin{eqnarray}I_{1B} & \approx & -2\pi \rho _s d^3\Bigg (g^{HS}_d(d)\int _{1}^{x_0}u^{Mie}(xd)x^{2} {\rm d}x\nonumber\\[4pt]&&+\, \left(\frac{{\rm d} g^{HS}_d(xd)}{{\rm d}x}\right)_{x=1}\int _{1}^{x_0}u^{Mie}(xd)x^{2}(x-1) {\rm d}x\Bigg ), \nonumber\!\!\!\!\! \\\end{eqnarray}
I1B2πρsd3(gdHS(d)1x0uMie(xd)x2dx+dgdHS(xd)dxx=11x0uMie(xd)x2(x1)dx),
(26)

with

\begin{equation}u^{Mie}(xd)=\mathcal {C}\epsilon \left(\left( \frac{x_0}{x}\right)^{\lambda _{r}} - \left(\frac{x_0}{x} \right)^{\lambda _{a}}\right) .\end{equation}
uMie(xd)=Cεx0xλrx0xλa.
(27)

To simplify Eq. (26) we define

\begin{equation}I_{\lambda }(\lambda )=\int _{1}^{x_0}\frac{x^2}{x^{\lambda }} {\rm d}x=-\frac{\left(x_0 \right)^{3-\lambda }-1}{\lambda -3}\end{equation}
Iλ(λ)=1x0x2xλdx=x03λ1λ3
(28)

and

\begin{eqnarray}J_{\lambda }(\lambda )&=&\int _{1}^{x_0}\frac{(x^3-x^2)}{x^{\lambda }} {\rm d}x\nonumber \\[4pt]&=&-\frac{\left(x_0 \right)^{4-\lambda }(\lambda -3)-\left(x_0 \right)^{3-\lambda }(\lambda -4)-1}{(\lambda -3)(\lambda -4)} . \qquad\end{eqnarray}
Jλ(λ)=1x0(x3x2)xλdx=x04λ(λ3)x03λ(λ4)1(λ3)(λ4).
(29)

By then using the Carnahan and Starling expression for

$g^{HS}_d(d)$
gdHS(d) (from
$Z^{HS}=1+ 4 \eta g^{HS}_d(d)$
ZHS=1+4ηgdHS(d)
) and the Percus-Yevick approximation for the derivative at contact,

\begin{equation}g^{HS}_d(d) =\frac{1-\eta /2}{(1-\eta )^3},\end{equation}
gdHS(d)=1η/2(1η)3,
(30)
\begin{equation}\left(\frac{{\rm d} g^{HS}_d(xd)}{{\rm d}x}\right)_{x=1} =-\frac{9\eta (1+\eta )}{2(1-\eta )^3} ,\end{equation}
dgdHS(xd)dxx=1=9η(1+η)2(1η)3,
(31)

a simple analytical expression for I1B can be obtained as

\begin{equation}I_{1B}=\mathcal {C}\big(x_0^{\lambda _{a}}B(\eta ;\lambda _{a})-x_0^{\lambda _{r}}B(\eta ;\lambda _{r}) \big) ,\end{equation}
I1B=Cx0λaB(η;λa)x0λrB(η;λr),
(32)

where

\begin{equation}B(\eta ;\lambda )=12\eta \epsilon \left(\frac{1-\eta /2}{(1-\eta )^3} I_{\lambda }(\lambda )-\frac{9\eta (1+\eta )}{2(1-\eta )^3}J_{\lambda }(\lambda ) \right).\end{equation}
B(η;λ)=12ηε1η/2(1η)3Iλ(λ)9η(1+η)2(1η)3Jλ(λ).
(33)

As a result the first-order term a1 for the Mie fluid can be approximated with the following algebraic expression:

\begin{eqnarray}a_{1}&=&\mathcal {C} \big [ x_0^{\lambda _{a}} \big ( a_{1}^{S}(\eta ;\lambda _{a}) + B(\eta ;\lambda _{a}) \big ) \nonumber \\[4pt]&& -x_0^{\lambda _{r}} \big ( a_{1}^{S}(\eta ;\lambda _{r}) +B(\eta ;\lambda _{r} )\big )\big ] .\end{eqnarray}
a1=C[x0λaa1S(η;λa)+B(η;λa)x0λra1S(η;λr)+B(η;λr)].
(34)

The same procedure is used to develop an algebraic expression for the second-order term a2:

\begin{eqnarray}a_{2}&=&-\pi \rho _s K^{HS}(1+\chi )\int _{\sigma }^{\infty }g^{HS}_d(r) \big(u^{Mie}(r)\big)^{2}r^{2} {\rm d}r \nonumber \\[6pt]&=&-\, \pi \rho _s K^{HS}(1+\chi )\int _{d}^{\infty }g^{HS}_d(r)\big(u^{Mie}(r)\big)^{2}r^{2} {\rm d}r \nonumber \\[6pt]&& +\, \pi \rho _s K^{HS}(1+\chi )\int _{d}^{\sigma }g^{HS}_d(r) \big(u^{Mie}(r)\big)^{2}r^{2} {\rm d}r \nonumber \\[6pt]&=&I_{2A}+I_{2B} .\end{eqnarray}
a2=πρsKHS(1+χ)σgdHS(r)uMie(r)2r2dr=πρsKHS(1+χ)dgdHS(r)uMie(r)2r2dr+πρsKHS(1+χ)dσgdHS(r)uMie(r)2r2dr=I2A+I2B.
(35)

The integrals I2A and I2B can be written in terms of the first-order perturbation terms

$a_1^S$
a1S of the corresponding hard-core Sutherland potentials. It can be easily shown that

\begin{eqnarray}a_{2}&=&\frac{1}{2}K^{HS}(1+\chi )\epsilon \ \mathcal {C}^{2} \bigg [ x_0^{2\lambda _{a}} \Big ( a_{1}^{S}(\eta ;2\lambda _{a})+ B(\eta ;2\lambda _{a})\Big ) \nonumber \\[6pt]&&-2x_0^{\lambda _{a}+\lambda _{r}} \Big ( a_{1}^{S}(\eta ;\lambda _{a}+\lambda _{r})+ B(\eta ;\lambda _{a}+\lambda _{r})\Big ) \nonumber \\[6pt]&&+\, x_0^{2\lambda _{r}} \Big ( a_{1}^{S}(\eta ;2\lambda _{r})+ B(\eta ;2\lambda _{r})\Big ) \bigg ] .\end{eqnarray}
a2=12KHS(1+χ)εC2[x02λaa1S(η;2λa)+B(η;2λa)2x0λa+λra1S(η;λa+λr)+B(η;λa+λr)+x02λra1S(η;2λr)+B(η;2λr)].
(36)

Finally, we recall that the expression for the third-order perturbation contribution a3 proposed in Eq. (19) does not require a knowledge of the RDF of the reference fluid. Its evaluation is straightforward using Eqs. (19) and (20).

3. Improved mean-value theorem mapping for the first-order perturbation term of the Sutherland potential |$a_{1} ^{s}$|a1s

As will have become apparent from Sec. III A 2, the development of an analytical expression for the Helmholtz free energy of a fluid of Mie spherical segments within the SAFT-VR approach56 involves the evaluation of the first-order terms

$a_{1} ^{S}(\eta ;\lambda )$
a1S(η;λ) of systems of hard-core Sutherland potentials of range λ (cf. Eq. (23)), where

\begin{equation}a_{1} ^{S}(\eta ;\lambda ) = 12\epsilon \eta \int _{1}^{\infty } \left(-\frac{1}{x^{\lambda }} \right) g^{HS}_d(xd) x^{2} {\rm d}x.\end{equation}
a1S(η;λ)=12εη11xλgdHS(xd)x2dx.
(37)

In the SAFT-VR treatment,56 the RDF for the hard-sphere reference

$g^{HS}_d(xd)$
gdHS(xd) is calculated numerically by making use of an integral-equation theory. In our current work the RHNC theory with the accurate hard-sphere bridge function of Malijevsky and Labik90 is used as it is known to provide a hard-sphere RDF which is in very close agreement with computer-simulation data. The integral in Eq. (37) is then evaluated as a function of the exponent λ and packing fraction η. Once the range and density dependence of
$a_{1} ^{S}(\eta ;\lambda )$
a1S(η;λ)
are known numerically the MVT can be applied to determine an analytical form for the integral by expressing Eq. (37) as

\begin{eqnarray}a_{1} ^{S}(\eta ;\lambda ) &=& 12\epsilon \eta g^{HS}_d(\xi )\int _{1}^{\infty } \left(-\frac{1}{x^{\lambda }} \right) x^{2} {\rm d}x \nonumber \\[12pt]\nonumber\\&=& -12\epsilon \eta g^{HS}_d(\xi )\left(\frac{1}{\lambda -3} \right) ,\end{eqnarray}
a1S(η;λ)=12εηgdHS(ξ)11xλx2dx=12εηgdHS(ξ)1λ3,
(38)

where ξ is a distance satisfying the MVT equality. As shown in the original SAFT-VR treatment56 it is possible to map

$g^{HS}_d(\xi )$
gdHS(ξ) to the contact RDF
$g^{HS}_d(d)$
gdHS(d)
evaluated at an effective density ηeff. By using the expression for the contact value of the hard-sphere RDF which is consistent with the Carnahan and Starling88 EOS (cf. Eq. (30)), a compact relation for
$a_{1} ^{S}(\eta ;\lambda )$
a1S(η;λ)
can be obtained as56 

\begin{equation}a_{1} ^{S}(\eta ;\lambda ) = -12\epsilon \eta \left(\frac{1}{\lambda -3} \right)\frac{1-\eta _{\scriptsize\textit{eff}}(\eta ; \lambda )/2}{(1-\eta _{\scriptsize\textit{eff}}(\eta ; \lambda ))^{3}} .\end{equation}
a1S(η;λ)=12εη1λ31ηeff(η;λ)/2(1ηeff(η;λ))3.
(39)

An algebraic expression for the integral in Eq. (37) can be obtained once the dependence of the effective packing fraction ηeff on λ and η is established. The original parameterization56 for ηeff(η; λ) was derived for exponents of the Sutherland potential in the range 3 < λ ⩽ 12. We now propose a new mapping in order to account for a broader range of exponents, particularly for the more-repulsive systems, i.e., 5 < λ ⩽ 100:

\begin{equation}\eta _{\scriptsize\textit{eff}}(\eta ; \lambda )=c_{1}(\lambda )\eta +c_{2}(\lambda )\eta ^{2}+c_{3}(\lambda )\eta ^{3}+c_{4}(\lambda )\eta ^{4} ,\end{equation}
ηeff(η;λ)=c1(λ)η+c2(λ)η2+c3(λ)η3+c4(λ)η4,
(40)

with

\begin{equation}\left( \begin{array}{c}c_{1} \\[8pt]c_{2} \\[8pt]c_{3} \\[8pt]c_{4} \\[8pt]\end{array} \right) = \left( \begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0.81096 & 1.7888 & -37.578 & 92.284 \\[8pt]1.0205 & -19.341 & 151.26 & -463.50 \\[8pt]-1.9057 & 22.845 & -228.14 & 973.92 \\[8pt]1.0885 & -6.1962 & 106.98 & -677.64 \\[8pt]\end{array} \right) . \left( \begin{array}{c}1 \\[8pt]1/\lambda \\[8pt]1/\lambda ^{2} \\[8pt]1/\lambda ^{3} \\[8pt]\end{array} \right) .\end{equation}
c1c2c3c4=0.810961.788837.57892.2841.020519.341151.26463.501.905722.845228.14973.921.08856.1962106.98677.64.11/λ1/λ21/λ3.
(41)

The resulting values of the first-order Sutherland term

$a_1^s(\eta ;\lambda )$
a1s(η;λ) obtained with Eq. (39) along with the numerical results calculated from Eq. (37) with the hard-sphere RDF from the RHNC integral-equation theory of Malijesvsky and Labik90 are shown in Fig. 3. It can be seen that the algebraic expression provides an excellent representation of
$a_1^S(\eta ,\lambda )$
a1S(η,λ)
over the entire range of λ for which the effective packing fraction has been parameterized.

FIG. 3.

The density

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V dependence of the first
$a_1^{S*}=a_1^S/\epsilon$
a1S*=a1S/ε
Barker and Henderson perturbation term of the Sutherland fluid with varying attractive exponent λ. The circles represent the essentially exact values of numerical quadrature with the radial distribution function of the hard-sphere reference system obtained using the RHNC integral equation. The continuous curves correspond to the analytical correlation determined by employing the mean-value theorem, cf. Eqs. (39)–(41).

FIG. 3.

The density

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V dependence of the first
$a_1^{S*}=a_1^S/\epsilon$
a1S*=a1S/ε
Barker and Henderson perturbation term of the Sutherland fluid with varying attractive exponent λ. The circles represent the essentially exact values of numerical quadrature with the radial distribution function of the hard-sphere reference system obtained using the RHNC integral equation. The continuous curves correspond to the analytical correlation determined by employing the mean-value theorem, cf. Eqs. (39)–(41).

Close modal

The calculation of the RDF of Mie fluids of variable ranges evaluated at contact is revisited and derived through a rigorous application of the perturbation theory of Barker and Henderson.85 It is important to recall that an accurate representation of the contact value of the RDF of the reference monomer fluid is a prerequisite in a Wertheim TPT1 description of the thermodynamic properties of chain fluids formed from Mie segments.

1. New SAFT-VR Mie formulation for gMie(σ)

Just as for the Helmholtz free energy, a high-temperature perturbation expansion can be used to represent the RDF for Mie fluids:

\begin{equation}g^{Mie}(r)=g_{0}(r)+\beta \epsilon g_{1}(r)+\left(\beta \epsilon \right)^{2}g_{2}(r)+ ..... .\end{equation}
gMie(r)=g0(r)+βεg1(r)+βε2g2(r)+......
(42)

Another possibility is to expand the logarithm of the RDF, ln gMie(r). When the expansion is taken to second order the RDF can be expressed as85 

\begin{eqnarray}g^{Mie}(r)&\approx & g_{0}(r) {\rm exp} \big( \beta \epsilon g_{1}(r)/g_{0}(r)+ (\beta \epsilon )^{2} g_{2}(r)/g_{0}(r)\big). \nonumber\!\!\!\!\! \\\end{eqnarray}
gMie(r)g0(r) exp βεg1(r)/g0(r)+(βε)2g2(r)/g0(r).
(43)

The latter representation is preferred as it allows for a faster convergence of the series for the contact value gMie(σ) than the high-temperature expression of Eq. (42). Another attractive aspect of expanding ln gMie(r) is that the resulting expression for gMie(r) will never give rise to negative values. This turns out to be very important in the Wertheim TPT1 implementation, since the logarithm of the RDF at contact is required for the calculation of the contribution to the Helmholtz free energy due to the formation of chains of Mie segments. One should note that a second-order expansion for the structure of the fluid is consistent with a third-order expansion for the free energy as the RDF essentially corresponds to a variation of the latter with respect to density.

In order to calculate the zeroth-order term g0(r), the approximation of Barker and Henderson is used so that the RDF for the soft-repulsive potential u0 defined in Eq. (4) is approximated as85 

\begin{equation}g_0(r) \approx g^{HS}_d(r) ,\end{equation}
g0(r)gdHS(r),
(44)

where

$g^{HS}_d$
gdHS is the RDF of the reference system of hard spheres of diameter d defined in Eq. (7), which accounts for the softness of u0. For the contact value at r = σ, the zeroth-order term can therefore be approximated as

\begin{equation}g_0(\sigma ) \approx g^{HS}_d(\sigma ) .\end{equation}
g0(σ)gdHS(σ).
(45)

It is important to recognize that even though σ is the contact distance of the Mie fluid, it is not the contact distance of the hard-sphere reference system; the contact distance of the HS fluid is d which is implicitly temperature dependent (cf. Eq. (7)). As σ > d the calculation of the term

$g^{HS}_d$
gdHS on the right-hand side of Eq. (45) therefore requires a knowledge of the RDF of the HS fluid at reduced distances x0 = σ/d > 1. Various semi-empirical representations are available for
$g^{HS}_d(r)$
gdHS(r)
. In our current work we use the expression proposed by Boublík98 which is valid for
$1<x_0<\root \of {2}$
1<x0<2
:

\begin{eqnarray}g^{HS}_d(\sigma )=g^{HS}_d(x_{0}d)={\rm exp} \big(k_{0}+k_{1}x_{0}+k_{2}x_{0}^{2}+k_{3}x_0^{3}\big) , \nonumber\!\!\!\!\!\!\! \\\end{eqnarray}
gdHS(σ)=gdHS(x0d)= exp k0+k1x0+k2x02+k3x03,
(46)

where the density-dependent coefficients ki are given by98 

\begin{equation}k_{0}=-\ln (1-\eta )+\frac{42\eta -39\eta ^{2}+9\eta ^{3}-2\eta ^{4}}{6(1-\eta )^{3}} ,\end{equation}
k0=ln(1η)+42η39η2+9η32η46(1η)3,
(47)
\begin{equation}k_{1}=\frac{(\eta ^{4}+6\eta ^{2}-12\eta )}{2(1-\eta )^{3}} ,\end{equation}
k1=(η4+6η212η)2(1η)3,
(48)
\begin{equation}k_{2}=\frac{-3\eta ^{2}}{8(1-\eta )^{2}} ,\end{equation}
k2=3η28(1η)2,
(49)

and

\begin{equation}k_{3}=\frac{(-\eta ^{4}+3\eta ^{2}+3\eta )}{6(1-\eta )^{3}} .\end{equation}
k3=(η4+3η2+3η)6(1η)3.
(50)

This representation for g0(σ), cf. Eqs. (45) and (46), differs from the approximation made in previous work65,72 where the zeroth-order RDF was simply evaluated with the Carnahan and Starling88 expression for the contact value at the packing fraction η, essentially corresponding to the approximation that d = σ. The values of our new zeroth-order term are always smaller than those of the previous formulation but both tend to the correct zero-density limit,

${\rm lim}_{\rho \rightarrow 0}\ g^{HS}_d(\sigma )=1$
lim ρ0gdHS(σ)=1⁠. In order to determine the contact value of the first-order term g1(σ) we follow the usual SAFT-VR recipe56 employing a self-consistent representation for the pressure P from the Clausius virial theorem85 and from the density derivative of the Helmholtz free energy, P = −(∂A/∂V)N, T = (∂A/∂ρs)N, Ts/V). The two expressions must provide the same values for the compressibly factor Z = P/(ρskT) of the Mie fluid:

\begin{equation}Z= 1 - \frac{4\eta }{d^{3}} \beta \int _{0}^{\infty } r \frac{{\rm d} u^{Mie}(r)}{{\rm d} r} g^{Mie}(r)r^{2} {\rm d}r \\end{equation}
Z=14ηd3β0rduMie(r)drgMie(r)r2dr
(51)

for the virial relation, and

\begin{equation}Z=Z^{HS}+\beta \eta \frac{\partial a_{1}}{\partial \eta }+(\beta )^{2}\eta \frac{\partial a_{2}}{\partial \eta } \\end{equation}
Z=ZHS+βηa1η+(β)2ηa2η
(52)

for the thermodynamic derivative, where ZHS is the compressibility factor of the HS system of diameter d. When one approximates the reference potential u0(r) ≈ uHS(r) as that of hard spheres of diameter d, noticing that the exponential function,

\begin{equation}{\rm exp}\left(\frac{-u^{HS}(r)}{kT} \right)= \left\lbrace \begin{array}{c@{\quad}c}0 & r \le d \\[4pt]1 & r > d \\[4pt]\end{array} \right. ,\end{equation}
exp uHS(r)kT=0rd1r>d,
(53)

is a step function so that its derivative with respect to r is a δ function at r = d, it then follows that Eq. (51) can be rewritten as

\begin{equation}Z=1 + 4\eta g^{Mie}(d)-\frac{4\eta }{d^{3}} \beta \int _{\sigma }^{\infty } \frac{{\rm d} u_{1}(r)}{{\rm d} r} g^{Mie}(r)r^{3} {\rm d}r.\end{equation}
Z=1+4ηgMie(d)4ηd3βσdu1(r)drgMie(r)r3dr.
(54)

After replacing gMie(r) by Eqs. (43) and (44), and expanding the exponential function one can write

\begin{eqnarray}Z&=&Z^{HS} + 4\eta \beta \epsilon g_{1}(d)+ 4\eta (\beta \epsilon )^2 g_{2}(d)\nonumber \\[4pt]&&-\frac{4\eta }{d^{3}}\beta \int _{\sigma }^{\infty } \frac{{\rm d} u_{1}(r)}{{\rm d} r} g^{HS}_d(r)r^{3} {\rm d}r \nonumber \\[4pt]&& - \frac{4\eta }{d^{3}}\beta ^2 \epsilon \int _{\sigma }^{\infty } \frac{{\rm d} u_{1}(r)}{{\rm d} r} g_{1}(r)r^{3} {\rm d}r+ \cdots \nonumber \\[4pt]&\approx &Z^{HS} + 4\eta \beta \epsilon g_{1}(d)+ 4\eta (\beta \epsilon )^2 g_{2}(d)\nonumber \\[4pt]&&+\, \frac{4\eta }{d^{3}}\beta \epsilon \mathcal {C}\lambda _{r}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{\lambda _{r}}g^{HS}_d(r) {\rm d}r \nonumber \\[4pt]&& - \frac{4\eta }{d^{3}}\beta \epsilon \mathcal {C}\lambda _{a}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{\lambda _{a}}g^{HS}_d(r) {\rm d}r\nonumber \\[4pt]&&+\, \frac{4\eta }{d^{3}}(\beta \epsilon )^2 \mathcal {C}\lambda _{r}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{\lambda _{r}}g_{1}(r) {\rm d}r \nonumber \\[4pt]&& - \frac{4\eta }{d^{3}}(\beta \epsilon )^2 \mathcal {C}\lambda _{a}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{\lambda _{a}}g_{1}(r) {\rm d}r.\end{eqnarray}
Z=ZHS+4ηβεg1(d)+4η(βε)2g2(d)4ηd3βσdu1(r)drgdHS(r)r3dr4ηd3β2εσdu1(r)drg1(r)r3dr+ZHS+4ηβεg1(d)+4η(βε)2g2(d)+4ηd3βεCλrσr2σrλrgdHS(r)dr4ηd3βεCλaσr2σrλagdHS(r)dr+4ηd3(βε)2Cλrσr2σrλrg1(r)dr4ηd3(βε)2Cλaσr2σrλag1(r)dr.
(55)

By comparing Eq. (55) with Eq. (52), compact expressions for the first-order and second-order terms of the RDF at the distance r = d are obtained:

\begin{eqnarray}g_{1}(d)&=&\frac{1}{4\epsilon } \frac{\partial a_{1}}{\partial \eta } - \frac{\mathcal {C}\lambda _{r}}{d^{3}}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{\lambda _{r}}g^{HS}_d(r) {\rm d}r\nonumber \\[4pt]&& +\, \frac{\mathcal {C}\lambda _{a}}{d^{3}}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{\lambda _{a}}g^{HS}_d(r){\rm d}r ,\end{eqnarray}
g1(d)=14εa1ηCλrd3σr2σrλrgdHS(r)dr+Cλad3σr2σrλagdHS(r)dr,
(56)

and

\begin{eqnarray}g_{2}(d)&=&\frac{1}{4\epsilon ^2} \frac{\partial a_{2}}{\partial \eta } - \frac{\mathcal {C}\lambda _{r}}{d^{3}}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{\lambda _{r}}g_{1}(r) {\rm d}r\nonumber \\[4pt]&& +\, \frac{\mathcal {C}\lambda _{a}}{d^{3}}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{\lambda _{a}}g_{1}(r) {\rm d}r.\end{eqnarray}
g2(d)=14ε2a2ηCλrd3σr2σrλrg1(r)dr+Cλad3σr2σrλag1(r)dr.
(57)

At this stage it is useful to make some observations about these expressions. It would appear that by making use of this self-consistent method, the evaluation of the perturbation terms of the RDF at the distance r = σ is not possible. We note, however, that the approximations g1(σ) ≈ g1(d) and g2(σ) ≈ g2(d) yield a good representation of the full gMie(σ) as long as the zeroth-order contribution

$g_0(\sigma ) \approx g^{HS}_d(\sigma )$
g0(σ)gdHS(σ) is calculated rigorously from BH theory at the contact distance r = σ with the Boublík expression.98 The crude approximation that gMie(σ) ≈ gMie(d) does not yield an accurate representation due to the rapid increase in gMie(r) between d and σ. It is important to emphasize that when the approximated perturbation potential as described in Eq. (13) is used in our current development then one recovers the expression for g1(d) used in the previous SAFT-VR Mie 2006 EOS.72 

While the first-order term is easily evaluated with a knowledge of the hard-sphere RDF of the reference HS fluid, the second-order term requires an intricate calculation of g1(r) for the Mie potential over all distances. Nevertheless, Eq. (57) can be simplified further by making use of the expression for the configurational energy together with the MCA approximation. After inserting the high-temperature expansion of the RDF in the energy, a correspondence between the first-order term of the RDF and the second-order contribution to the Helmholtz free energy can be obtained as99 

\begin{eqnarray}a_{2} &=& \frac{6\eta }{d^3}\int _{\sigma }^{\infty }g_1(r)\frac{u^{Mie}(r)}{\epsilon }r^{2} {\rm d}r \nonumber \\[4pt]&=& \frac{6\eta }{d^3}\int _{\sigma }^{\infty }g_1(r)\frac{u_{1}(r)}{\epsilon }r^{2} {\rm d}r .\end{eqnarray}
a2=6ηd3σg1(r)uMie(r)εr2dr=6ηd3σg1(r)u1(r)εr2dr.
(58)

If one then considers the MCA approximation, which relates a2 to the RDF of the HS fluid,

\begin{eqnarray}{a}_{2}^{\scriptsize\textit{MCA}}=-\frac{6\eta }{d^3} K^{HS}\int _{\sigma }^{\infty }g^{HS}_d(r) \big(u_{1}(r)\big)^{2}r^{2} {\rm d}r ,\end{eqnarray}
a2MCA=6ηd3KHSσgdHS(r)u1(r)2r2dr,
(59)

and then equates Eqs. (58) and (59), the corresponding MCA approximation for g1(r) for separations r > σ is obtained as

\begin{eqnarray}g_1(r)\approx {g}_1^{\scriptsize\textit{MCA}}(r)=-K^{HS}g^{HS}_d(r)\frac{u_{1}(r)}{\epsilon } .\end{eqnarray}
g1(r)g1MCA(r)=KHSgdHS(r)u1(r)ε.
(60)

Note that the MCA approximation for g1(r) at the contact distance σ is not useful since it leads to the incorrect value,

$g_1^{\scriptsize\textit{MCA}}(\sigma) = 0$
g1MCA(σ)=0⁠. By substituting Eq. (60) into Eq. (57) we obtain the MCA approximation for the second-order RDF of the Mie potential at the distance d which is now written in terms of the known RDF of the HS fluid:

\begin{eqnarray}{g}_{2}^{\scriptsize\textit{MCA}}(d)&=&\frac{1}{4\epsilon ^2} \frac{\partial a_{2}^{\scriptsize\textit{MCA}}}{\partial \eta } \nonumber \\[4pt]&&+\, \frac{\mathcal {C}^2\lambda _{r} K^{HS}}{d^{3}}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{2\lambda _{r}}g^{HS}_d(r) {\rm d}r \nonumber \\[4pt]&&-\, \frac{\mathcal {C}^2\lambda _{r}K^{HS}}{d^{3}}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{\lambda _{r}+\lambda _{a}}g^{HS}_d(r) {\rm d}r \nonumber \\[4pt]&&-\, \frac{\mathcal {C}^2\lambda _{a}K^{HS}}{d^{3}}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{\lambda _{r}+\lambda _{a}}g^{HS}_{d}(r) {\rm d}r \nonumber \\[4pt]&&+\, \frac{\mathcal {C}^2\lambda _{a}K^{HS}}{d^{3}}\int _{\sigma }^{\infty } r^{2} \left( \frac{\sigma }{r}\right)^{2\lambda _{a}}g^{HS}_d(r) {\rm d}r . \nonumber\!\!\!\!\! \\\end{eqnarray}
g2MCA(d)=14ε2a2MCAη+C2λrKHSd3σr2σr2λrgdHS(r)drC2λrKHSd3σr2σrλr+λagdHS(r)drC2λaKHSd3σr2σrλr+λagdHS(r)dr+C2λaKHSd3σr2σr2λagdHS(r)dr.
(61)

In order to validate the approximations made in our new approach, we combine the zeroth-, first- and second-order contributions given by Eqs. (46), (56), and (61), respectively, into the perturbation expression, Eq. (43), to evaluate the contact value of the RDF of the Mie fluid, assuming that g1(d) ≈ g1(σ) and g2(d) ≈ g2(σ). The resulting values of gLJ(σ) for the LJ (12-6) fluid are shown in Fig. 4(a) as function of density at a reduced temperature T* = 1, making a comparison with the “exact” contact values of the RDF obtained from correlations to simulation data.63 To gain insight into the influence of truncating the series expansion, we also display the resulting values obtained with the zeroth- and first-order perturbation expansions. The introduction of a second-order term in the perturbation is seen to enhance the accuracy of the RDF at low temperatures (cf. Fig. 4(b)), although the MCA approximation for g2 still fails to capture the slope of the contact value of the RDF at zero density for T* = 1 (cf. Fig. 4(a)). The positive slope obtained for the temperature dependence of the contact value of LJ fluids is also seen for square-well fluids, as described theoretically by Tavares et al.,100 who demonstrated that for a SW interaction with a range λ ⩾ 1.5, a positive slope is obtained for any temperature. It has been argued that it is not possible to capture the subtle temperature dependence of the contact value with high-temperature perturbation expansions.100 It is, however, reasonable to assume that, if no approximations (such as the MCA) were made in calculating the second- and higher-order terms, then the inadequacies at low density and temperature would not be present.

FIG. 4.

The density

$\rho^\ast _s=N_s\sigma^3/V$
ρs*=Nsσ3/V and temperature T* = kT/ε dependence of the contact value of the radial distribution function g(σ) of the monomer fluid: (a) Lennard-Jones (12-6) fluid at fixed temperature T* = 1; (b) Lennard-Jones (12-6) fluid at fixed density
$\rho^\ast_s = 0.15$
ρs*=0.15
; and (c) Lennard-Jones (12-6) and Mie (30-6) fluids at a fixed temperature of T* = 1. The curves correspond to calculations using Eq. (43) with different truncations of the perturbation series: zeroth-order (dotted), first-order (dashed), second-order MCA (dashed-dotted), corrected second-order expansion (continuous). The circles represent the essentially exact correlated simulation data of Johnson et al.78 for the Lennard-Jones fluid and our simulation data for the Mie (30-6) fluid.

FIG. 4.

The density

$\rho^\ast _s=N_s\sigma^3/V$
ρs*=Nsσ3/V and temperature T* = kT/ε dependence of the contact value of the radial distribution function g(σ) of the monomer fluid: (a) Lennard-Jones (12-6) fluid at fixed temperature T* = 1; (b) Lennard-Jones (12-6) fluid at fixed density
$\rho^\ast_s = 0.15$
ρs*=0.15
; and (c) Lennard-Jones (12-6) and Mie (30-6) fluids at a fixed temperature of T* = 1. The curves correspond to calculations using Eq. (43) with different truncations of the perturbation series: zeroth-order (dotted), first-order (dashed), second-order MCA (dashed-dotted), corrected second-order expansion (continuous). The circles represent the essentially exact correlated simulation data of Johnson et al.78 for the Lennard-Jones fluid and our simulation data for the Mie (30-6) fluid.

Close modal

Here we propose the following simple empirical correction to the second-order term in order to make the contact value of the RDF of the Mie potential accurate at low temperature:

\begin{eqnarray}g_2(d)= (1+\gamma _c){g}_{2}^{\scriptsize\textit{MCA}}(d) ,\end{eqnarray}
g2(d)=(1+γc)g2MCA(d),
(62)

where γc is a correction factor that is a function of both density and temperature,

\begin{eqnarray}\gamma _c&=&\phi _{7,0} \big( {-}\tanh \big (\phi _{7,1}(\phi _{7,2}-\alpha )\big )+1 \big)\eta x_0^3 \theta \nonumber \\[4pt]&& \times\, {\rm exp}\left(\phi _{7,3}\eta x_0^3 + \phi _{7,4}\eta ^2 x_0^6 \right) ,\end{eqnarray}
γc=ϕ7,0tanhϕ7,1(ϕ7,2α)+1ηx03θ× exp ϕ7,3ηx03+ϕ7,4η2x06,
(63)

where θ = exp(βε) − 1, and α is the reduced constant defined in Eq. (18). The coefficients ϕ7, i  (i = 0, 3, and 4) are obtained by correlating the contact values of the RDF of the LJ fluid, as described by the expression of Johnson et al.,63 at low temperatures and densities. The functional form of the factor is chosen such that both Eqs. (61) and (62) provide an equivalent description at high density and high temperature. We find that a temperature correction, characterized by the parameter θ, is necessary to capture the RDF of LJ fluids accurately at very low temperature (cf. Fig. 4(b)). In this sense the expression for g2 given by Eq. (62) represents an effective term which accounts for the contributions of the higher-order terms in the perturbation expansion. A hyperbolic tangent is used in order to bring the correction γc to zero for very short-ranged potentials. This allows one to capture the decreasing magnitude of the density dependence of the slope of the RDF at zero density for Mie potentials with increasing repulsive exponents. We illustrate this effect in Fig. 4(c) where our corrected second-order RDF is compared with contact values of the RDF of LJ (12-6) and Mie (30-6) fluids obtained by molecular simulation.63 Overall the agreement is satisfactory in both cases, with a slight over-estimation of the contact values of the RDF at higher densities; nevertheless, the theory enables one to capture accurately the slope of the RDF at high density and its dependence on the repulsive part of the potential. To assess further the accuracy of the corrected second-order contact value of the RDF, our new description is compared in Fig. 5 with the empirical correlation to simulation data for the LJ fluid63 as well as with the corresponding approximate expressions used in the development of the SAFT-VR Mie 2006 EOS72 and the SAFT-VR LJC EOS.65 It can be seen that our corrected second-order RDF is in much better agreement with the simulation data than the other approximations, and is now seen to approach the correct limit at high temperatures; this is a result of using of the correct contact distance σ in the zeroth-order term instead of the approximate diameter of the reference HS system d employed previously.72 Overall, the reformulation proposed here is accurate over the entire fluid region including low temperatures. A key aspect of this work is that the new expression presented for the contact value of the RDF is generic and applicable for any value of the repulsive and attractive exponents of the Mie potential (so long as λa > 4), as opposed to the correlation of Johnson et al.63 which is appropriate only for the LJ fluid.

FIG. 5.

(a)–(f) The density

$\rho _s^\ast =N_s\sigma^3/V$
ρs*=Nsσ3/V and temperature T* = kT/ε dependence of the contact value of the radial distribution function gLJ(σ) of the monomer Lennard-Jones (12-6) fluid. The dashed and continuous curves represent the description with the previous72 and our new SAFT-VR Mie approaches, respectively; the dotted curves correspond to the description with SAFT-VR LJC65 approach. The circles represent the essentially exact simulation data of Johnson et al.78 

FIG. 5.

(a)–(f) The density

$\rho _s^\ast =N_s\sigma^3/V$
ρs*=Nsσ3/V and temperature T* = kT/ε dependence of the contact value of the radial distribution function gLJ(σ) of the monomer Lennard-Jones (12-6) fluid. The dashed and continuous curves represent the description with the previous72 and our new SAFT-VR Mie approaches, respectively; the dotted curves correspond to the description with SAFT-VR LJC65 approach. The circles represent the essentially exact simulation data of Johnson et al.78 

Close modal

2. Algebraic expression for gMie(σ)

It is now straightforward to obtain an algebraic expression for the first-order contribution to the RDF at contact as the integral in Eq. (56) can easily be expressed as function of the first-order terms of the Sutherland fluid

$a_{1}^{S}(\lambda )$
a1S(λ)⁠. The treatment outlined in Sec. III B 1 is followed whereby

\begin{eqnarray}g_{1}(\sigma ) \approx g_{1}(d)&=&\frac{1}{2\pi \epsilon d^{3}}\biggl [ 3\frac{\partial a_{1}}{\partial \rho_s } \nonumber\\&& -\, \mathcal {C}\lambda _{a}x_0^{\lambda _{a}} \frac{ a_{1}^{S}(\eta ;\lambda _{a})+ B(\eta ;\lambda _{a}) }{\rho_s } \nonumber \\[4pt]&& +\, \mathcal {C}\lambda _{r}x_0^{\lambda _{r}} \frac{ a_{1}^{S}(\eta ;\lambda _{r})+ B(\eta ;\lambda _{r}) }{\rho_s } \biggr ] .\end{eqnarray}
g1(σ)g1(d)=12πεd3[3a1ρsCλax0λaa1S(η;λa)+B(η;λa)ρs+Cλrx0λra1S(η;λr)+B(η;λr)ρs].
(64)

The second-order term requires a somewhat lengthier derivation. The details are omitted and we give only the final expression,

\begin{eqnarray}g_{2}(\sigma ) \approx g_{2}(d)= (1+\gamma _c){g}_{2}^{\scriptsize\textit{MCA}}(d) ,\end{eqnarray}
g2(σ)g2(d)=(1+γc)g2MCA(d),
(65)

where

{\fontsize{9.5}{10.5}{\begin{eqnarray}&&\hspace*{-6pt} g_{2}^{\scriptsize\textit{MCA}}(\sigma )\nonumber\\&&\hspace*{-6pt} \quad = \frac{1}{2\pi \epsilon ^2 d^{3}}\biggl [ 3\frac{\partial \frac{a_{2}}{1+\chi }}{\partial \rho_s} - \epsilon K^{HS} \mathcal {C}^2 \lambda _{r}x_0^{2\lambda _{r}} \frac{ a_{1}^{S}(\eta ;2\lambda _{r})+ B(\eta ;2\lambda _{r}) }{\rho_s } \nonumber \\[4pt]&&\hspace*{-6pt} \qquad +\, \epsilon K^{HS} \mathcal {C}^2 (\lambda _{r}+\lambda _{a})x_0^{\lambda _{r}+\lambda _{a}} \frac{ a_{1}^{S}(\eta ;\lambda _{r}+\lambda _{a})+ B(\eta ;\lambda _{r}+\lambda _{a}) }{\rho_s }\nonumber \\[4pt]&&\hspace*{-6pt} \qquad -\, \epsilon K^{HS} \mathcal {C}^2 \lambda _{a}x_0^{2\lambda _{a}} \frac{ a_{1}^{S}(\eta ;2\lambda _{a})+ B(\eta ;2\lambda _{a}) }{\rho_s } \biggr ] .\end{eqnarray}}}
g2MCA(σ)=12πε2d3[3a21+χρsεKHSC2λrx02λra1S(η;2λr)+B(η;2λr)ρs+εKHSC2(λr+λa)x0λr+λaa1S(η;λr+λa)+B(η;λr+λa)ρsεKHSC2λax02λaa1S(η;2λa)+B(η;2λa)ρs].
(66)

The expression for the contact value of the RDF of the Mie fluid is then obtained by substituting Eq. (46), (64), and (65) into Eq. (43). One should note, however, that this algebraic expression is valid only for Mie potentials with values of the repulsive and attractive exponents over the interval which was used in the parameterization of

$a^S_1(\eta ;\lambda )$
a1S(η;λ) [5 ⩽ λ ⩽ 100].

With the Helmholtz free energy and the RDF at contact of the Mie fluids in hand, the traditional SAFT41,42 formalism for the equation of state of fluids of associating chain molecules formed from Mie segments can be implemented in a straightforward manner. We summarize here the expressions for pure substances. The extension to mixtures with the appropriate mixing rules is described in the Appendix.

For a fluid of N chains formed from ms Mie segments there are a total of Ns = msN segments, and the number densities of chain molecules ρ = N/V and of segments ρs = Ns/V are related simply through ρs = msρ. The dimensionless Helmholtz free energy of the fluid a = A/(NkT) can be expressed in terms of the ideal, monomer segment, chain, and associative contributions as

\begin{equation}a=a^{\scriptsize\textit{IDEAL}}+a^{\scriptsize\textit{MONO}}+a^{\scriptsize\textit{CHAIN}}+a^{\scriptsize\textit{ASSOC}}.\end{equation}
a=aIDEAL+aMONO+aCHAIN+aASSOC.
(67)

The ideal-gas term is given by85 

\begin{equation}a^{\scriptsize\textit{IDEAL}}={\rm ln}(\rho \Lambda ^{3})-1 ,\end{equation}
aIDEAL= ln (ρΛ3)1,
(68)

where Λ3 is taken to represent the de Broglie volume which incorporates all of the kinetic (translational, rotational, and vibrational) contributions to the partition function.

The contribution to the Helmholtz free energy for the reference monomeric system (corresponding to the atomized chains) is

\begin{equation}a^{\scriptsize\textit{MONO}}=m_{s}a^{M} ,\end{equation}
aMONO=msaM,
(69)

where aM = AM/(NskT) is the corresponding residual Helmholtz free energy per monomer. As outlined in Sec. III A we have developed a perturbation series expansion in the inverse of the temperature β = 1/(kT) up to third order:

\begin{equation}a^{M}=a^{HS}+\beta a_1+(\beta )^{2}a_2+(\beta )^{3}a_3 ,\end{equation}
aM=aHS+βa1+(β)2a2+(β)3a3,
(70)

where aHS is given by Eq. (11), a1 by Eq. (34), a2 by Eq. (36), and a3 by Eq. (19).

Assuming that the Mie monomers are tangentially bonded at r = σ to form chains of ms segments, then the chain contribution to the Helmholtz free energy can be expressed in the usual Wertheim TPT1 form as47,101

\begin{equation}a^{\scriptsize\textit{CHAIN}}=-(m_{s}-1){\rm ln}g^{Mie}(\sigma ) ,\end{equation}
aCHAIN=(ms1) ln gMie(σ),
(71)

where the RDF at contact is obtained from

\begin{eqnarray}g^{Mie}(\sigma )&=&g^{HS}_d(\sigma ){\rm exp} \big[ \beta \epsilon g_{1}(\sigma )/g^{HS}_d(\sigma ) \nonumber \\[4pt]&&+\, (\beta \epsilon )^2 g_{2}(\sigma )/g^{HS}_d(\sigma ) \big], \\end{eqnarray}
gMie(σ)=gdHS(σ) exp [βεg1(σ)/gdHS(σ)+(βε)2g2(σ)/gdHS(σ)],
(72)

with the zeroth-order term

$g^{HS}_d(\sigma )$
gdHS(σ) given by Eq. (46), the first-order term g1(σ) by Eq. (64), and the second-order term g2(σ) is obtained from Eqs. (65) and (66).

The contribution to the Helmholtz free energy due to association between chains of Mie segments with s bonding site types, na sites of type a, a = 1, …, s, is obtained from the usual Wertheim TPT1 expression:101–103 

\begin{equation}a^{\scriptsize\textit{ASSOC}}=\sum _{a=1}^{s}n_a\left(\ln X_{a}-\frac{1}{2}X_{a}+\frac{1}{2} \right) .\end{equation}
aASSOC=a=1snalnXa12Xa+12.
(73)

Here Xa is the fraction of molecules not bonded at sites of type a, which is obtained from the solution of the following mass-action relation:

\begin{equation}1=X_{a}+\rho X_{a}\sum _{b=1}^{s}n_bX_{b}\Delta _{ab}\ \ \ \ a,b =1,2,\ldots,s.\end{equation}
1=Xa+ρXab=1snbXbΔaba,b=1,2,...,s.
(74)

The function Δab, which characterizes the strength and bonding volume of the association interaction, is given by

\begin{equation}\Delta _{ab}=\int g^{Mie}(r)\big\langle f_{ab}(\mathbf {r}_{ab})\big \rangle {\rm d} r ,\end{equation}
Δab=gMie(r)fab(rab)dr,
(75)

where fab(rab) = exp( − βϕab(rab)) − 1 is the Mayer-f function of the association potential described here with short-ranged square-well potentials:

\begin{equation}\phi _{ab}(\mathbf {r}_{ab})= \left\lbrace \begin{array}{c @{} c}-\epsilon ^{HB}_{ab} &\hspace{28.45274pt} |\mathbf {r}_{ab}| < r^c_{ab} \\[6pt]0 &\hspace{28.45274pt} |\mathbf {r}_{ab}| \ge r^c_{ab}. \\\end{array} \right.\end{equation}
ϕab(rab)=εabHB|rab|<rabc0|rab|rabc.
(76)

The vector rab denotes the centre-centre distance between the association sites a and b, and

$r^{c}_{ab}$
rabc is the range of the association interaction; the square-well bonding sites are placed at a distance
$r^d_{ab}$
rabd
from the centres of the Mie segments. By introducing the square-well site-site interaction into Eq. (75) and performing the angle average of the Mayer-f function one can express Δab as

\begin{equation}\Delta _{ab}= \sigma ^{3} F_{ab} I_{ab} ,\end{equation}
Δab=σ3FabIab,
(77)

where

$F_{ab}={\rm exp} ( \beta \epsilon ^{HB}_{ab} )-1$
Fab= exp (βεabHB)1 and Iab is a dimensionless association integral defined as

\begin{eqnarray}I_{ab} &=& \frac{\pi }{6 \sigma ^3 {r^{d}_{ab}}^{2}} \int _{2r^d_{ab}-r^c_{ab}}^{2r^d_{ab}+r^c_{ab}}g^{Mie}(r)\left( r^c_{ab}+2r^d_{ab}-r \right)^2\nonumber \\[4pt]&&\times \, \big( 2r^c_{ab}-2r^d_{ab}+r \big)r {\rm d}r .\end{eqnarray}
Iab=π6σ3rabd22rabdrabc2rabd+rabcgMie(r)rabc+2rabdr2×2rabc2rabd+rrdr.
(78)

The determination of this integral therefore requires a knowledge of the RDF of the reference Mie fluid over a range of separations around the contact distance σ, so that the expression for the RDF at contact derived in Sec. III B is not sufficient. Even in the specific case of the Lennard-Jones potential, an accurate analytical expression for this function is challenging to obtain.99,104,105

A compromise involves the use of empirical or semi-empirical correlations of the corresponding computer-simulation data for fluids with association sites. The bonding integral has been evaluated numerically for LJ fluids with single SW bonding sites for a specific site position and range of interaction (bonding geometry) by Müller and Gubbins106 using accurate tabulated values for the LJ RDF at different particle separations, the latter obtained from molecular-dynamics simulation. The results of the quadrature were then correlated to a 25-parameter empirical relation in terms of reduced density and temperature. This is, however, only appropriate for the LJ system with the specific bonding geometry, i.e., for the fixed value of range

$r^c_{ab}$
rabc and position
$r^d_{ab}$
rabd
of the association sites relative to the centre of the LJ segment employed in the integration. Another sound treatment of the association in such models has been reported by Tang and Lu107 who proposed a new method to solve the Ornstein-Zernike integral equation in order to obtain the RDF analytically around r = σ. This methodology has proven to be as accurate as the empirical correlation of Ref. 106 but suffers from two limitations: the approach cannot be applied directly in the case of the generic Mie potential; the procedure leads to a complicated expression for the RDF which does not allow for the integral in Eq. (78) to be evaluated analytically. A simpler route is to use a perturbation theory such as WCA83 which leads to a simple expression for the LJ RDF while retaining a satisfactory description of the association integral Iab for a broad range of thermodynamic conditions. This type of analytical treatment can be improved by adding the correct low-density limit (LWCA) in an ad hoc manner.108 Though the LWCA approximation is less accurate than the integral-equation method of Tang and Lu,107 it can easily be applied to any spherically symmetrical segment and for any geometry of the association site.

As a complement to the full numerical solution of the integral in Eq. (78) with an accurate representation of the RDF of the reference Mie fluid (obtained from integral-equation theories or molecular simulation), we also propose a simple approximate expression for the algebraic determination of Iab. Our approach is based on the zeroth-order high-temperature expansion of the RDF, whereby gMie(r) is related to the cavity function

$y^{HS}_d(r)$
ydHS(r) of a hard-sphere fluid of diameter d through

\begin{equation}g^{Mie}(r) \approx g^{HS}_d(r)={\rm exp} \left( -\beta u_{0}(r) \right)y^{HS}_d(r),\end{equation}
gMie(r)gdHS(r)= exp βu0(r)ydHS(r),
(79)

with the reference potential u0(r) given by Eq. (4). By substituting this expression into Eq. (78) we obtain a tractable approximation for the bonding integral Iab of the Mie reference fluid with association sites:

\begin{eqnarray}I_{ab}&=& \frac{\pi }{6 \sigma ^3 {r^d_{ab}}^{2}} \int _{2r^d_{ab}-r^c_{ab}}^{2r^d_{ab}+r^c_{ab}} {\rm exp} \left( -\beta u_{0}(r) \right)y^{HS}_d(r)\nonumber \\[4pt]&&\times \left( r^c_{ab}+2r^d_{ab}-r \right)^2 \left( 2r^c_{ab}-2r^d_{ab}+r \right)r {\rm d}r .\end{eqnarray}
Iab=π6σ3rabd22rabdrabc2rabd+rabc exp βu0(r)ydHS(r)×rabc+2rabdr22rabc2rabd+rrdr.
(80)

When u0(r) is approximated by that of hard spheres of diameter d, uHS(r), the integral can be simplified as

\begin{eqnarray}I_{ab}&=& \frac{\pi }{6 \sigma ^3 {r^d_{ab}}^{2}} \int _{d}^{2r^d_{ab}+r^c_{ab}} g^{HS}_d(r) \nonumber \\[4pt]&&\times\, \big ( r^c_{ab}+ 2r^d_{ab}-r \big )^2 \left( 2r^c_{ab}-2r^d_{ab}+r \right)r {\rm d}r .\end{eqnarray}
Iab=π6σ3rabd2d2rabd+rabcgdHS(r)×rabc+2rabdr22rabc2rabd+rrdr.
(81)

As can be seen from Eq. (81) the evaluation of Iab requires the hard-sphere RDF over a relatively short range. In order to keep a simple analytical form for the association contribution one can assume that

$r^2 g^{HS}_d(r) \approx d^2 g^{HS}_d(d)$
r2gdHS(r)d2gdHS(d) over the interval [d;
$2r^d_{ab}+r^c_{ab}$
2rabd+rabc
]. The approximation is expected to be reasonable for a broad range of densities.102 The bonding integral can then simply be expressed as

\begin{eqnarray}I_{ab} &=& \frac{\pi d^2}{6 \sigma ^3 {r^d_{ab}}^{2}}g^{HS}_d(d) \nonumber \\[4pt]&& \times \, \int _{d}^{2{r^d_{ab}}+{r^c_{ab}}} \frac{ \left( {r^c_{ab}}+2{r^d_{ab}}-r \right)^2 \left( 2{r^c_{ab}}-2{r^d_{ab}}+r \right)}{r} {\rm d}r . \nonumber\!\!\!\!\! \\\end{eqnarray}
Iab=πd26σ3rabd2gdHS(d)×d2rabd+rabcrabc+2rabdr22rabc2rabd+rrdr.
(82)

It is important to note that an implicit dependence on temperature is retained in this expression for Iab through the temperature dependent hard-sphere diameter d, cf. Eq. (7). One can perform the integration to give

\begin{equation}I_{ab} = g^{HS}_d(d)K_{ab} ,\end{equation}
Iab=gdHS(d)Kab,
(83)

where the dimensionless bonding volume is102 

{\fontsize{9.7}{10.7}{\begin{eqnarray}K_{ab}&=&4\pi d^2\big [ \ln \big ( \big(r^c_{ab}+2r^d_{ab}\big)\big/d \big ) \big(6{r^c_{ab}}^3+18 {r^c_{ab}}^{2}r^d_{ab}-24{r^d_{ab}}^3\big) \nonumber \\[4pt]&&+\, \big(r^c_{ab}+2r^d_{ab}-d\big) \big(22{r^d_{ab}}^2-5r^c_{ab}r^d_{ab}-7r^d_{ab}d \nonumber \\[4pt]&&-\, 8{r^c_{ab}}^2+r^c_{ab}d+d^2\big) \big ] \big/ \big(72{r^d_{ab}}^2\sigma ^3\big)\end{eqnarray}}}
Kab=4πd2[lnrabc+2rabd/d6rabc3+18rabc2rabd24rabd3+rabc+2rabdd(22rabd25rabcrabd7rabdd8rabc2+rabcd+d2)]/72rabd2σ3
(84)

The theoretical treatment of the association contribution described in this section complies with the usual Wertheim TPT1 scheme43–48 requiring that: bonding at a given site is independent of bonding at another; multiple bonding at a given site is forbidden; and only dimeric, chain, and network aggregates are formed (with no ring-like or double bonded structures).

The preceding exposition has centred on the development of an improved theory for both the EOS and the contact RDF of monomer Mie fluids as well as the chain contribution and the association integral, collectively referred to as SAFT-VR Mie. In this section we assess and discuss the accuracy of the new theory in describing the thermodynamic properties of systems of chain molecules formed from Mie segments. At this stage the focus is on pure fluids by making use of the expressions developed in Sec. IV. A detailed analysis for mixtures is left for future work.

1. Vapour-liquid equilibria of the Lennard-Jones (12-6) fluid

The vapour-liquid coexistence and vapour pressure curves of the ubiquitous Lennard-Jones fluid determined by solving for the equality of pressure, temperature and chemical potential in the two phases with our reformulated SAFT-VR Mie EOS are shown in Fig. 6 and compared with the essentially exact data obtained by computer simulation. The description obtained with our third-order perturbation expansion of the Helmholtz free energy is also compared with the corresponding results with the expansion truncated at first and second order. As can be seen from the figure, the coexistence densities and vapour pressures obtained with the new SAFT-VR Mie EOS are in excellent agreement with molecular simulation. It is worth noting that the SAFT-VR Mie description is almost of identical quality to that obtained with the EOS of Johnson et al.63 which is based on a correlation of simulation data. One of the key findings of our current work is the striking impact of the third-order perturbation term on the accurate representation of the near-critical region of the VLE, allowing for a significant decrease in the characteristic overshoot of the critical temperature which is obtained with the lower-order expansions. This feature may at first appear rather surprising as it is widely acknowledged that any analytical EOS (including those obtained with a traditional Barker and Henderson84,85 perturbation theory) will exhibit a classical quadratic convergence of the difference in the liquid and vapour coexistence densities near the critical temperature, in disagreement with the experimental near-cubic critical behaviour: the exponent characterizing the classical van der Waals critical behaviour is 1/2 while the universal critical exponent is ∼ 0.325.109–112 Despite this inherent deficiency of classical perturbation expansions, we note here that our third-order expansion yields a quadratic convergence only very close to the critical temperature thus allowing for a very good description over the rest of the near-critical region without the need of a cross-over treatment. We also emphasize that experimentally the cross-over between the classic and universal behaviour is seen only within a tenth of a kelvin of the critical point.

FIG. 6.

The vapour-liquid equilibria of the Lennard-Jones (12-6) fluid: (a) the temperature-density coexistence envelope, where T* = kT/ε and

$\rho _s^\ast\break =N_s\sigma^{3}/V$
ρs*=Nsσ3/V⁠; (b) the vapour-pressure curve, where P* = Pσ2/ε. Comparison of the description with our new SAFT-VR Mie EOS with the Barker and Henderson expansion of the Helmholtz free energy taken to first order (dotted curves), second order (dashed-dotted curves), and third order (continuous curve) terms. The circles represent the molecular-simulation data of Lofti et al.,113 and the dashed curves the description with the equation of state of Johnson et al.63 

FIG. 6.

The vapour-liquid equilibria of the Lennard-Jones (12-6) fluid: (a) the temperature-density coexistence envelope, where T* = kT/ε and

$\rho _s^\ast\break =N_s\sigma^{3}/V$
ρs*=Nsσ3/V⁠; (b) the vapour-pressure curve, where P* = Pσ2/ε. Comparison of the description with our new SAFT-VR Mie EOS with the Barker and Henderson expansion of the Helmholtz free energy taken to first order (dotted curves), second order (dashed-dotted curves), and third order (continuous curve) terms. The circles represent the molecular-simulation data of Lofti et al.,113 and the dashed curves the description with the equation of state of Johnson et al.63 

Close modal

2. Vapour-liquid equilibria and isotherms of Mie (λr-6) and (2λaa) fluids

In Fig. 7(a) the simulated VLE for Mie (λr-6) fluids obtained by Okumura and Yonezawa95 for λr = 8, 10, 15, and 20, by Lofti et al.113 for λr = 12, and by Potoff and Bernard-Brunel96 for λr = 36 are compared with the description with our new SAFT-VR Mie EOS and with that of the previous version of the theory from 2006.72 In the first comparison, only the repulsive exponent of the potential is varied in order to assess the capability of the theory in capturing the non-conformal behaviour that is due to the steepness of the potential. As shown in Fig. 7(a), the SAFT-VR Mie EOS of 2006 provides an excellent description of the coexistence curve of the Mie (10-6) fluid, a fair description for the Mie (12-6) and (8-6) fluids, but poor agreement for the softer potentials such as the Mie (8-6) fluid, especially for the liquid branch. It is also important to point out that this previous version of the theory72 cannot be applied to potentials with exponents that are higher than λr = 12; the reader is referred to Sec. III A 3 for a detailed explanation of the source of this limitation. Our new SAFT-VR Mie EOS is seen to provide a very satisfactory agreement for the full range of repulsive exponents considered, as well as an excellent description of the VLE near the critical region. The SAFT-VR Mie EOS allows for an accurate description of very steeply repulsive potentials such as the Mie (36-6) fluid; one should note that this potential is very short ranged and the VLE envelope appears at very low temperatures. The high-fidelity agreement obtained with the SAFT-VR Mie EOS supports the validity of the form of the third-order expansion proposed in our current work. A slight deterioration in the theory is still, however, apparent for very soft potentials such as the Mie (8-6) fluid. We attribute this discrepancy to the representation of the reference Helmholtz free energy a0 of the Mie potential which is mapped on to that of the hard-sphere reference aHS through the use of a simple density-independent effective diameter d; this approximation is known to lead to some shortcomings for very soft potentials.114 

FIG. 7.

The temperature-density coexistence envelopes of Mie (λra) fluids, where T* = kT/ε and

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V⁠: (a) Mie (λr-6) potentials; and (b) Mie (2λaa) potentials. The description with our SAFT-VR Mie equation of state (continuous curves) is compared to that obtained with the previous version of the theory72 (dashed curves), and to the corresponding molecular-simulation data of Lofti et al.,113 Okumura and Yonezawa,95 Orea et al.,115 and Potoff and Bernard-Brunel96 (circles).

FIG. 7.

The temperature-density coexistence envelopes of Mie (λra) fluids, where T* = kT/ε and

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V⁠: (a) Mie (λr-6) potentials; and (b) Mie (2λaa) potentials. The description with our SAFT-VR Mie equation of state (continuous curves) is compared to that obtained with the previous version of the theory72 (dashed curves), and to the corresponding molecular-simulation data of Lofti et al.,113 Okumura and Yonezawa,95 Orea et al.,115 and Potoff and Bernard-Brunel96 (circles).

Close modal

We also consider the VLE of the (2λaa) Mie fluids. The theoretical description is displayed in Fig. 7(b) along with the Monte Carlo simulation data of Orea et al.115 Four different Mie potentials are considered (with λa = 6, 7, 9, and 12). It can be seen that the new SAFT-VR Mie EOS provides a very accurate description of the coexistence densities, and can be applied simultaneously to very high values of both the repulsive and attractive exponents.

The accuracy of the theory in capturing the vapour pressure of the Mie fluids for different forms of the repulsive interactions is also analysed. We compare the available molecular-simulation data of Okumura and Yonezawa95 (for λr = 10), Lofti et al.113 (for λr = 12), and Potoff and Bernard-Brunel96 (for λr = 14 and 36) with the vapour pressure obtained using our new SAFT-VR Mie EOS (and with the previous version of the theory72). Excellent agreement is found in each case with a slight discrepancy for the Mie (10-6) potential, cf. Fig. 8. The new theory can be applied without loss of accuracy up to extremely steep potentials. From this analysis it is clearly apparent that the present form of the theory with the third-order expansion greatly enhances the accuracy of the description of the exact molecular-simulation results compared to the previous version,72 particularly in the vicinity of the critical point. The third-order term is useful in allowing one to correctly capture the properties of very short-ranged potentials (e.g., with large values of the repulsive exponent) for which the reduced critical temperature is located below T* = 1.

FIG. 8.

The vapour-pressure of the Mie (λr-6) fluids, where T* = kT/ε and P* = Pσ2/ε. The description with our SAFT-VR Mie equation of state (continuous curves) is compared to that obtained with the previous version of the theory72 (dashed curves), and to the corresponding molecular-simulation data of Okumura and Yonezawa,95 Lofti et al.,113 and Potoff and Bernard-Brunel96 (circles).

FIG. 8.

The vapour-pressure of the Mie (λr-6) fluids, where T* = kT/ε and P* = Pσ2/ε. The description with our SAFT-VR Mie equation of state (continuous curves) is compared to that obtained with the previous version of the theory72 (dashed curves), and to the corresponding molecular-simulation data of Okumura and Yonezawa,95 Lofti et al.,113 and Potoff and Bernard-Brunel96 (circles).

Close modal

We have also determined pressure-density isotherms using our new theory for Mie (2λaa) and (λr-6) fluids to compare with the Monte Carlo simulation data reported by Orea et al.115 As can be seen from Fig. 9 the new theory provides an excellent description of the simulation data in all cases, even for “unusual” forms of the potential such as the Mie (32-6) and Mie (24-12) fluids. It is therefore apparent that our SAFT-VR Mie EOS could be used in principle to model the thermodynamic properties of the Girifalco potential114,116–118 (Mie, 29.21-11.83) which has been developed for Buckminsterfullerene, C60.

FIG. 9.

Pressure-density isotherms for Mie (λra) fluids, where T* = kT/ε, P* = Pσ2/ε, and

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V⁠: (a)–(c) Mie (λr-6) potentials; and (d)–(f) Mie (2λaa) potentials. The description obtained with our new SAFT-VR Mie equation of state (continuous curves) is compared to the corresponding molecular-simulation data of Orea et al.115 (circles).

FIG. 9.

Pressure-density isotherms for Mie (λra) fluids, where T* = kT/ε, P* = Pσ2/ε, and

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V⁠: (a)–(c) Mie (λr-6) potentials; and (d)–(f) Mie (2λaa) potentials. The description obtained with our new SAFT-VR Mie equation of state (continuous curves) is compared to the corresponding molecular-simulation data of Orea et al.115 (circles).

Close modal

In this overall assessment we have demonstrated the robustness of our new SAFT-VR Mie EOS in accurately describing the thermodynamic and structural properties of monomeric Mie fluids and the capability of the theory in capturing the non-conformal behaviour induced by changing the repulsive and/or the attractive range of the Mie potential.

3. Vapour-liquid equilibria of Lennard-Jones (12-6) chain fluids

In this section the VLE of fluids of fully flexible chain molecules formed from tangent Lennard-Jones segments as described using our new SAFT-VR Mie EOS is presented and compared to the description obtained with the previous version of the theory,72 as well as the SAFT-VR LJC,65 and the EOSs of Johnson et al.63 The vapour-liquid coexistence densities and vapour pressures for LJ chains formed from ms = 2, 4, and 8 segments are shown in Fig. 10. The coexistence envelopes are compared with the corresponding molecular-simulation data of Lofti et al.113 for the monomer (ms = 1) and of Escobedo and de Pablo119 for the chains, except for the dimer fluid where the more-recent simulation data of Vega et al.120 are used; the simulated vapour pressures for chains of ms = 4 and 8 LJ segments are taken from the work of MacDowell and Blas.121 It is interesting to point out the marked discrepancies with the previous SAFT-VR Mie 2006 and SAFT-VR LJC approaches. This is an important comparison as these EOSs differ in the level at which the perturbation theory is truncated and most critically the expression used for the RDF at contact to describe the chain contribution of the Helmholtz free energy. Such an analysis indicates that the Wertheim TPT1 treatment is particularly sensitive to approximations employed in estimating the pair correlation function for the reference fluid (the monomer Mie fluid in our case). This feature of Wertheim's theory has already been pointed out80 where it was shown that much of the disagreement observed between a TPT1 representation and molecular-simulation data can be attributed to the approximations used to derive the RDF at contact. By contrast, the EOS63 of Johnson et al. involves an empirical correlation of the EOS and contact value of the RDF obtained from exact simulation data for the monomer LJ fluid, and as a consequence represents the best agreement that one could obtain using a SAFT-like equation for these systems. Our new SAFT-VR EOS is derived with a chain term using a RDF treated up to second order and is therefore expected to provide a marked improvement in the description of LJ chains. Essentially, no overshoot of the critical temperature of the LJ dimer is observed. The progressive increase in the over-prediction of the critical temperature found for longer chains can be attributed to the use of the linear approximation to evaluate the many-body distribution function at contact in the chain contribution to the Helmholtz free energy.122 

FIG. 10.

Vapour-liquid equilibria of the fluids of chain molecules formed from ms Lennard-Jones (12-6) segments: (a) the temperature-density coexistence envelope, where T* = kT/ε and

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V⁠; (b) the vapour-pressure curve, where P* = Pσ2/ε. The description with our new SAFT-VR Mie equation of state (continuous curves) is compared to that obtained with the previous version of the theory72 (dashed curves), the SAFT-VR LJC approach65 (dotted curves), the empirical equation of state of Johnson et al.63 (dotted-dashed curves), and to the corresponding molecular-simulation data (circles): ms = 1 (Lofti et al.113), ms = 2 (Vega et al.120), ms = 4 and ms = 8 (Escobedo et al.,119 MacDowell and Blas121).

FIG. 10.

Vapour-liquid equilibria of the fluids of chain molecules formed from ms Lennard-Jones (12-6) segments: (a) the temperature-density coexistence envelope, where T* = kT/ε and

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V⁠; (b) the vapour-pressure curve, where P* = Pσ2/ε. The description with our new SAFT-VR Mie equation of state (continuous curves) is compared to that obtained with the previous version of the theory72 (dashed curves), the SAFT-VR LJC approach65 (dotted curves), the empirical equation of state of Johnson et al.63 (dotted-dashed curves), and to the corresponding molecular-simulation data (circles): ms = 1 (Lofti et al.113), ms = 2 (Vega et al.120), ms = 4 and ms = 8 (Escobedo et al.,119 MacDowell and Blas121).

Close modal

4. Vapour-liquid equilibria of long Lennard-Jones (12-6) chain fluids

The description with our new SAFT-VR Mie EOS for the coexistence densities and vapour pressures of fully flexible chains of ms = 20, 50, and 100 tangent LJ segments is compared in Fig. 11 with the corresponding VLE obtained using the semi-empirical EOS of Johnson et al.63 As can be seen, both EOSs provide equivalent satisfactory descriptions of the fluid-phase equilibria of these long chains. Discrepancies between the two approaches are apparent near the critical region, where a slightly larger over-prediction of the critical temperature is observed with the SAFT-VR Mie EOS, together with a deterioration of the vapour pressure with increasing chain length. These deviations are expected for very long chains and are a natural consequence of the linear approximation used to represent the many-body correlations (in terms of the pair distribution function) within the Wertheim TPT1 treatment of the chain contribution. In order to address this deficiency of the TPT1 approximation, Blas and Vega123 have built on the work of Ghonasgi and Chapman124 and Chang and Sandler,125 demonstrating that an improved description of the coexistence densities and vapour pressure of chain fluids of moderate length can be obtained by considering a dimer (instead of a monomer) reference fluid, thereby including more structural information in the chain contribution to the Helmholtz free energy.

FIG. 11.

Vapour-liquid equilibria of fluids of long chains formed from ms Lennard-Jones (12-6) segments: (a) the temperature-density coexistence envelope, where T* = kT/ε and

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V⁠; (b) the vapour-pressure curve, where P* = Pσ2/ε. The description with our new SAFT-VR Mie equation of state (continuous curves) is compared to that obtained with the empirical approach of Johnson et al.63 (dotted-dashed curves), and to the corresponding molecular-simulation data (circles): ms = 2 (Vega et al.120); ms = 4 and ms = 8 (Escobedo et al.119); ms = 20, ms = 50, and ms = 100 (Sheng et al.148).

FIG. 11.

Vapour-liquid equilibria of fluids of long chains formed from ms Lennard-Jones (12-6) segments: (a) the temperature-density coexistence envelope, where T* = kT/ε and

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V⁠; (b) the vapour-pressure curve, where P* = Pσ2/ε. The description with our new SAFT-VR Mie equation of state (continuous curves) is compared to that obtained with the empirical approach of Johnson et al.63 (dotted-dashed curves), and to the corresponding molecular-simulation data (circles): ms = 2 (Vega et al.120); ms = 4 and ms = 8 (Escobedo et al.119); ms = 20, ms = 50, and ms = 100 (Sheng et al.148).

Close modal

5. Vapour-liquid equilibria of Mie (λr-6) chain fluids

As no molecular-simulation data have to our knowledge been reported for fluids of chain molecules formed from Mie segments, we have performed direct molecular-dynamics simulations of the VLE for two model Mie (λr-6) flexible chain molecules: three tangent segments (ms = 3) interacting through a Mie (22-6) potential; and six tangent segments (ms = 6) interacting through a Mie (16-6) potential. For details of the simulations the reader is directed to our series of papers126–128 on the use of the Mie potential in the development of coarse-grained models. A comparison of the description obtained for VLE of these Mie chain fluids with our SAFT-VR EOS and the simulated values is shown in Fig. 12. The good agreement obtained in both cases provides a strong confirmation of the adequacy of the expression derived for the contact value of the RDF of Mie fluids when considering potentials of variable repulsive ranges. No deterioration of the description is apparent for chains with steep repulsive segment interactions such as the Mie (22-6) potential. We emphasize again that this system is a challenge for any high-temperature perturbation theory because as the repulsive exponent λr is increased (i.e., as the range of the interaction is shortened), the VLE envelope is found at progressively lower reduced temperatures. It is therefore necessary to make use of a very accurate description of both EOS and the RDF of the reference Mie monomer fluid, which our SAFT-VR Mie treatment gratifyingly provides. As will be shown in Sec. V B an interaction potential with moderately large values of the repulsive exponent turns out to be very useful in describing simultaneously the VLE and the second-derivative thermodynamic properties of long n-alkanes.

FIG. 12.

The vapour-liquid temperature-density coexistence envelope of fluids of chains formed from ms Mie (λra) segments, where T* = kT/ε and

$\rho _s^\ast = N_s\sigma^3/V$
ρs*=Nsσ3/V⁠. The description with our new SAFT-VR Mie equation of state (continuous curves) is compared to the corresponding molecular-simulation data obtained in our current work (circles).

FIG. 12.

The vapour-liquid temperature-density coexistence envelope of fluids of chains formed from ms Mie (λra) segments, where T* = kT/ε and

$\rho _s^\ast = N_s\sigma^3/V$
ρs*=Nsσ3/V⁠. The description with our new SAFT-VR Mie equation of state (continuous curves) is compared to the corresponding molecular-simulation data obtained in our current work (circles).

Close modal

6. Vapour-liquid equilibria of associating LJ (12-6) fluids

When a soft-core reference interaction is used to model the particles, association sites can interact even when fully embedded within the molecular core. An original feature of our current approach lies in the approximate algebraic expression for the association integral (cf. Eqs. (77), (83), and (84)) which allows one to take into account different geometries of the sites, as opposed to the fixed position (

$r^d_{ab} = 0.4\sigma$
rabd=0.4σ and
$r^c_{ab} = 0.2\sigma$
rabc=0.2σ
) of the semi-empirical correlated expressions.106 Unfortunately, molecular-simulation results for different geometries of the sites are not readily available so we will assess only the association geometry used in the existing correlation.106 In Fig. 13 we display the vapour-liquid coexistence densities for LJ particles with a single embedded square-well association site for two different values of the site-site association energy (
$\epsilon ^{HB}_{ab}=10\epsilon$
εabHB=10ε
and
$\epsilon ^{HB}_{ab}=20\epsilon$
εabHB=20ε
). For the dimerizing LJ system with an association energy of
$\epsilon ^{HB}_{ab}=10\epsilon$
εabHB=10ε
, the agreement in the coexistence envelope between theory and simulation is very good, while small discrepancies are apparent for the strongly associating system with
$\epsilon ^{HB}_{ab}\break =20\epsilon$
εabHB=20ε
. For these associating systems the description with the SAFT-VR Mie EOS turns out to be of similar accuracy to that obtained with other perturbation theories for the reference radial distribution function, such as the low density WCA recipe.80 

FIG. 13.

The vapour-liquid temperature-density coexistence envelope of fluids of associating Lennard-Jones (12-6) molecules with one association site, where T* = kT/ε and

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V⁠. The sites are placed a distance
$r^d_{ab} = 0.4\sigma$
rabd=0.4σ
interacting with a energy of
$\epsilon ^{HB}_{ab}$
εabHB
when they are within a range of
$r^c_{ab} = 0.2\sigma$
rabc=0.2σ
. The description with our new SAFT-VR Mie equation of state (continuous curves) is compared to the corresponding molecular-simulation data (circles).

FIG. 13.

The vapour-liquid temperature-density coexistence envelope of fluids of associating Lennard-Jones (12-6) molecules with one association site, where T* = kT/ε and

$\rho _s^\ast=N_s\sigma^3/V$
ρs*=Nsσ3/V⁠. The sites are placed a distance
$r^d_{ab} = 0.4\sigma$
rabd=0.4σ
interacting with a energy of
$\epsilon ^{HB}_{ab}$
εabHB
when they are within a range of
$r^c_{ab} = 0.2\sigma$
rabc=0.2σ
. The description with our new SAFT-VR Mie equation of state (continuous curves) is compared to the corresponding molecular-simulation data (circles).

Close modal

In our current formulation of the SAFT-VR Mie EOS, molecules are treated in a homonuclear fashion (i.e., all of the segments in the chain model representing the molecule are taken to be identical). We apply the theory to selected compounds, including members of three homologous series (n-alkanes, n-perfluoroalkanes, and n-alkanols), carbon dioxide, benzene, and toluene. The methodology commonly applied to the determination of the molecular parameters of SAFT-like EOS is to estimate them from experimental saturation pressures Psat and liquid densities ρsat. The task involves a minimization of a relative least-squares objective function consisting of the appropriate sums of individual residuals:

\begin{eqnarray}\min _{\theta } F &=&\frac{w_1}{N_{P_{{\rm{sat}}}}}\sum _{i=1}^{N_{P_{{\rm{sat}}}}} \left[\frac{P^{{{\rm exp}}}_{{\rm{sat}},i}(T_i) -P^{{\rm{calc}}}_{{\rm{sat}},i}(T_i;\theta )}{P^{{{\rm exp}}}_{{\rm{sat}},i}(T_i)} \right]^2 \nonumber \\[4pt]&&+\, \frac{w_2}{N_{\rho _{{\rm{sat}}}}}\sum _{i=1}^{N_{\rho _{{\rm{sat}}}}} \left[\frac{\rho ^{{{\rm exp}}}_{{\rm{sat}},i}(T_i) -\rho ^{{\rm{calc}}}_{{\rm{sat}},i}(T_i;\theta )}{\rho ^{{{\rm exp}}}_{{\rm{sat}},i}(T_i)} \right]^2 , \nonumber \!\!\!\!\! \\\end{eqnarray}
minθF=w1NP sat i=1NP sat P sat ,i exp (Ti)P sat ,i calc (Ti;θ)P sat ,i exp (Ti)2+w2Nρ sat i=1Nρ sat ρ sat ,i exp (Ti)ρ sat ,i calc (Ti;θ)ρ sat ,i exp (Ti)2,
(85)

where θ is the vector of parameters of the intermolecular model, and

$N_{P_{{\rm{sat}}}}$
NP sat and
$N_{\rho _{{\rm{sat}}}}$
Nρ sat
are the numbers of vapour pressure and saturated liquid density data points, respectively; the superscripts exp and calc refer to experimental data points and calculated values. The saturation values are calculated with the EOS by determining the vapour and liquid densities that satisfy the equality in the pressure P = −(∂A/∂V)N, T and chemical potential μ = (∂A/∂N)V, T in both phases at a given temperature T. Equal weights are generally chosen for both properties in this objective function so that w1 = w2 = 1. Other properties can be obtained from the standard thermodynamic relations involving the appropriate derivatives of the Helmholtz free energy.

It has, however, been demonstrated72 that the inclusion of single-phase condensed-liquid density ρ and speed of sound u data often leads to a better balance of the different thermophysical properties. This has proven to be very useful in determining intermolecular potential parameters, and in the particular case of the Mie interaction, physically relevant values of the repulsive exponent λr are obtained. The objective function is expressed as

\begin{eqnarray}\hspace*{-10pt}\min _{\theta } F&=&\frac{w_1}{N_{P_{{\rm{sat}}}}}\sum _{i=1}^{N_{P_{{\rm{sat}}}}} \left[\frac{P^{{{\rm exp}}}_{{\rm{sat}},i}(T_i) -P^{{\rm{calc}}}_{{\rm{sat}},i}(T_i;\theta )}{P^{{{\rm exp}}}_{{\rm{sat}},i}(T_i)} \right]^2 \nonumber \\[4pt]&&+ \, \frac{w_2}{N_{\rho _{{\rm{sat}}}}}\sum _{i=1}^{N_{\rho _{{\rm{sat}}}}} \left[\frac{\rho ^{{{\rm exp}}}_{{\rm{sat}},i}(T_i) -\rho ^{{\rm{calc}}}_{{\rm{sat}},i}(T_i;\theta )}{\rho ^{{{\rm exp}}}_{{\rm{sat}},i}(T_i)} \right]^2 \nonumber \\[4pt]&&+\, \frac{w_3}{N_{\rho }}\sum _{i=1}^{N_{\rho }} \left[\frac{\rho ^{{{\rm exp}}}_{i}(T_i,P_i) -\rho ^{{\rm{calc}}}_{i}(T_i,P_i;\theta )}{\rho ^{{{\rm exp}}}_{i}(T_i,P_i)} \right]^2 \nonumber \\[4pt]&&+\, \frac{w_4}{N_{u}}\sum _{i=1}^{N_{u}} \left[\frac{u^{{{\rm exp}}}_{i}(T_i,P_i) -u^{{\rm{calc}}}_{i}(T_i,P_i;\theta )}{u^{{{\rm exp}}}_{i}(T_i,P_i)} \right]^2,\hspace*{15pt}\end{eqnarray}
minθF=w1NP sat i=1NP sat P sat ,i exp (Ti)P sat ,i calc (Ti;θ)P sat ,i exp (Ti)2+w2Nρ sat i=1Nρ sat ρ sat ,i exp (Ti)ρ sat ,i calc (Ti;θ)ρ sat ,i exp (Ti)2+w3Nρi=1Nρρi exp (Ti,Pi)ρi calc (Ti,Pi;θ)ρi exp (Ti,Pi)2+w4Nui=1Nuui exp (Ti,Pi)ui calc (Ti,Pi;θ)ui exp (Ti,Pi)2,
(86)

where Nρ and Nu are the number of condensed-liquid density and speed of sound data points, and w3 and w4 are the corresponding weights. Several values of the different weight factors have been tested, and a good balance in the accuracies of the description of the different estimated properties is obtained using w1 = 1, w2 = 1, w3 = 1, and w4 = 1/4. The resulting Mie intermolecular parameters for the substances examined in our current study with the SAFT-VR Mie EOS are reported in Table III.

Table III.

The Mie (λra) molecular parameters of some selected compounds estimated from the experimental vapour-liquid equilibria and single-phase data with our new SAFT-VR Mie equation of state. The value of ms denotes the number of Mie segments making up the chains, σ characterizes the diameter of the segment, ε the strength of the interaction, λr and λa the repulsive and attractive exponents; for molecules with association sites,

$r^c_{ab}$
rabc characterizes the range of the site-site interaction,
$\epsilon ^{HB}_{ab}$
εabHB
the strength of the bonding, and
$r^d_{ab}= 0.4\sigma$
rabd=0.4σ
is the position of the site in all cases. The n-alcohols are modeled with a two-site association scheme corresponding to one donor and one acceptor site (denoted as 2B in the classification scheme of Huang and Radosz136).

Substancemsσ/Å(ε/k)/Kλrλa
$r^c_{ab}/\sigma$
rabc/σ
$(\epsilon ^{HB}_{ab}/k)$
(εabHB/k)
/K
Methane 1.0000 3.7412 153.36 12.650 … … 
Ethane 1.4373 3.7257 206.12 12.400 … … 
Propane 1.6845 3.9056 239.89 13.006 … … 
n-butane 1.8514 4.0887 273.64 13.650 … … 
n-pentane 1.9606 4.2928 321.94 15.847 … … 
n-hexane 2.1097 4.4230 354.38 17.203 … … 
n-heptane 2.3949 4.4282 358.51 17.092 … … 
n-octane 2.6253 4.4696 369.18 17.378 … … 
n-nonane 2.8099 4.5334 387.55 18.324 … … 
n-decane 2.9976 4.5890 400.79 18.885 … … 
n-dodecane 3.2519 4.7484 437.72 20.862 … … 
n-pentadecane 3.9325 4.7738 444.51 20.822 … … 
n-eicosane 4.8794 4.8788 475.76 22.926 … … 
Methanol 1.5283 3.3063 167.72 8.6556 0.41314 2904.7 
Ethanol 1.9600 3.4914 168.15 7.6134 0.34558 2833.7 
Propan-1-ol 2.3356 3.5612 227.66 10.179 0.35377 2746.2 
n-butan-1-ol 2.4377 3.7856 278.92 11.660 0.32449 2728.1 
Perfluoromethane 1.0000 4.3372 232.62 42.553 5.1906 … … 
Perfluoroethane 1.8529 3.9336 211.46 19.192 5.7506 … … 
Perfluoropropane 1.9401 4.2983 263.26 22.627 5.7506 … … 
n-perfluorobutane 2.1983 4.4495 290.49 24.761 5.7506 … … 
n-perfluoropentane 2.3783 4.6132 328.56 29.750 5.7506 … … 
n-perfluorohexane 2.5202 4.7885 349.30 30.741 5.7506 … … 
Fluorine 1.3211 2.9554 96.268 11.606 … … 
Carbon dioxide 1.5000 3.1916 231.88 27.557 5.1646 … … 
Benzene 1.9163 4.0549 372.59 14.798 … … 
Toluene 1.9977 4.2777 409.73 16.334 … … 
Substancemsσ/Å(ε/k)/Kλrλa
$r^c_{ab}/\sigma$
rabc/σ
$(\epsilon ^{HB}_{ab}/k)$
(εabHB/k)
/K
Methane 1.0000 3.7412 153.36 12.650 … … 
Ethane 1.4373 3.7257 206.12 12.400 … … 
Propane 1.6845 3.9056 239.89 13.006 … … 
n-butane 1.8514 4.0887 273.64 13.650 … … 
n-pentane 1.9606 4.2928 321.94 15.847 … … 
n-hexane 2.1097 4.4230 354.38 17.203 … … 
n-heptane 2.3949 4.4282 358.51 17.092 … … 
n-octane 2.6253 4.4696 369.18 17.378 … … 
n-nonane 2.8099 4.5334 387.55 18.324 … … 
n-decane 2.9976 4.5890 400.79 18.885 … … 
n-dodecane 3.2519 4.7484 437.72 20.862 … … 
n-pentadecane 3.9325 4.7738 444.51 20.822 … … 
n-eicosane 4.8794 4.8788 475.76 22.926 … … 
Methanol 1.5283 3.3063 167.72 8.6556 0.41314 2904.7 
Ethanol 1.9600 3.4914 168.15 7.6134 0.34558 2833.7 
Propan-1-ol 2.3356 3.5612 227.66 10.179 0.35377 2746.2 
n-butan-1-ol 2.4377 3.7856 278.92 11.660 0.32449 2728.1 
Perfluoromethane 1.0000 4.3372 232.62 42.553 5.1906 … … 
Perfluoroethane 1.8529 3.9336 211.46 19.192 5.7506 … … 
Perfluoropropane 1.9401 4.2983 263.26 22.627 5.7506 … … 
n-perfluorobutane 2.1983 4.4495 290.49 24.761 5.7506 … … 
n-perfluoropentane 2.3783 4.6132 328.56 29.750 5.7506 … … 
n-perfluorohexane 2.5202 4.7885 349.30 30.741 5.7506 … … 
Fluorine 1.3211 2.9554 96.268 11.606 … … 
Carbon dioxide 1.5000 3.1916 231.88 27.557 5.1646 … … 
Benzene 1.9163 4.0549 372.59 14.798 … … 
Toluene 1.9977 4.2777 409.73 16.334 … … 

In the case of the n-alkane series, we find rather different values for the repulsive exponent λr with our new SAFT-VR Mie EOS than were found with the previous version of the theory.72 The new SAFT-VR Mie models are consistent with the force fields developed in the simulation studies of Potoff and Bernard-Brunel96 who found that the non-bonded interactions between two CH2 groups in n-alkanes can be accurately modeled with a steeper repulsion than the classical Lennard-Jones (12-6) potential. The SAFT-VR Mie and Potoff and Bernard-Brunel96 models are not, however, identical: in SAFT-VR Mie the bending and torsional potentials are not included explicitly and the non-bonded interactions are the same between all the spherical segments while a heteronuclear model is used by Potoff and Bernard-Brunel;96 we also note that the theory is developed for a tangent-chain model, while the simulation force field is for highly overlapping segments. The difference in the repulsive exponent obtained with the SAFT-VR Mie EOS for eicosane (λr = 22.926) and the value of Potoff and Bernard-Brunel96r = 16) can also be attributed to the optimization procedure used by these authors who did not include vapour-pressure data in the development of the force field. The percentage average absolute deviation AAD(%) of the theoretical description from the experimental data for the vapour pressure Psat, saturated-liquid densities ρsat, enthalpy of vaporization ΔHv, as well as properties in the condensed-liquid region, such as the density ρ, speed of sound u, and isobaric heat capacity Cp, are reported in Table IV; the reader is referred to Table V for details of the thermodynamic conditions considered for each substance and the sources of the experimental data.129–131 It is worth noting the very small overall deviation obtained for both vapour-liquid equilibria and second-derivative thermodynamic properties. A separate comparison of the corresponding values of the critical temperature, pressure, and density is made in Table VI. A key advantage of the formulation of the theory to third order in the perturbation expansion is the significant reduction of the overshoot of the critical temperature. This is illustrated in Fig. 14 where the experimental vapour pressures, saturated densities, and enthalpies of vaporization of the n-alkanes from methane to n-decane are compared to the description with the SAFT-VR Mie EOS. The new theory is also significantly more accurate at low temperatures than the previous SAFT-VR Mie EOS developed in 2006.72 We should note that for some low-temperature isotherms, more than three volume roots can sometimes exist at a given pressure as a result of the non-cubic nature of the perturbation expansion in density. These additional roots do not compromise the adequacy of the EOS as they are found inside the phase coexistence regions and correspond to unstable/metastable maxima/local minima in the Gibbs free energy of the system. Care should however be taken to ensure that the numerical solver used to identify volume roots converges to the global solution of the Gibbs free energy of the system (corresponding to the smallest and largest density roots). In our experience with the SAFT-VR Mie EOS, we have not found any instance of unphysical behaviour, such as the two separate regions of vapour-liquid equilibrium found with PC-SAFT.132–134 

Table IV.

Average absolute deviation AAD(%) of the description obtained with our new SAFT-VR MIE equation of state compared to experimental data for the vapour pressure Psat, the saturated-liquid density ρsat, the enthalpy of vaporization ΔHv, single-phase condensed-liquid density ρ, speed of sound u, and isobaric heat capacity Cp. The definition of the AAD(%) is given in the caption of Table I.

 AAD(%) 
SubstancePsatρsatΔHvρuCp
Methane 0.63 0.78 2.86 1.02 1.29 4.73 
Ethane 0.25 0.83 2.09 0.21 1.84 1.37 
Propane 0.14 0.73 1.97 0.14 2.42 1.40 
n-butane 0.25 0.57 1.87 0.13 2.60 1.07 
n-pentane 0.84 0.40 2.07 0.17 1.56 0.25 
n-hexane 1.12 0.25 2.44 0.15 0.88 1.02 
n-heptane 0.72 0.43 2.24 0.15 1.35 0.79 
n-octane 0.82 0.47 1.72 0.11 0.41 1.17 
n-nonane 0.76 1.46 1.87 0.22 0.59 0.76 
n-decane 0.99 0.49 2.07 1.65 2.48 0.79 
n-dodecane 0.71 0.50 1.89 0.09 0.82 1.37 
n-pentadecane 1.76 0.60 4.09 0.18 0.58 1.55 
n-eicosane 1.86 0.96 1.20 0.36 1.33 0.74 
Methanol 1.05 0.11 1.48 0.21 2.23 8.82 
Ethanol 1.56 0.29 2.51 0.19 3.35 2.69 
Propan-1-ol 0.25 0.13 1.26 0.17 1.68 1.76 
n-butan-1-ol 1.36 0.29 1.40 0.22 4.35 2.09 
Perfluoromethane 0.92 0.49 1.48 1.08 7.38 0.92 
Perfluoroethane 0.32 0.46 1.32 0.95 1.93 0.61 
Perfluoropropane 0.86 0.3- 0.93 0.24 6.64 0.82 
n-perfluorobutane 0.50 1.75 0.69 1.34 5.79 0.85 
n-perfluoropentane 0.22 0.78 0.86 0.39 2.90 0.42 
n-perfluorohexane 0.69 1.45 … 0.73 … … 
Fluorine 0.47 0.52 1.65 0.31 4.45 1.33 
Carbon dioxide 0.40 1.18 3.26 1.01 10.70 3.49 
Benzene 0.95 0.37 2.58 0.35 1.41 3.03 
Toluene 1.44 0.49 2.41 1.87 2.19 2.15 
 AAD(%) 
SubstancePsatρsatΔHvρuCp
Methane 0.63 0.78 2.86 1.02 1.29 4.73 
Ethane 0.25 0.83 2.09 0.21 1.84 1.37 
Propane 0.14 0.73 1.97 0.14 2.42 1.40 
n-butane 0.25 0.57 1.87 0.13 2.60 1.07 
n-pentane 0.84 0.40 2.07 0.17 1.56 0.25 
n-hexane 1.12 0.25 2.44 0.15 0.88 1.02 
n-heptane 0.72 0.43 2.24 0.15 1.35 0.79 
n-octane 0.82 0.47 1.72 0.11 0.41 1.17 
n-nonane 0.76 1.46 1.87 0.22 0.59 0.76 
n-decane 0.99 0.49 2.07 1.65 2.48 0.79 
n-dodecane 0.71 0.50 1.89 0.09 0.82 1.37 
n-pentadecane 1.76 0.60 4.09 0.18 0.58 1.55 
n-eicosane 1.86 0.96 1.20 0.36 1.33 0.74 
Methanol 1.05 0.11 1.48 0.21 2.23 8.82 
Ethanol 1.56 0.29 2.51 0.19 3.35 2.69 
Propan-1-ol 0.25 0.13 1.26 0.17 1.68 1.76 
n-butan-1-ol 1.36 0.29 1.40 0.22 4.35 2.09 
Perfluoromethane 0.92 0.49 1.48 1.08 7.38 0.92 
Perfluoroethane 0.32 0.46 1.32 0.95 1.93 0.61 
Perfluoropropane 0.86 0.3- 0.93 0.24 6.64 0.82 
n-perfluorobutane 0.50 1.75 0.69 1.34 5.79 0.85 
n-perfluoropentane 0.22 0.78 0.86 0.39 2.90 0.42 
n-perfluorohexane 0.69 1.45 … 0.73 … … 
Fluorine 0.47 0.52 1.65 0.31 4.45 1.33 
Carbon dioxide 0.40 1.18 3.26 1.01 10.70 3.49 
Benzene 0.95 0.37 2.58 0.35 1.41 3.03 
Toluene 1.44 0.49 2.41 1.87 2.19 2.15 
Table V.

The experimental data considered in the determination of the AAD(%) (cf. Table I for its definition) for each compound reported in Table IV.

 Vapour-liquid equilibriumCompressed liquid phase 
 T range/KT range/KP range/MPa 
SubstanceρsatPsatΔHvρuCpρuCpRef.
Methane 95–171 95–171 95–171 100–160 100–170 100–150 10–100 10–50 10–50 129  
Ethane 137–274 137–274 137–274 174–220 174–220 174–220 10–60 10–60 10–60 129  
Propane 166–332 166–332 166–332 210–290 210–290 210–290 10–60 10–60 10–60 129  
n-butane 191–382 191–382 191–382 240–340 240–340 240–340 10–60 10–60 10–60 129  
n-pentane 211–422 211–422 211–422 280–360 280–360 280–360 11–80 11–80 10–60 129  
n-hexane 228–456 228–456 228–456 293–373 293–373 293–373 0.1–150 0.1–150 0.1–90 129  
n-heptane 243–486 243–486 243–486 293–373 293–373 293–373 0.1–150 0.1–150 0.1–90 129  
n-octane 255–511 255–511 255–511 293–373 293–373 293–373 5–150 0.1–100 10–100 129  
n-nonane 267–535 267–535 267–535 293–373 293–413 293–413 5–150 0.1–140 10–150 129  
n-decane 277–555 277–555 277–555 293–373 303–413 293–413 5–150 0.1–60 10–150 129  
n-dodecane 296–592 296–592 296–592 293–373 293–473 193–413 5–150 0.1–140 10–150 129  
n-pentadecane 318–637 318–637 298–373 293–383 293–383 313–373 0.1–150 0.1–150 0.1–100 130  
n-eicosane 318–637 318–637 298–373 293–383 293–383 313–373 0.1–150 0.1–150 0.1–100 130  
Methanol 256–461 256–461 298–458 310–430 303–373 248–423 5–100 0.1–50 0.5–12 129  
Ethanol 256–462 256–462 298–458 298–460 273–333 325–460 2–50 0.1–100 4.9–4.9 131  
Propan-1-ol 268–483 268–483 298–383 298–423 273–333 325–470 2–50 0.1–100 6.4–6.4 131  
n-butan-1-ol 281–506 296–476 298–418 298–448 303–453 278–338 1–50 0.1–100 0.1–0.1 131  
Perfluoromethane 113–230 113–230 113–230 140–400 140–400 140–400 0.1–50 0.1–50 0.1–50 129  
Perfluoroethane 173–263 173–263 113–263 180–380 180–380 180–380 0.1–50 0.1–50 0.1–50 129  
Perfluoropropane 130–310 130–310 130–310 150–300 150–300 150–300 0.1–20 0.1–20 0.1–20 129  
n-perfluorobutane 190–350 190–350 190–350 200–320 200–320 200–320 0.1–30 0.1–30 0.1–30 129  
n-perfluoropentane 200–380 200–380 200–380 270–350 270–350 270–350 0.1–30 0.1–30 0.1–30 129  
n-perfluorohexane 273–371 304–430 … 280–313 … … 0.1–40 … … 129  
Fluorine 54–129 54–129 54–129 63–133 63–133 63–133 0.1–20 0.1–20 0.1–20 129  
Carbon dioxide 218–273 218–273 218–273 231–400 231–400 231–400 0.1–50 0.1–50 0.1–50 129  
Benzene 279–514 279–514 279–514 293–453 293–453 293–453 0.1–78 0.1–78 0.1–78 129  
Toluene 178–538 178–538 178–538 293–453 293–453 293–453 0.1–100 0.1–100 0.1–100 129  
 Vapour-liquid equilibriumCompressed liquid phase 
 T range/KT range/KP range/MPa 
SubstanceρsatPsatΔHvρuCpρuCpRef.
Methane 95–171 95–171 95–171 100–160 100–170 100–150 10–100 10–50 10–50 129  
Ethane 137–274 137–274 137–274 174–220 174–220 174–220 10–60 10–60 10–60 129  
Propane 166–332 166–332 166–332 210–290 210–290 210–290 10–60 10–60 10–60 129  
n-butane 191–382 191–382 191–382 240–340 240–340 240–340 10–60 10–60 10–60 129  
n-pentane 211–422 211–422 211–422 280–360 280–360 280–360 11–80 11–80 10–60 129  
n-hexane 228–456 228–456 228–456 293–373 293–373 293–373 0.1–150 0.1–150 0.1–90 129  
n-heptane 243–486 243–486 243–486 293–373 293–373 293–373 0.1–150 0.1–150 0.1–90 129  
n-octane 255–511 255–511 255–511 293–373 293–373 293–373 5–150 0.1–100 10–100 129  
n-nonane 267–535 267–535 267–535 293–373 293–413 293–413 5–150 0.1–140 10–150 129  
n-decane 277–555 277–555 277–555 293–373 303–413 293–413 5–150 0.1–60 10–150 129  
n-dodecane 296–592 296–592 296–592 293–373 293–473 193–413 5–150 0.1–140 10–150 129  
n-pentadecane 318–637 318–637 298–373 293–383 293–383 313–373 0.1–150 0.1–150 0.1–100 130  
n-eicosane 318–637 318–637 298–373 293–383 293–383 313–373 0.1–150 0.1–150 0.1–100 130  
Methanol 256–461 256–461 298–458 310–430 303–373 248–423 5–100 0.1–50 0.5–12 129  
Ethanol 256–462 256–462 298–458 298–460 273–333 325–460 2–50 0.1–100 4.9–4.9 131  
Propan-1-ol 268–483 268–483 298–383 298–423 273–333 325–470 2–50 0.1–100 6.4–6.4 131  
n-butan-1-ol 281–506 296–476 298–418 298–448 303–453 278–338 1–50 0.1–100 0.1–0.1 131  
Perfluoromethane 113–230 113–230 113–230 140–400 140–400 140–400 0.1–50 0.1–50 0.1–50 129  
Perfluoroethane 173–263 173–263 113–263 180–380 180–380 180–380 0.1–50 0.1–50 0.1–50 129  
Perfluoropropane 130–310 130–310 130–310 150–300 150–300 150–300 0.1–20 0.1–20 0.1–20 129  
n-perfluorobutane 190–350 190–350 190–350 200–320 200–320 200–320 0.1–30 0.1–30 0.1–30 129  
n-perfluoropentane 200–380 200–380 200–380 270–350 270–350 270–350 0.1–30 0.1–30 0.1–30 129  
n-perfluorohexane 273–371 304–430 … 280–313 … … 0.1–40 … … 129  
Fluorine 54–129 54–129 54–129 63–133 63–133 63–133 0.1–20 0.1–20 0.1–20 129  
Carbon dioxide 218–273 218–273 218–273 231–400 231–400 231–400 0.1–50 0.1–50 0.1–50 129  
Benzene 279–514 279–514 279–514 293–453 293–453 293–453 0.1–78 0.1–78 0.1–78 129  
Toluene 178–538 178–538 178–538 293–453 293–453 293–453 0.1–100 0.1–100 0.1–100 129  
Table VI.

Critical temperature, pressure, and density of selected compounds (cf. Table III for the molecular parameters) obtained with our new SAFT-VR MIE equation of state compared to experimental data.129 

 SAFT-VR MieExperimental data
SubstanceTc/KPc/MPaρc/kg m−3Tc/KPc /MPaρc /kg m−3
Methane 195.30 5.15 154.45 190.6 4.6 162.7 
Ethane 311.38 5.49 205.84 305.3 4.9 206.6 
Propane 376.20 4.77 219.98 369.8 4.2 220.5 
n-butane 432.68 4.27 228.04 425.1 3.8 228.0 
n-pentane 476.44 3.81 238.14 469.7 3.4 232.0 
n-hexane 515.29 3.44 241.16 507.8 3.0 233.2 
n-heptane 547.33 3.06 233.76 540.1 2.7 232.0 
n-octane 576.72 2.77 227.83 569.3 2.5 235.0 
n-nonane 602.20 2.52 224.89 594.6 2.3 232.1 
n-decane 626.37 2.31 219.85 617.7 2.1 233.3 
n-dodecane 668.75 1.99 214.26 658.2 1.8 221.4 
n-pentadecane 720.98 1.61 194.23 708.0 1.5 212.4 
n-eicosane 786.33 1.21 174.91 768.0 1.1 241.6 
Methanol 547.73 12.21 260.59 512.6 8.1 328.4 
Ethanol 554.41 8.83 247.96 514.0 6.3 276.4 
Propan-1-ol 560.43 6.89 256.34 536.9 5.2 275.2 
n-butan-1-ol 583.97 5.50 253.86 562.0 4.5 270.5 
Perfluoromethane 232.77 4.14 644.29 227.5 3.8 625.6 
Perluoroethane 295.46 3.24 634.96 293.0 3.0 613.3 
Perfluoropropane 347.88 2.79 648.71 345.0 2.6 628.0 
n-perfluorobutane 386.86 2.38 635.61 386.3 2.3 599.8 
n-perfluoropentane 421.36 2.13 634.72 420.6 2.0 609.5 
Fluorine 146.20 5.66 559.47 144.4 5.2 592.9 
Carbon dioxide 307.00 7.86 472.15 304.1 7.4 467.6 
Benzene 568.33 5.51 307.69 562.1 4.9 309.0 
Toluene 600.25 4.73 301.21 591.8 4.1 292.0 
 SAFT-VR MieExperimental data
SubstanceTc/KPc/MPaρc/kg m−3Tc/KPc /MPaρc /kg m−3
Methane 195.30 5.15 154.45 190.6 4.6 162.7 
Ethane 311.38 5.49 205.84 305.3 4.9 206.6 
Propane 376.20 4.77 219.98 369.8 4.2 220.5 
n-butane 432.68 4.27 228.04 425.1 3.8 228.0 
n-pentane 476.44 3.81 238.14 469.7 3.4 232.0 
n-hexane 515.29 3.44 241.16 507.8 3.0 233.2 
n-heptane 547.33 3.06 233.76 540.1 2.7 232.0 
n-octane 576.72 2.77 227.83 569.3 2.5 235.0 
n-nonane 602.20 2.52 224.89 594.6 2.3 232.1 
n-decane 626.37 2.31 219.85 617.7 2.1 233.3 
n-dodecane 668.75 1.99 214.26 658.2 1.8 221.4 
n-pentadecane 720.98 1.61 194.23 708.0 1.5 212.4 
n-eicosane 786.33 1.21 174.91 768.0 1.1 241.6 
Methanol 547.73 12.21 260.59 512.6 8.1 328.4 
Ethanol 554.41 8.83 247.96 514.0 6.3 276.4 
Propan-1-ol 560.43 6.89 256.34 536.9 5.2 275.2 
n-butan-1-ol 583.97 5.50 253.86 562.0 4.5 270.5 
Perfluoromethane 232.77 4.14 644.29 227.5 3.8 625.6 
Perluoroethane 295.46 3.24 634.96 293.0 3.0 613.3 
Perfluoropropane 347.88 2.79 648.71 345.0 2.6 628.0 
n-perfluorobutane 386.86 2.38 635.61 386.3 2.3 599.8 
n-perfluoropentane 421.36 2.13 634.72 420.6 2.0 609.5 
Fluorine 146.20 5.66 559.47 144.4 5.2 592.9 
Carbon dioxide 307.00 7.86 472.15 304.1 7.4 467.6 
Benzene 568.33 5.51 307.69 562.1 4.9 309.0 
Toluene 600.25 4.73 301.21 591.8 4.1 292.0 
FIG. 14.

Vapour-liquid equilibria for the homologous series of the n-alkanes from methane to n-decane: (a) temperature-density coexistence envelopes; (b) vapour-pressure curves; and (c) heats of vaporization ΔHv. Comparison of the experimental data129 (circles) with the description obtained using our new SAFT-VR Mie equation of state (continuous curves).

FIG. 14.

Vapour-liquid equilibria for the homologous series of the n-alkanes from methane to n-decane: (a) temperature-density coexistence envelopes; (b) vapour-pressure curves; and (c) heats of vaporization ΔHv. Comparison of the experimental data129 (circles) with the description obtained using our new SAFT-VR Mie equation of state (continuous curves).

Close modal

As a further critical assessment of the new theory we compare the representation of the VLE and thermophysical properties with the descriptions obtained with other commonly employed SAFT EOSs such as SAFT-VR SW, soft-SAFT, and PC-SAFT. The comparison with SAFT-VR SW is relevant since both equations are developed with a BH perturbation theory for the reference monomeric fluid. An improvement in the description can in a large part be attributed to the choice of the potential used to model the reference fluid, i.e., the more-realistic soft-core Mie potential rather than the hard-core square-well potential. The soft-SAFT EOS also affords a useful comparison as the theory is formulated using a semi-empirical EOS and RDF of the reference LJ fluid which are correlated to molecular-simulation data. An analysis of the relative performance of the SAFT-VR Mie and soft-SAFT EOSs is of interest since the deviations observed from the experimental data are therefore mainly a result of the choice of intermolecular potential model (a generalized LJ or a fixed LJ (12-6) form) and not to approximations used to evaluate the EOS and RDF of the LJ reference.

The AADs(%) of the theoretical values from the experimental data obtained for n-hexane, n-pentadecane, and n-eicosane for fluid-phase equilibria and selected thermodynamic properties are reported in Tables VII and VIII. The PC-SAFT, SAFT-VR SW and SAFT-VR LJC molecular parameters for these compounds are taken from previous work.56,65,71,72 The soft-SAFT parameters for n-eicosane are not available in the literature and are only reported as correlations in terms of the molar mass; furthermore, in order to present a fair comparison, parameters of the LJ chain segments considered in this approach are re-determined by estimation to experimental vapour pressures and saturated-liquid densities. It is important to recall that one should avoid using speed of sound data in the parameter estimation procedure with “traditional” SAFT approaches that incorporate a fixed model of the repulsive interactions. Indeed, it has been shown72 that one is unable to represent the speed of sound with these EOSs while retaining an accurate description of the vapour-liquid equilibrium. Our goal here is not to provide a comprehensive comparison for the entire n-alkane series but just to confirm the conclusions of the previous study72 with the new version of SAFT-VR Mie: the versatility of a potential model with a variable repulsive exponent is required to provide accurate simultaneous descriptions of the VLE and other thermodynamic properties, something which is not possible with a fixed LJ or SW form of interaction.

Table VII.

Vapour-liquid equilibria of selected n-alkanes. Comparison of the average absolute deviation AAD(%) of the theoretical description from the experimental data for the saturated liquid-density ρsat and the vapour pressure Psat obtained with the different SAFT equations of state. The definition of the AAD% is given in the caption of Table I.

 PC-SAFTSAFT-VR SWSoft-SAFTSAFT-VR Mie
SubstancePsatρsatPsatρsatPsatρsatPsatρsat
n-hexane 1.03 0.64 3.37 0.90 1.56 0.70 1.12 0.25 
n-pentadecane 2.26 0.77 9.37 0.98 1.96 2.90 1.76 0.60 
n-eicosane 4.38 0.62 5.33 0.47 3.56 3.06 1.86 0.96 
 PC-SAFTSAFT-VR SWSoft-SAFTSAFT-VR Mie
SubstancePsatρsatPsatρsatPsatρsatPsatρsat
n-hexane 1.03 0.64 3.37 0.90 1.56 0.70 1.12 0.25 
n-pentadecane 2.26 0.77 9.37 0.98 1.96 2.90 1.76 0.60 
n-eicosane 4.38 0.62 5.33 0.47 3.56 3.06 1.86 0.96 
Table VIII.

Second-derivative thermodynamic properties in the condensed-liquid phase of selected n-alkanes. Comparison of the average absolute deviation AAD(%) of the theoretical description from the experimental data for the speed of sound u and isobaric heat capacity Cp obtained with the different SAFT equations of state. The definition of the AAD(%) is given in the caption of Table I.

 PC-SAFTSAFT-VR SWsoft-SAFTSAFT-VR Mie
SubstanceuCpuCpuCpuCp
n-hexane 10.34 0.56 4.71 5.91 5.24 0.87 0.88 1.02 
n-pentadecane 16.00 0.76 12.60 8.60 10.27 3.36 0.58 1.55 
n-eicosane 16.90 2.55 15.90 10.00 11.23 3.74 1.33 0.74 
 PC-SAFTSAFT-VR SWsoft-SAFTSAFT-VR Mie
SubstanceuCpuCpuCpuCp
n-hexane 10.34 0.56 4.71 5.91 5.24 0.87 0.88 1.02 
n-pentadecane 16.00 0.76 12.60 8.60 10.27 3.36 0.58 1.55 
n-eicosane 16.90 2.55 15.90 10.00 11.23 3.74 1.33 0.74 

The second-derivative properties are closely interrelated through standard thermodynamic identities. For example, the isobaric and isochoric heat capacities, Cp and Cv, respectively, depend on the coefficient of thermal expansion αp = (1/V)(∂V/∂T)p and isothermal compressibility κT = −(1/V)(∂V/∂P)T through the relation CpCv = VT2T). We have already shown that the isobaric heat capacity, thermal expansivity, and isothermal compressibility are accurately represented with the SAFT-VR Mie EOS. In view of the thermodynamic interrelationships between these properties it follows that the isochoric heat capacity (a property that is often poorly represented with commonly employed EOSs) should also be accurately described with the theory. The temperature and pressure dependence predicted for Cv with the SAFT-VR Mie approach for selected hydrocarbons (n-hexane, n-octane, n-decane, and toluene) is depicted in Fig. 15, where the description is seen to be very satisfactory for the vapour and liquid branches, including the discontinuity at the vapour-liquid transition. We should point out that neither the heat capacity nor the thermal expansivity or isothermal compressibility are considered in the parameter estimation procedure so this represents a true indication of the predictive capability of our methodology.

FIG. 15.

The (a) temperature and (b) pressure dependence of the isochoric heat capacity Cv for selected hydrocabons: n-hexane (C6H14), n-octane (C8H18), n-decane (C10H22), and toluene (C7H8). Comparison of the experimental data129 (circles) with the description obtained using our new SAFT-VR Mie equation of state (continuous curves).

FIG. 15.

The (a) temperature and (b) pressure dependence of the isochoric heat capacity Cv for selected hydrocabons: n-hexane (C6H14), n-octane (C8H18), n-decane (C10H22), and toluene (C7H8). Comparison of the experimental data129 (circles) with the description obtained using our new SAFT-VR Mie equation of state (continuous curves).

Close modal

We now turn our attention to the n-perfluoroalkane homologous series. These fluorinated molecules are very challenging to describe with traditional intermolecular potential models due to their unique properties related partly to their high polarizability and repulsive nature of the fluorocarbon-fluorocarbon interactions.40,135 The molecular parameters for the first six n-perfluoroalkanes of the series are reported in Table III. In order to take into account the multipolar effects in an effective way for this class of substance, it is reasonable to relax the criterion on the value of the attractive exponent λa = 6, which can now be considered as an extra adjustable parameter; one should also recall that the segments forming the chains are coarse-grained representations of the molecules comprising a number of carbon and fluorine atoms so a simple London dispersion interaction is not necessarily the most-appropriate form. This is done in the case of perfluoromethane and perfluoroethane; for the sake of convenience the value of λa = 5.7506 obtained for perfluoroethane is transferred to higher members of the series. It is interesting to note that very repulsive Mie potentials (with repulsive exponents in the range λr ∼ 19 to 42) are required to correctly capture the VLE and the speed of sound of these compounds. This is a particular feature of the n-perfluoroalkanes which are characterized by very low values of the speed of sound in the liquid phase compared to their nonfluorinated n-alkane analogues. The theoretical description of the VLE and heat capacity of perfluoromethane with the PC-SAFT and SAFT-VR Mie EOSs are compared in Fig. 16. In the case of PC-SAFT, the model parameters with the number of segments fixed to ms = 1 (perfluoromethane is a near spherical molecule) have not previously been reported, and we have determined them in the usual way by estimation to the experimental saturated liquid densities and vapour pressures. The resulting values are σ = 4.1712 Å and ε/k = 194.51 K. It can be seen that one is not able to adequately capture the VLE of perfluoromethane with PC-SAFT (cf. Fig. 16(a)). The AAD(%) from the experimental data for the vapour pressure is 17.10% with PC-SAFT compared to 0.92% with SAFT-VR Mie. Three different isobars (10, 25, and 50 MPa) are also presented for the density and isobaric heat capacity in Figs. 16(a) and 16(c), respectively. A significantly improved description is achieved with the SAFT-VR Mie EOS. We also illustrate the excellent agreement between the experimental data and the SAFT-VR Mie approach for other perfluorinated compounds in Fig. 17. A key feature of our new theory is that it provides the capability of describing the VLE of pure substances even close to the critical region.

FIG. 16.

Selected thermodynamic properties of perfluoromethane CF4: (a) vapour-liquid temperature-density coexistence envelopes and density isobars of the condensed-liquid state; (b) vapour-pressure curves; and (c) the temperature dependence of the isobaric heat capacities Cp. The isobars in (a) and (c) are for pressures of P= 10, 25, and 50 MPa. Comparison of the experimental data (symbols) with the description obtained using the PC-SAFT (dashed curves) and our new SAFT-VR Mie equation of state (continuous curves).

FIG. 16.

Selected thermodynamic properties of perfluoromethane CF4: (a) vapour-liquid temperature-density coexistence envelopes and density isobars of the condensed-liquid state; (b) vapour-pressure curves; and (c) the temperature dependence of the isobaric heat capacities Cp. The isobars in (a) and (c) are for pressures of P= 10, 25, and 50 MPa. Comparison of the experimental data (symbols) with the description obtained using the PC-SAFT (dashed curves) and our new SAFT-VR Mie equation of state (continuous curves).

Close modal
FIG. 17.

The vapour-liquid equilibria for fluorine and the homologous series of the n-perfluoroalkanes from perfluoromethane to n-perfluoropentane: (a) temperature-density coexistence envelopes; (b) vapour-pressure curves; and (c) heats of vaporization ΔHv. Comparison of the experimental data129 (symbols) with the description obtained using our new SAFT-VR Mie equation of state (continuous curves).

FIG. 17.

The vapour-liquid equilibria for fluorine and the homologous series of the n-perfluoroalkanes from perfluoromethane to n-perfluoropentane: (a) temperature-density coexistence envelopes; (b) vapour-pressure curves; and (c) heats of vaporization ΔHv. Comparison of the experimental data129 (symbols) with the description obtained using our new SAFT-VR Mie equation of state (continuous curves).

Close modal

Associating substances have also been considered by assessing the adequacy of our approach for the first members of the n-alkanol series. In the current work we investigate a two-site association scheme corresponding to one donor and one acceptor site (denoted as 2B in the classification scheme of Huang and Radosz136). Seven molecular parameters need to be determined for each of these compounds: the number of segments ms, the diameter σ, the dispersion energy ε, the repulsive exponent λr, the attractive exponent λa, and the parameters associated with the site-site square-well potential including the association energy

$\epsilon ^{HB}_{ab}$
εabHB⁠, the distance to the centre of the Mie segment
$r^d_{ab}$
rabd
, and the range of interaction
$r^c_{ab}$
rabc
. The attractive range is kept at the London value of λa = 6 in order to reduce the number of adjustable parameters; the position of the site is also fixed to a constant value and set to
$r^d_{ab}/\sigma =0.4$
rabd/σ=0.4
, i.e., the sites are placed close to the surface of the Mie segment. This is the same site geometry that is used in the soft-SAFT approach and has proven to be very successful in modeling a wide variety of associating substances. The SAFT-VR Mie molecular parameters for methanol, ethanol, propan-1-ol, and n-butan-1-ol are reported in Table III. Once again, it can be seen (cf. Table IV) that one is able to accurately capture both the VLE and the second-derivative properties over a wide range of thermodynamic conditions with the SAFT-VR Mie EOS. A relatively soft Mie effective potential is found to be required to model these alkanols underlining the usefulness of developing an equation of state for a broad range of repulsive exponents, in accordance with previous findings.73 

A preliminary assessment of the adequacy of the new SAFT-VR Mie EOS in describing the vapour-liquid equilibria of binary mixtures is now presented. Mixtures of ethane + n-decane and carbon dioxide + n-decane are examined here, and the resulting description is compared with the existing studies based on a SAFT-VR SW description. In order to make a fair comparison of the two equations of state we first estimate the optimal unlike interaction dispersion energy parameters from the experimental VLE data. In the case of SAFT-VR SW, a deviation from the geometric mean of kij = 0.0204 is obtained for ethane + n-decane, and kij = 0.1203 for carbon dioxide + n-decane. The corresponding unlike energy parameters obtained with SAFT-VR Mie are kij = −0.0222 for ethane + n-decane and kij = 0.0500 for carbon dioxide + n-decane.

The resulting theoretical description with the SAFT-VR SW and SAFT-VR Mie EOSs are displayed in Fig. 18. While one is able to correlate the experimental data at low pressures with both approaches, it is striking to see the marked improvement obtained with the SAFT-VR Mie EOS at high pressures near the critical region. To better illustrate this point, we display the corresponding vapour-liquid critical lines of both mixtures in Fig. 19 together with the corresponding experimental data.137,138 The significant improvement seen near the critical region is not only a consequence of the choice of a more-realistic pair potential, but also to the fact that the new theory makes use of a perturbation expansion to higher-order than the SAFT-VR SW EOS providing a much better description particularly in the critical region. This is depicted in Figs. 18 and 19 where calculations obtained with the full SAFT-VR Mie perturbation expansion are compared with those obtained when only a first-order expansion of the RDF at contact is used to evaluate the chain term.

FIG. 18.

Isothermal pressure-composition slices of the vapour-phase equilibria for binary mixtures of: (a) ethane + n-decane at T = 444.15 K (circles) and T = 511.15 K (triangles); and (b) carbon dioxide + n-decane at T = 444.26 K (triangles) and T = 511.15 K (circles). Comparison of experimental data (symbols) with the corresponding description obtained with the SAFT-VR SW approach (dashed curves), our new SAFT-VR Mie equation of state (continuous curves), and the new SAFT-VR Mie with a chain term incorporating an expansion of the contact value of the radial distribution function truncated at first order (dotted curves).

FIG. 18.

Isothermal pressure-composition slices of the vapour-phase equilibria for binary mixtures of: (a) ethane + n-decane at T = 444.15 K (circles) and T = 511.15 K (triangles); and (b) carbon dioxide + n-decane at T = 444.26 K (triangles) and T = 511.15 K (circles). Comparison of experimental data (symbols) with the corresponding description obtained with the SAFT-VR SW approach (dashed curves), our new SAFT-VR Mie equation of state (continuous curves), and the new SAFT-VR Mie with a chain term incorporating an expansion of the contact value of the radial distribution function truncated at first order (dotted curves).

Close modal
FIG. 19.

Pressure-temperature projection of the vapour-liquid equilibria and critical lines for binary mixtures of: (a) ethane + n-decane; and (b) carbon dioxide + n-decane. Comparison of experimental data (symbols) with the corresponding description obtained with the SAFT-VR SW approach (dashed curves), our new SAFT-VR Mie equation of state (continuous curves), and the new SAFT-VR Mie with a chain term incorporating an expansion of the contact value of the radial distribution function truncated at first order (dotted curves).

FIG. 19.

Pressure-temperature projection of the vapour-liquid equilibria and critical lines for binary mixtures of: (a) ethane + n-decane; and (b) carbon dioxide + n-decane. Comparison of experimental data (symbols) with the corresponding description obtained with the SAFT-VR SW approach (dashed curves), our new SAFT-VR Mie equation of state (continuous curves), and the new SAFT-VR Mie with a chain term incorporating an expansion of the contact value of the radial distribution function truncated at first order (dotted curves).

Close modal

We have developed a new, accurate SAFT-VR Mie EOS based on a reference monomer system characterized by a Mie potential of variable repulsive and attractive range. A number of important modifications have been made in the theory for both the equation of state and radial distribution function of the reference fluid compared to previous work.72 In our current paper the Barker and Henderson84,85 high-temperature perturbation expansion is applied up to third order for the Helmholtz free energy of the reference Mie fluid, and up to second order for the corresponding RDF at contact. Analytical expressions for the second- and third-order contributions to the Helmholtz free energy are obtained by making use of Monte Carlo simulation data (to obtain a highly accurate fluctuation term) and VLE data up to the critical point (to determine the form of the third-order contribution). A key feature of the new SAFT-VR Mie EOS is the use of a third-order contribution which allows one to reproduce accurately the critical point of the Mie family of potentials. A more-rigorous derivation of the contact value of the RDF for Mie fluids is also presented. The high sensitivity of the Wertheim TPT1 description of the chain contribution to approximations that are made for the structure of the reference system is emphasized. The expression proposed here for the RDF is shown to capture accurately the non-conformal behaviour for chains of tangent Mie spheres for a broad range of values of the repulsive exponent. Good agreement is obtained with molecular simulation for chains of up to 8 segments; the description is found to become slightly less accurate for very long chains of up to 100 segments. This is not surprising however, and is mainly a consequence of the use of the linear approximation in the TPT1 treatment to evaluate the m-body distribution function at contact. A new parameterization of the mean-attractive energy of the hard-core Sutherland potential of variable range λ is developed which allows for the treatment of very short-ranged potentials with attractive/repulsive exponents spanning (5<λ<100). In view of the accuracy of the theory in the representation of the thermophysical properties of model systems, the adequacy of the SAFT-VR Mie EOS in describing the vapour-liquid equilibria and second-derivative properties of selected real substances (n-alkanes, fluorinated molecules, n-alcohols, and other compounds including carbon dioxide, benzene, and toluene) is investigated. In particular we demonstrate the versatility of varying the steepness of the Mie potential to represent correctly the softness/harness of the force field, e.g., the steep repulsive interactions between perfluorinated hydrocarbons. As Mie19 and others had recognized over a century earlier, the second-derivative thermodynamic properties (in Mie's case the compressibility) are found to be very sensitive to the precise form of the repulsive interactions. Our SAFT-VR Mie approach based on potentials of a Mie (rather than a fixed LJ) form therefore offers an enhanced capability in the simultaneous description of both VLE and derivative properties. It should also be noted that, in general, the new theory allows for a much better description of the critical region of non-associating molecules compared to previous work. We highlight this significant improvement by examining the fluid-phase behaviour of two prototypical binary mixtures which exhibit large regions of fluid-fluid criticality, namely, ethane + n-decane and carbon dioxide + n-decane. In both cases, the SAFT-VR Mie EOS provides a very satisfactory description of the entire phase envelope up to the critical point, despite the fact that no special treatment of the critical behaviour is considered here. An in-depth study of the performance of the theory for more-complex mixtures, including challenging associating substances, will be presented in future publications.

The main aims that we set ourselves in Sec. I have now been achieved: to develop a highly accurate analytical description of the thermodynamic and structural properties of associating chain molecules formed from Mie segments; to use the resulting SAFT-VR Mie EOS to describe accurately the properties of real substances and thus uncover the empirical form of the intermolecular potential which arises from subtle quantum-mechanical effects for a number of representative examples. The homonuclear formulation inherent in the SAFT-VR Mie approach is being extended to heteronuclear molecules (formed from segments of different types) to represent chemically distinct groups, which is the basis of the SAFT-γ Mie group contribution approach139 originally developed for heteronuclear molecules formed from SW segments.140,141 Our accurate SAFT-VR Mie approach is also proving invaluable in the development of coarse-grained intermolecular potential models for a wide variety of apolar and polar systems.126–128 The use of a versatile potential based on the Mie interaction of variable form offers a very significant advantage in the description of coarse-grained interactions. The repulsive and attractive exponents obtained for coarse-grained (multi-atomic) models turn out to be very different to those characterizing the smaller united-atom segments (traditionally represented with a LJ form), something that Maxwell13 had alluded to as early as 1867 in his analysis of the form of the repulsive interactions from studies of the transport properties of gases.

We are very grateful to Alejandro Gil-Villegas and Fernando del Río for useful input into aspects of the historical perspective, to Felix Llovell for his assistance with the soft-SAFT calculations, and Patrice Paricaud for helpful discussions regarding the fluctuation contribution. T.L. and C.A. thank the Engineering and Physical Sciences Research Council (EPSRC) of the UK for the award of postdoctoral fellowships. Additional funding to the Molecular Systems Engineering Group from the EPSRC (Grant Nos. GR/T17595, GR/N35991, EP/E016340, and EP/J014958), the Joint Research Equipment Initiative (JREI) (GR/M94426), and the Royal Society-Wolfson Foundation refurbishment scheme is also gratefully acknowledged.

1. Ideal contribution

The ideal contribution to the Helmholtz free energy is expressed as the usual sum over all species i in the mixture142 

\begin{equation}a^{\scriptsize\textit{IDEAL}}=\biggl (\sum _{i=1}^{n}x_{i}{\rm ln}\rho _{i}\Lambda _{i}^{3}\biggr )-1 ,\end{equation}
aIDEAL=i=1nxi ln ρiΛi31,
(A1)

where xi = Ni/N is the mole fraction, ρi = Ni/V is the molecular number density, Ni is the number of molecules,

$\Lambda _{i}^3$
Λi3 is the thermal de Broglie volume of species i, and n is the number of components in the mixture.

2. Monomer contribution for Mie fluid mixtures

The interactions between the monomers of species i and j in the mixture are modeled with the Mie potential

\begin{equation}u_{ij}^{Mie}=\mathcal {C}_{ij}\epsilon _{ij}\left( \left( \frac{\sigma _{ij}}{r_{ij}}\right)^{\lambda _{r,ij}} - \left( \frac{\sigma _{ij}}{r_{ij}}\right)^{\lambda _{a,ij}} \right) ,\end{equation}
uijMie=Cijεijσijrijλr,ijσijrijλa,ij,
(A2)

where

\begin{equation}\mathcal {C}_{ij}=\frac{\lambda _{r,ij}}{\lambda _{r,ij}-\lambda _{a,ij}}\left(\frac{\lambda _{r,ij}}{\lambda _{a,ij}} \right)^{\frac{\lambda _{a,ij}}{\lambda _{r,ij}-\lambda _{a,ij}}} .\end{equation}
Cij=λr,ijλr,ijλa,ijλr,ijλa,ijλa,ijλr,ijλa,ij.
(A3)

The monomer contribution to the Helmholtz free energy of the mixture can be expressed as70 

\begin{equation}a^{\scriptsize\textit{MONO}}= \biggl (\sum _{i=1}^{n}x_{i}m_{i}\biggr )a^{M},\end{equation}
aMONO=i=1nximiaM,
(A4)

where mi is the number of spherical segments making up the chains of component i. The expression for aM is obtained by applying the MX1b mixing rule70 to the expressions derived in Sec. III A. The MX1b prescription has been found to provide a consistent representation of both vapour-liquid and liquid-liquid critical behaviour.143 The monomer Helmholtz free energy aM is expressed as a third-order high-temperature expansion:

\begin{equation}a^{M}=a^{HS}+\beta a_1+(\beta )^{2}a_2+(\beta )^{3}a_3 .\end{equation}
aM=aHS+βa1+(β)2a2+(β)3a3.
(A5)

The expression of Boublík144 and Mansoori et al.145 for a multi-component mixture of hard spheres is used for the reference term,

\begin{eqnarray}a^{HS}&=&\frac{6}{\pi \rho _{s}}\Big [ \left(\frac{\zeta _{2}^{3}}{\zeta _{3}^{2}}-\zeta _{0} \right){\rm ln}(1-\zeta _{3})+ \frac{3\zeta _{1}\zeta _{2}}{(1-\zeta _{3})}\nonumber \\[4pt]&&+\, \frac{\zeta _{2}^{3}}{\zeta _{3}(1-\zeta _{3})^{2}} \Big ] ,\end{eqnarray}
aHS=6πρs[ζ23ζ32ζ0 ln (1ζ3)+3ζ1ζ2(1ζ3)+ζ23ζ3(1ζ3)2],
(A6)

where ρs is now the total number density of spherical segments, which is related to the total number density of the mixture ρ through

$\rho _{s}=\rho \sum _{i}^{n}x_{i}m_{i}$
ρs=ρinximi⁠. The moments of the number density ζl are defined as

\begin{equation}\zeta _{l}=\frac{\pi \rho _{s}}{6} \left(\sum _{i=1}^{n}x_{s,i}d_{ii}^{l} \right), \ \ l=0,1,2,3,\end{equation}
ζl=πρs6i=1nxs,idiil,l=0,1,2,3,
(A7)

where xs, i is the mole fraction of segments of component i,

\begin{equation}x_{s,i}=\frac{m_{i}x_{i}}{\sum _{k=1}^{n}m_{k}x_{k}} .\end{equation}
xs,i=mixik=1nmkxk.
(A8)

The effective diameter dii for the segments of each component is obtained from

\begin{equation}d_{ii}=\int _{0}^{\sigma _{ii}}\left( 1-{\rm exp}\left(-\beta u^{Mie}_{ii}(r) \right)\right) {\rm d}r.\end{equation}
dii=0σii1 exp βuiiMie(r)dr.
(A9)

The first-order mean-attractive energy a1 of the mixture is determined by summing the terms for each type of pair interaction,

\begin{equation}a_{1}=\sum _{i=1}^{n} \sum _{j=1}^{n} x_{s,i}x_{s,j}a_{1,ij},\end{equation}
a1=i=1nj=1nxs,ixs,ja1,ij,
(A10)

where the analytical expression for a1,ij is used,

\begin{eqnarray}\hspace*{-10pt}{a_{1,ij}}&=&{\mathcal {C}_{ij}\big [ x_{0,ij}^{\lambda _{a,ij}} \left( a_{1,ij}^{S}(\rho _s;\lambda _{a,ij}) + B_{ij}(\rho _s;\lambda _{a,ij}) \right)} \nonumber \\[4pt]&& -\, {x_{0,ij}^{\lambda _{r,ij}} \left( a_{1,ij}^{S}(\rho _s;\lambda _{r,ij}) +B_{ij}(\rho _s;\lambda _{r,ij}) \right) \big ]},\hspace*{10pt}\end{eqnarray}
a1,ij=Cij[x0,ijλa,ija1,ijS(ρs;λa,ij)+Bij(ρs;λa,ij)x0,ijλr,ija1,ijS(ρs;λr,ij)+Bij(ρs;λr,ij)],
(A11)

with x0,ij = σij/dij,

\begin{eqnarray}B_{ij}(\rho _s;\lambda _{ij}) &=& 2\pi \rho _s d_{ij}^3 \epsilon _{ij} \bigg [\! \frac{1-\zeta _{x}/2}{(1-\zeta _{x})^3} I_{\lambda ,ij}\nonumber\\&& -\, \frac{9\zeta _{x}(1+\zeta _{x})}{2(1-\zeta _{x})^3}J_{\lambda ,ij} \!\bigg ] ,\end{eqnarray}
Bij(ρs;λij)=2πρsdij3εij[1ζx/2(1ζx)3Iλ,ij9ζx(1+ζx)2(1ζx)3Jλ,ij],
(A12)

and

\begin{equation}\zeta _{x}=\frac{\pi \rho _{s}}{6}\sum _{i=1}^{n}\sum _{j=1}^{n}x_{s,i}x_{s,j}d_{ij}^{3} .\end{equation}
ζx=πρs6i=1nj=1nxs,ixs,jdij3.
(A13)

The quantities Iλ,ij and Jλ,ij are given by

\begin{equation}I_{\lambda ,ij}=-\frac{x_{0,ij}^{3-\lambda _{ij}}-1}{\lambda _{ij}-3},\end{equation}
Iλ,ij=x0,ij3λij1λij3,
(A14)
\begin{equation}J_{\lambda ,ij}=-\frac{ (x_{0,ij} )^{4-\lambda _{ij}}(\lambda _{ij}-3)-(x_{0,ij})^{3-\lambda _{ij}}(\lambda _{ij}-4)-1}{(\lambda _{ij}-3)(\lambda _{ij}-4)} ,\end{equation}
Jλ,ij=(x0,ij)4λij(λij3)(x0,ij)3λij(λij4)1(λij3)(λij4),
(A15)

and the perturbation terms

$a_{1,ij}^{S}$
a1,ijS for the corresponding Sutherland potentials are obtained from

\begin{equation}a_{1,ij}^{S}(\rho _s;\lambda _{ij})=-2\rho _s \left(\frac{\pi \epsilon _{ij} d_{ij}^3}{\lambda _{ij}-3} \right)\frac{1-\zeta _{x}^{eff}(\lambda _{ij})/2}{\bigl (1-\zeta _{x}^{eff}(\lambda _{ij})\bigr )^{3}} ,\end{equation}
a1,ijS(ρs;λij)=2ρsπεijdij3λij31ζxeff(λij)/21ζxeff(λij)3,
(A16)

where the effective packing fraction

$\zeta _{x}^{eff}(\lambda _{ij})$
ζxeff(λij) is obtained within a van der Waals one-fluid (VDW-1) approximation70 from the corresponding packing fraction of the mixture ζx:

\begin{eqnarray}\zeta _{x}^{eff}(\lambda _{ij})=c_{1}(\lambda _{ij})\zeta _{x}+c_{2}(\lambda _{ij})\zeta _{x}^{2}+c_{3}(\lambda _{ij})\zeta _{x}^{3}+c_{4}(\lambda _{ij})\zeta _{x}^{4} . \nonumber\!\!\!\!\! \\\end{eqnarray}
ζxeff(λij)=c1(λij)ζx+c2(λij)ζx2+c3(λij)ζx3+c4(λij)ζx4.
(A17)

The coefficients c1, c2, c3, and c4 are approximated by those of the pure fluid (cf. Eq. (41)),

\begin{equation}\left( \begin{array}{c}c_{1} \\[4pt]c_{2} \\[4pt]c_{3} \\[4pt]c_{4} \\[4pt]\end{array} \right) = \left( \begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0.81096 & 1.7888 & -37.578 & 92.284 \\[4pt]1.0205 & -19.341 & 151.26 & -463.50 \\[4pt]-1.9057 & 22.845 & -228.14 & 973.92 \\[4pt]1.0885 & -6.1962 & 106.98 & -677.64 \\[4pt]\end{array} \right) . \left( \begin{array}{c}1 \\[4pt]1/\lambda _{ij} \\[4pt]1/\lambda _{ij}^{2} \\[4pt]1/\lambda _{ij}^{3} \\[4pt]\end{array} \right) .\end{equation}
c1c2c3c4=0.810961.788837.57892.2841.020519.341151.26463.501.905722.845228.14973.921.08856.1962106.98677.64.11/λij1/λij21/λij3.
(A18)

The same treatment is applied to the second-order fluctuation term a2 of the mixture:

\begin{equation}a_{2}=\sum _{i=1}^{n} \sum _{j=1}^{n} x_{s,i}x_{s,j}a_{2,ij} ,\end{equation}
a2=i=1nj=1nxs,ixs,ja2,ij,
(A19)

where

\begin{eqnarray}a_{2,ij} &=& \frac{1}{2}K^{HS}(1+\chi _{ij})\epsilon _{ij} \mathcal {C}_{ij}^{2} \nonumber \\[4pt]&& \times\, \big [x_{0,ij}^{2\lambda _{a,ij}} \left( a_{1,ij}^{S}(\rho _s;2\lambda _{a,ij})+B_{ij}(\rho _s;2\lambda _{a,ij})\right) \nonumber \\[4pt]&&-\, 2x_{0,ij}^{\lambda _{a,ij}+\lambda _{r,ij}} \big( a_{1,ij}^{S}(\rho _s;\lambda _{a,ij}+\lambda _{r,ij})\nonumber\\&& +\, B_{ij}(\rho _s;\lambda _{a,ij}+\lambda _{r,ij})\big) \nonumber \\[4pt]&&+\, x_{0,ij}^{2\lambda _{r,ij}} \left( a_{1,ij}^{S}(\rho _s;2\lambda _{r,ij})+B_{ij}(\rho _s;2\lambda _{r,ij})\right) \big ] . \qquad\end{eqnarray}
a2,ij=12KHS(1+χij)εijCij2×[x0,ij2λa,ija1,ijS(ρs;2λa,ij)+Bij(ρs;2λa,ij)2x0,ijλa,ij+λr,ij(a1,ijS(ρs;λa,ij+λr,ij)+Bij(ρs;λa,ij+λr,ij))+x0,ij2λr,ija1,ijS(ρs;2λr,ij)+Bij(ρs;2λr,ij)].
(A20)

In this expression KHS is the isothermal compressibility of the mixture of hard spheres (cf. Eq. (16)),

\begin{equation}K^{HS}=\frac{\left(1-\zeta _{x} \right)^{4}}{1+4\zeta _{x}+4\zeta _{x}^{2}-4\zeta _{x}^{3}+\zeta _{x}^{4}} ,\end{equation}
KHS=1ζx41+4ζx+4ζx24ζx3+ζx4,
(A21)

and χij is given by

\begin{equation}\chi _{ij}=f_1(\alpha _{ij}) \bar{\zeta }_{x}+f_2(\alpha _{ij}){\bar{\zeta }_{x}}^5+f_3(\alpha _{ij}){\bar{\zeta }_{x}}^8 ,\end{equation}
χij=f1(αij)ζ¯x+f2(αij)ζ¯x5+f3(αij)ζ¯x8,
(A22)

with

\begin{equation}\bar{\zeta }_{x}=\frac{\pi \rho _{s}}{6}\sum _{i=1}^{n}\sum _{j=1}^{n}x_{s,i}x_{s,j}\sigma _{ij}^{3} ,\end{equation}
ζ¯x=πρs6i=1nj=1nxs,ixs,jσij3,
(A23)

and

\begin{equation}\alpha _{ij}=\mathcal {C}_{ij} \left(\frac{1}{\lambda _{a,ij}-3}-\frac{1}{\lambda _{r,ij}-3} \right) .\end{equation}
αij=Cij1λa,ij31λr,ij3.
(A24)

Finally, the third-order perturbation term is obtained from the following simple expression:

\begin{equation}a_{3,ij}=-\epsilon _{ij}^3f_4(\alpha _{ij})\bar{\zeta }_{x} {\rm exp} \big(f_5(\alpha _{ij}) \bar{\zeta }_{x} +f_6(\alpha _{ij} )\bar{\zeta }_{x}^2 \big) .\end{equation}
a3,ij=εij3f4(αij)ζ¯x exp f5(αij)ζ¯x+f6(αij)ζ¯x2.
(A25)

The functions fk (k = 1, …, 6) in Eqs. (A22) and (A25) are obtained with the generalized form of Eq. (20) using the coefficients defined in Table II:

\begin{equation}f_{k}(\alpha _{ij})=\left( \sum _{n=0}^{n=3} \phi _{k,n} \alpha _{ij}^n \right) \ \bigg / \left(1+\sum _{n=4}^{n=6} \phi _{k,n} \alpha _{ij}^{n-3}\right) .\end{equation}
fk(αij)=n=0n=3ϕk,nαijn/1+n=4n=6ϕk,nαijn3.
(A26)
3. Chain contribution for the formation of mixtures of chains of Mie segments

The Helmholtz free energy due to chain formation can be obtained as a sum of contributions from the RDF of the Mie fluid evaluated at contact,101 corresponding to tangentially bonded segments:

\begin{equation}a^{\scriptsize\textit{CHAIN}}=-\sum _{i=1}^{n}x_{i}\left(m_{i}-1 \right){\rm ln}g_{ii}^{Mie}(\sigma _{ii}) .\end{equation}
aCHAIN=i=1nximi1 ln giiMie(σii).
(A27)

An analytical expression for

$g^{Mie}_{ij}(\sigma _{ij})$
gijMie(σij) can be obtained by applying the MX1b70 mixing rule to the perturbation expansion, Eq. (72), described in Sec. III B, so that

\begin{eqnarray}g_{ij}^{Mie}(\sigma _{ij})&=&g^{HS}_{d,ij}(\sigma _{ij}){\rm exp}\big [ \beta \epsilon g_{1,ij}(\sigma _{ij})/g^{HS}_{d,ij}(\sigma _{ij}) \nonumber \\[4pt]&&+\, (\beta \epsilon )^2 g_{2,ij}(\sigma _{ij})/g^{HS}_{d,ij}(\sigma _{ij}) \big ] .\end{eqnarray}
gijMie(σij)=gd,ijHS(σij) exp [βεg1,ij(σij)/gd,ijHS(σij)+(βε)2g2,ij(σij)/gd,ijHS(σij)].
(A28)

The expression for

$g^{HS}_{d,ij}(\sigma _{ij})$
gd,ijHS(σij) is obtained from

\begin{equation}g^{HS}_{d,ij}(\sigma _{ij})={\rm exp} \big(k_{0}+k_{1}x_{0,ij}+k_{2}x_{0,ij}^{2}+k_{3}x_{0,ij}^{3}\big) ,\end{equation}
gd,ijHS(σij)= exp k0+k1x0,ij+k2x0,ij2+k3x0,ij3,
(A29)

where the density-dependent coefficients ki (for i = 0, …, 3) are calculated from the pure-component relations given in Eqs. (47)–(50) using a VDW-1 treatment:

\begin{equation}k_{0}=-\ln (1-\zeta _{x})+\frac{42\zeta _{x}-39\zeta _{x}^{2}+9\zeta _{x}^{3}-2\zeta _{x}^{4}}{6(1-\zeta _{x})^{3}} ,\end{equation}
k0=ln(1ζx)+42ζx39ζx2+9ζx32ζx46(1ζx)3,
(A30)
\begin{equation}k_{1}=\frac{(\zeta _{x}^{4}+6\zeta _{x}^{2}-12\zeta _{x})}{2(1-\zeta _{x})^{3}} ,\end{equation}
k1=(ζx4+6ζx212ζx)2(1ζx)3,
(A31)
\begin{equation}k_{2}=\frac{-3\zeta _{x}^{2}}{8(1-\zeta _{x})^{2}} ,\end{equation}
k2=3ζx28(1ζx)2,
(A32)
\begin{equation}k_{3}=\frac{\big(-\zeta _{x}^{4}+3\zeta _{x}^{2}+3\zeta _{x}\big)}{6(1-\zeta _{x})^{3}} .\end{equation}
k3=ζx4+3ζx2+3ζx6(1ζx)3.
(A33)

The first-order contribution to the RDF is obtained in terms of the mean-attractive energy of the corresponding Sutherland term by expressing Eq. (64) in its generalized VDW-1 form:

\begin{eqnarray}g_{1,ij}(\sigma _{ij})&=&\frac{1}{2\pi \epsilon _{ij} d_{ij}^{3}}\biggl [ 3\frac{\partial a_{1,ij}}{\partial \rho _{s}} \nonumber \\[4pt]&&-\, \mathcal {C}_{ij}\lambda _{a,ij}x_{0,ij}^{\lambda _{a,ij}} \frac{ a_{1,ij}^{s}(\rho _s;\lambda _{a,ij})+B_{ij}(\rho _s;\lambda _{a,ij}) }{\rho _{s}} \nonumber \\[4pt]&&+\, \mathcal {C}_{ij}\lambda _{r,ij}x_{0,ij}^{\lambda _{r,ij}} \frac{ a_{1,ij}^{s}(\rho _s;\lambda _{r,ij})+B_{ij}(\rho _s;\lambda _{r,ij}) }{\rho _{s}} \biggr ] . \nonumber\!\!\!\!\! \\\end{eqnarray}
g1,ij(σij)=12πεijdij3[3a1,ijρsCijλa,ijx0,ijλa,ija1,ijs(ρs;λa,ij)+Bij(ρs;λa,ij)ρs+Cijλr,ijx0,ijλr,ija1,ijs(ρs;λr,ij)+Bij(ρs;λr,ij)ρs].
(A34)

The second-order RDF is similarly obtained (cf. Eqs. (65) and (66)) in terms of the fluctuation of Sutherland potentials as

\begin{eqnarray}g_{2,ij}(\sigma _{ij})= (1+\gamma _{c,ij}){g_{2,ij}^{\scriptsize\textit{MCA}}}(\sigma _{ij}) ,\end{eqnarray}
g2,ij(σij)=(1+γc,ij)g2,ijMCA(σij),
(A35)

where

\begin{eqnarray}g_{2,ij}^{\scriptsize\textit{MCA}}(\sigma _{ij})&=& \frac{1}{2\pi {\epsilon _{ij}}^2 {d_{ij}}^{3}}\biggl [ 3\frac{\partial \frac{a_{2}}{1+\chi _{ij}}}{\partial \rho _s} \nonumber \\[4pt]&&-\, \epsilon K^{HS}{\mathcal {C}_{ij}}^2 \lambda _{r,ij}x_{0,ij}^{2\lambda _{r,ij}} \nonumber \\[4pt]&& \times\, \frac{ a_{1,ij}^{S}(\rho _s;2\lambda _{r,ij})+ B_{ij}(\rho _s;2\lambda _{r,ij}) }{\rho _s} \nonumber \\[4pt]&&+\, \epsilon _{ij} K^{HS} {\mathcal {C}_{ij}}^2(\lambda _{r,ij}+\lambda _{a,ij})x_{0,ij}^{\lambda _{r,ij}+\lambda _{a,ij}} \nonumber \\[4pt]&& \times\, \frac{ a_{1,ij}^{S}(\rho _s;\lambda _{r,ij}+\lambda _{a,ij})+ B_{ij}(\rho _s;\lambda _{r,ij}+\lambda _{a,ij}) }{\rho _s} \nonumber \\[4pt]&&- \, \epsilon _{ij} K^{HS} {\mathcal {C}_{ij}}^2 \lambda _{a,ij}x_{0,ij}^{2\lambda _{a,ij}}\nonumber \\[4pt]&& \times\, \frac{ a_{1,ij}^{S}(\rho _s;2\lambda _{a,ij})+ B_{ij}(\rho _s;2\lambda _{a,ij}) }{\rho _s} \biggr ] ,\end{eqnarray}
g2,ijMCA(σij)=12πεij2dij3[3a21+χijρsεKHSCij2λr,ijx0,ij2λr,ij×a1,ijS(ρs;2λr,ij)+Bij(ρs;2λr,ij)ρs+εijKHSCij2(λr,ij+λa,ij)x0,ijλr,ij+λa,ij×a1,ijS(ρs;λr,ij+λa,ij)+Bij(ρs;λr,ij+λa,ij)ρsεijKHSCij2λa,ijx0,ij2λa,ij×a1,ijS(ρs;2λa,ij)+Bij(ρs;2λa,ij)ρs],
(A36)

and the correction factor is obtained from the VDW-1 form of Eq. (63) as

\begin{eqnarray}\gamma _{c,ij}& =&\phi _{7,0} \big( {-}\tanh \big (\phi _{7,1}(\phi _{7,2}-\alpha _{ij})\big )+1 \big )\bar{\zeta }_x \theta _{ij}\nonumber \\[4pt]&& \times\, {\rm exp}\left(\phi _{7,3}\bar{\zeta }_x + \phi _{7,4}{\bar{\zeta }_x}^2 \right) ,\end{eqnarray}
γc,ij=ϕ7,0tanhϕ7,1(ϕ7,2αij)+1ζ¯xθij× exp ϕ7,3ζ¯x+ϕ7,4ζ¯x2,
(A37)

where θij = exp (βεij) − 1.

4. Association contribution for bonding between Mie segments in mixtures

Finally, the contribution to the Helmholtz free energy due to the association of the si site types on molecules of species i, with nai sites of type a on a particular molecule i is obtained from the Wertheim TPT1 formalism:101 

\begin{equation}a^{\it ASSOC}=\sum _{i=1}^{n}x_{i}\sum _{a=1}^{s_{i}}n_{ai}\left[ \ln X_{ai}-\frac{X_{ai}}{2} +\frac{1}{2} \right].\end{equation}
aASSOC=i=1nxia=1sinailnXaiXai2+12.
(A38)

The fraction Xai of component i that is not bonded at site a is obtained from the mass-action relations of the form

\begin{equation}X_{ai}=\frac{1}{1+\rho \sum _{j=1}^{n}x_{j}\sum _{b=1}^{s_j}n_{b,j} X_{bj}\Delta _{abij}} ,\end{equation}
Xai=11+ρj=1nxjb=1sjnb,jXbjΔabij,
(A39)

where Δabij characterizes the association strength between sites of type a on component i and sites of type b on component j. In the case of the Mie potential this association integral can be approximated as (cf. Sec. IV D)

\begin{equation}\Delta _{abij}=\sigma _{ij}^{3} F_{abij} I_{abij} ,\end{equation}
Δabij=σij3FabijIabij,
(A40)
\begin{equation}F_{abij}={\rm exp}\left(\beta \epsilon _{abij}^{HB} \right)-1,\end{equation}
Fabij= exp βεabijHB1,
(A41)

with

$\epsilon _{abij}^{HB}$
εabijHB being the site-site association energy, and Iabij a dimensionless quantity characterizing the density dependence and the volume available for bonding between two sites. In the case of mixtures we can write

\begin{equation}I_{abij}= g^{HS}_{d,ij}(d)K_{abij} ,\end{equation}
Iabij=gd,ijHS(d)Kabij,
(A42)

where

\begin{eqnarray}K_{abij}&=&4\pi {d_{ij}}^2\big [ \ln \big ( \big(r^c_{abij}+2r^d_{abij}\big)\big/d_{ij} \big ) \nonumber \\[4pt]&&\times \, \big(6{r^c_{abij}}^3+18 {r^c_{abij}}^2r^d_{abij}-24{r^d_{abij}}^3\big) \nonumber \\[4pt]&&+\, (r^c_{abij}+2r^d_{abij}-d_{ij}) \big(22{r^d_{abij}}^2-5r^c_{abij}r^d_{abij} \nonumber \\[4pt]&&-\, 7r^d_{abij}d_{ij} -8{r^c_{abij}}^2+r^c_{abij}d_{ij}+{d_{ij}}^2\big) \big ]\nonumber \\[4pt]&&\big/\, \big(72{r^d_{abij}}^2\sigma _{ij}^3\big),\end{eqnarray}
Kabij=4πdij2[lnrabijc+2rabijd/dij×6rabijc3+18rabijc2rabijd24rabijd3+(rabijc+2rabijddij)(22rabijd25rabijcrabijd7rabijddij8rabijc2+rabijcdij+dij2)]/72rabijd2σij3,
(A43)

and

$g^{HS}_{d,ij}(d)$
gd,ijHS(d) is evaluated with the Boublík144 expression as

\begin{eqnarray}g^{HS}_{d,ij}(d)&=&\frac{1}{1-\zeta _{3}}+3\frac{d_{ii}d_{jj}}{d_{ii}+d_{jj}}\frac{\zeta _2}{{(1-\zeta _3)}^2} \nonumber \\[4pt]&&+\, 2{\left (\frac{d_{ii}d_{jj}}{d_{ii}+d_{jj}}\right )}^2\frac{{\zeta _2}^2}{{(1-\zeta _3)}^3} .\end{eqnarray}
gd,ijHS(d)=11ζ3+3diidjjdii+djjζ2(1ζ3)2+2diidjjdii+djj2ζ22(1ζ3)3.
(A44)
5. Combining rules for the unlike intermolecular potential parameters in mixtures of associating Mie chains

A number of unlike intermolecular parameters also have to be defined in order to treat mixtures of the associating Mie chains. The unlike size parameters are obtained using the usual combining rule146 

\begin{equation}\sigma _{ij}=\frac{\sigma _{ii}+\sigma _{jj}}{2} .\end{equation}
σij=σii+σjj2.
(A45)

For computational efficiency, the same combining rule is used for the temperature dependent BH diameter dij,

\begin{equation}d_{ij}=\frac{d_{ii}+d_{jj}}{2} ,\end{equation}
dij=dii+djj2,
(A46)

instead of the more-rigorous expression,

\begin{equation}d_{ij}=\int _{0}^{\sigma _{ij}}\left( 1-{\rm exp}\left(\frac{-u^{Mie}_{ij}(r)}{kT} \right)\right) {\rm d}r .\end{equation}
dij=0σij1 exp uijMie(r)kTdr.
(A47)

The algebraic form of Eq. (A46) allows for a very similar representation of the thermodynamic properties of the mixture while keeping the numerical calculations to a minimum. Though the values of the other unlike interactions are commonly estimated from the experimental fluid-phase equilibria of mixtures, a number of additional combining rules can also be used in the case of limited or unavailable experimental data. To determine the unlike repulsive and attractive exponents, the following relationship can be used:

\begin{equation}\lambda _{k,ij}-3=\root \of {(\lambda _{k,ii}-3)(\lambda _{k,jj}-3)} , \ k=a,r.\end{equation}
λk,ij3=(λk,ii3)(λk,jj3),k=a,r.
(A48)

This is consistent with unlike Mie attractive interaction energies εij obtained using a Berthelot-like geometric criterion for the van der Waals attractive constants

$\alpha _{ij,\lambda _{a,ij}}^{VDW}$
αij,λa,ijVDW of a Sutherland λa, ij potential:

\begin{equation}\alpha _{ij,\lambda _{a,ij}}^{VDW} \; =\; \root \of{ \alpha _{ii,\lambda _{a,ii}}^{VDW}\alpha _{jj,\lambda _{a,jj}}^{VDW}} ,\end{equation}
αij,λa,ijVDW=αii,λa,iiVDWαjj,λa,jjVDW,
(A49)

where

\begin{equation}\alpha _{ij,\lambda _{a,ij}}^{VDW}=2\pi \epsilon _{ij}\int _{\sigma _{ij}}^{\infty }(\sigma _{ij}/r)^{\lambda _{a,ij}}r^{2} {\rm d}r .\end{equation}
αij,λa,ijVDW=2πεijσij(σij/r)λa,ijr2dr.
(A50)

In order to satisfy Eqs. (A48) and (A49), a specific geometric relation for εij has to be imposed:

\begin{equation}\epsilon _{ij}=\frac{\root \of {\sigma _{ii}^3\sigma _{jj}^3}}{\sigma _{ij}^3}\root \of {\epsilon _{ii}\epsilon _{jj}} .\end{equation}
εij=σii3σjj3σij3εiiεjj.
(A51)

In the case of the association contribution, combining rules can also be defined for the energy of the site-site interaction

$\epsilon ^{{{HB}}}_{abij}$
εabijHB⁠, and the parameters related to the geometry of the sites
$r^c_{abij}$
rabijc
and
$r^d_{abij}$
rabijd
. A Berthelot-like geometric mean can be used for the unlike site-site interaction energies,

\begin{equation}\epsilon ^{{{HB}}}_{abij}= \sqrt{\epsilon ^{{{HB}}}_{abii}\epsilon ^{{{HB}}}_{abjj} } ,\end{equation}
εabijHB=εabiiHBεabjjHB,
(A52)

and simple arithmetic combining rules are used for the unlike size parameters:

\begin{equation}r^c_{abij}=\frac{r^c_{abii}+r^c_{abjj}}{2} ,\end{equation}
rabijc=rabiic+rabjjc2,
(A53)

and

\begin{equation}r^d_{abij}=\frac{r^d_{abii}+r^d_{abjj}}{2} .\end{equation}
rabijd=rabiid+rabjjd2.
(A54)

Though our theory allows for a generic treatment of the association with different positions of the associating site

$r^d_{ab}$
rabd⁠, its value is commonly fixed in practice.

It should be emphasized that real substances often exhibit large deviations from the Lorentz-Berthelot combining rules, especially when the components of the mixtures are chemically very distinct.147 As we have already mentioned it is usually necessary to estimate the unlike parameters using experimental data. The unlike attractive dispersive energy εij is generally characterized by introducing a binary correction coefficient kij defined here from the relation

\begin{equation}\epsilon _{ij}=(1-k_{ij})\frac{\root \of {\sigma _{ii}^3\sigma _{jj}^3}}{\sigma _{ij}^3}\root \of {\epsilon _{ii}\epsilon _{jj}} .\end{equation}
εij=(1kij)σii3σjj3σij3εiiεjj.
(A55)

In a similar manner, a correction parameter γij for the combining rule for the unlike repulsive exponent can also be introduced as

\begin{equation}\lambda _{r,ij}-3=(1-\gamma _{ij})\,\, \root \of {(\lambda _{r,ii}-3)(\lambda _{r,jj}-3)} ,\end{equation}
λr,ij3=(1γij)(λr,ii3)(λr,jj3),
(A56)

while the attractive exponent is commonly taken as the London value, λa, ii = λa, jj = λa, ij = 6, or with the combining rule defined by Eq. (A48).

1.
J. S.
Rowlinson
,
Cohesion: A Scientific History of Intermolecular Forces
(
Cambridge University Press
,
Cambridge
,
2002
).
2.
I.
Newton
, Philosophiae Naturalis Principia Mathematica (London,
1687
) [English translation of 3rd ed. by A. Motte, The Mathematical Principles of Natural Philosophy (London, 1729)].
3.
R. J.
Bošković
, Theoria Philosophiae Naturalis, 2nd ed. (Venice,
1763
)
[English translation by
J. M.
Child
, A Theory of Natural Philosophy (Chicago,
1922
)].
4.
J.
Priestley
, The History and Present State of Electricity, with Original Experiments (London,
1767
).
5.
H.
Cavendish
,
Philos. Trans.
61
,
584
(
1771
).
6.
C.-A.
de Coulomb
, Histoire de l'Académie Royale des Sciences (
1784
), pp.
229
269
.
7.
T.
Young
,
Philos. Trans. R. Soc. London
95
,
65
(
1805
).
8.
P. S.
Laplace
,
Traité de Mécanique Céleste; Supplément au Dixième Livre, sur L'Action Capillaire
(
Courcier
,
Paris
,
1806
).
9.
O. F.
Mosotti
,
Sur les Forces qui Régissent la Constitution Intérieure des Corps, Aperçu pour Servir à la Détermination de la Cause des Lois de l' Action Moléculaire
(
Turin
,
1836
)
[English translation in Taylor's Scientific Memoirs (
1836
), Vol. 1, pp.
448
469
].
10.
H.
Yukawa
,
Proc. Phys. Math. Soc. Jpn.
17
,
48
(
1935
).
11.
J. S.
Rowlinson
,
Physica A
156
,
15
(
1989
).
12.
R.
Clausius
,
Ann. Phys. (Berlin)
176
,
353
(
1857
)
R.
Clausius
[
Philos. Mag.
14
,
108
(
1857
)].
13.
J. C.
Maxwell
,
Philos. Trans. R. Soc. London
157
,
49
(
1867
).
14.
L.
Boltzmann
,
Sitz. Akad. Naturwiss. Classe Kaiser Akad. Wissen Wien
66
,
275
(
1872
).
15.
J. D.
van der Waals
, “
Over de Continuiteit van den Gas-en Vloeistoftoestand
,” Ph.D. thesis (
University of Leiden
,
1873
);
J. D.
van der Waals
,
On the Continuity of the Gaseous and Liquid States
, edited by
J. S.
Rowlinson
(
North-Holland
,
Amsterdam
,
1988
).
16.
W.
Sutherland
,
Philos. Mag.
22
,
81
(
1886
);
W.
Sutherland
,
Philos. Mag.
24
,
113
(
1887
);
W.
Sutherland
,
Philos. Mag.
24
,
168
(
1887
).
17.
W.
Sutherland
,
Philos. Mag.
27
,
305
(
1889
).
18.
W.
Sutherland
,
Philos. Mag.
36
,
507
(
1893
).
19.
G.
Mie
,
Ann. Phys. (Berlin)
316
,
657
(
1903
).
20.
E. A.
Grüneisen
,
Z. Elektrochem. Angew. Phys. Chem.
17
,
737
(
1911
).
21.
E. A.
Grüneisen
,
Ann. Phys. (Berlin)
344
,
257
(
1912
).
22.
W. H.
Keesom
,
Commun. Phys. Lab. Univ. Leiden
32
(
Suppl. 24B
),
6
(
1912
).
23.
W. H.
Keesom
,
Proc. Sect. Sci. Konink. Akad. Weten. Amsterdam
15
,
643
(
1912
).
24.
J. E.
Jones
,
Proc. R. Soc. London, Ser. A
106
,
441
(
1924
).
25.
J. E.
Jones
,
Proc. R. Soc. London, Ser. A
106
,
463
(
1924
).
26.
J. E.
Lennard-Jones
,
Proc. R. Soc. London
43
,
461
(
1931
).
27.
J. C.
Slater
,
Phys. Rev.
32
,
349
(
1928
).
28.
J. C.
Slater
and
J. G.
Kirkwood
,
Phys. Rev.
37
,
682
(
1931
).
29.
S. C.
Wang
,
Phys. Z.
28
,
663
(
1927
).
30.
R.
Eisenschitz
and
F.
London
,
Z. Phys.
60
,
491
(
1930
).
31.
F.
London
,
Z. Phys. Chem. B
11
,
222
(
1930
).
32.
H.
Margenau
,
Phys. Rev.
38
,
747
(
1931
).
33.
H.
Margenau
,
Rev. Mod. Phys.
11
,
1
(
1939
).
34.
R. A.
Buckingham
,
Proc. R. Soc. London, Ser. A
168
,
264
(
1938
).
35.
R. A.
Buckingham
and
J.
Corner
,
Proc. R. Soc. London
189
,
118
(
1947
).
36.
A. J.
Stone
,
The Theory of Intermolecular Forces
,
International Series of Monographs on Chemistry
No.
32
(
Clarendon Press
,
Oxford
,
1996
).
37.
T.
Kihara
,
Nippon Sugaku-Buturegakukaisi
17
,
11
(
1943
).
38.
T.
Kihara
,
Rev. Mod. Phys.
25
,
831
(
1953
).
39.
J. O.
Hirschfelder
,
C. F.
Curtiss
, and
R. B.
Bird
,
Molecular Theory of Gases and Liquids
(
Wiley
,
New York
,
1954
).
40.
G. C.
Maitland
,
M.
Rigby
,
E. B.
Smith
, and
W. A.
Wakeham
,
Intermolecular Forces: Their Origin and Determination
,
International Series of Monographs on Chemistry No. 3
(
Clarendon Press
,
Oxford
,
1981
).
41.
W. G.
Chapman
,
K. E.
Gubbins
,
G.
Jackson
, and
M.
Radosz
,
Fluid Phase Equilib.
52
,
31
(
1989
).
42.
W. G.
Chapman
,
K. E.
Gubbins
,
G.
Jackson
, and
M.
Radosz
,
Ind. Eng. Chem. Res.
29
,
1709
(
1990
).
43.
M. S.
Wertheim
,
J. Stat. Phys.
35
,
19
(
1984
).
44.
M. S.
Wertheim
,
J. Stat. Phys.
35
,
35
(
1984
).
45.
M. S.
Wertheim
,
J. Stat. Phys.
42
,
459
(
1986
).
46.
M. S.
Wertheim
,
J. Stat. Phys.
42
,
477
(
1986
).
47.
M. S.
Wertheim
,
J. Chem. Phys.
85
,
2929
(
1986
).
48.
M. S.
Wertheim
,
J. Chem. Phys.
87
,
7323
(
1987
).
49.
E. A.
Müller
and
K. E.
Gubbins
, “
Associating fluids and fluid mixtures
,” in
Equations of State for Fluids and Fluid Mixtures
, edited by
J. V.
Sengers
,
R. F.
Kayser
,
C. J.
Peters
, and
H. J.
White
 Jr.
(
Elsevier
,
New York
,
2000
), Vol.
2
, p.
435
.
50.
E. A.
Müller
and
K. E.
Gubbins
,
Ind. Eng. Chem. Res.
40
,
2193
(
2001
).
51.
I. G.
Economou
,
Ind. Eng. Chem. Res.
41
,
953
(
2002
).
52.
P.
Paricaud
,
A.
Galindo
, and
G.
Jackson
,
Fluid Phase Equilib.
194
,
87
(
2002
).
53.
S. P.
Tan
,
H.
Adidharma
, and
M.
Radosz
,
Ind. Eng. Chem. Res.
47
,
8063
(
2008
).
54.
G.
Kontogeorgis
and
G. K.
Folas
,
Thermodynamic Models For Industrial Applications
(
Wiley-VCH
,
Weinheim
,
2010
).
55.
C.
McCabe
, and
A.
Galindo
, “
SAFT associating fluids and fluid mixtures
,” in
Applied Thermodynamics of Fluids
, edited by
A.
Goodwin
,
J. V.
Sengers
, and
C. J.
Peters
(
Royal Society of Chemistry
,
UK
,
2010
), Chap. XI.
56.
A.
Gil-Villegas
,
A.
Galindo
,
P. J.
Whitehead
,
S. J.
Mills
,
G.
Jackson
, and
A. N.
Burgess
,
J. Chem. Phys
106
,
4168
(
1997
).
57.
H.
Adidharma
and
M.
Radosz
,
Ind. Eng. Chem. Res.
37
,
4453
(
1998
).
58.
B. H.
Patel
,
H.
Docherty
,
S.
Varga
,
A.
Galindo
, and
G. C.
Maitland
,
Mol. Phys.
103
,
129
(
2005
).
59.
J.
Li
,
H.
He
,
C.
Peng
,
H.
Liu
, and
Y.
Hu
,
Fluid Phase Equilib.
276
,
57
(
2009
).
60.
H.
Guérin
,
J. Mol. Liq.
156
,
179
(
2010
).
61.
W. G.
Chapman
,
J. Chem. Phys.
93
,
4299
(
1990
).
62.
D.
Ghonasgi
and
W. G.
Chapman
,
AIChE J.
40
,
878
(
1994
).
63.
J. K.
Johnson
,
E. A.
Müller
, and
K. E.
Gubbins
,
J. Phys. Chem.
98
,
6413
(
1994
).
64.
T.
Kraska
and
K. E.
Gubbins
,
Ind. Eng. Chem. Res.
35
,
4727
(
1996
).
65.
L. A.
Davies
,
A.
Gil-Villegas
, and
G.
Jackson
,
Int. J. Thermophys.
19
,
675
(
1998
).
66.
F. J.
Blas
and
L. F.
Vega
,
Mol. Phys.
92
,
135
(
1997
).
67.
F. J.
Blas
and
L. F.
Vega
,
Ind. Eng. Chem. Res.
37
,
660
(
1998
).
68.
L. A.
Davies
,
A.
Gil-Villegas
, and
G.
Jackson
,
J. Chem. Phys.
111
,
8659
(
1999
).
69.
I.
Nezbeda
,
R.
Melnyk
, and
A.
Trokhymchuk
,
Fluid Phase Equilib.
309
,
174
(
2011
).
70.
A.
Galindo
,
L. A.
Davies
,
A.
Gil-Villegas
, and
G.
Jackson
,
Mol. Phys.
93
,
241
(
1998
).
71.
J.
Gross
and
G.
Sadowski
,
Ind. Eng. Chem. Res.
41
,
1084
(
2002
).
72.
T.
Lafitte
,
D.
Bessières
,
M. M.
Piñeiro
, and
J. L.
Daridon
,
J. Chem. Phys.
124
,
024509
(
2006
).
73.
T.
Lafitte
,
M. M.
Piñeiro
,
J. L.
Daridon
, and
D.
Bessières
,
J. Phys. Chem. B.
111
,
3447
(
2007
).
74.
A. J.
de Villiers
,
C. E.
Schwarz
,
A. J.
Burger
, and
G.
Kontogeorgis
,
Fluid Phase Equilib.
338
,
1
(
2013
).
75.
D. Y.
Peng
and
D. B.
Robinson
,
Ind. Eng. Chem. Fundam.
15
,
59
(
1976
).
76.
G.
Soave
,
Chem. Eng. Sci.
27
,
1197
(
1972
).
77.
G.
Galliero
,
T.
Lafitte
,
D.
Bessières
, and
C.
Boned
,
J. Chem. Phys.
127
,
184506
(
2007
).
78.
J. K.
Johnson
,
J. A.
Zollweg
, and
K. E.
Gubbins
,
Mol. Phys.
78
,
591
(
1993
).
79.
J.
Kolafa
and
I.
Nezbeda
,
Fluid Phase Equilib.
100
,
1
(
1994
).
80.
E. A.
Müller
,
K. E.
Gubbins
,
D. M.
Tsangaris
, and
J. J.
de Pablo
,
J. Chem. Phys.
103
,
3868
(
1995
).
81.
E. A.
Müller
,
L. F.
Vega
, and
K. E.
Gubbins
,
Mol. Phys.
83
,
1209
(
1994
).
82.
E. A.
Müller
,
L. F.
Vega
, and
K. E.
Gubbins
,
Int. J. Thermophys.
16
,
705
(
1995
).
83.
J. D.
Weeks
,
D.
Chandler
, and
H. C.
Anderson
,
J. Chem. Phys.
54
,
5237
(
1971
).
84.
J. A.
Barker
and
D.
Henderson
,
J. Chem. Phys.
47
,
4714
(
1967
).
85.
J. A.
Barker
and
D.
Henderson
,
Rev. Mod. Phys.
48
,
587
(
1976
).
86.
D. M.
Tsangaris
and
J. J.
de Pablo
,
J. Chem. Phys.
101
,
1477
(
1994
).
87.
P.
Paricaud
,
J. Chem. Phys.
124
,
154505
(
2006
).
88.
N. F.
Carnahan
and
K. E.
Starling
,
J. Chem. Phys.
51
,
635
(
1969
).
89.
S.
Zhou
,
J. Chem. Phys.
130
,
054103
(
2009
).
90.
A.
Malijevsky
and
S.
Labik
,
Mol. Phys.
60
,
663
(
1987
).
91.
F. F.
Betancourt-Cárdenas
,
L. A.
Galicia-Luna
, and
S. I.
Sandler
,
Fluid Phase Equilib.
264
,
174
(
2008
).
92.
W. R.
Smith
,
D.
Henderson
, and
J. A.
Barker
,
J. Chem. Phys.
53
,
508
(
1970
).
93.
B. J.
Zhang
,
Fluid Phase Equilib.
154
,
1
(
1999
).
94.
R.
Espindola-Heredia
,
F.
del Río
, and
A.
Malijevsky
,
J. Chem. Phys.
130
,
024509
(
2009
).
95.
H.
Okumura
and
F.
Yonezawa
,
J. Chem. Phys.
113
,
9162
(
2000
).
96.
J. J.
Potoff
and
D. A.
Bernard-Brunel
,
J. Phys. Chem. B
113
,
14725
(
2009
).
97.
I.
Nezbeda
and
G.
Iglesias-Silva
,
Mol. Phys.
69
,
767
(
1990
).
98.
T.
Boublík
,
Mol. Phys.
59
,
775
(
1986
).
99.
Y.
Tang
and
B. C.-Y.
Lu
,
J. Chem. Phys.
100
,
3079
(
1994
).
100.
F. W.
Tavares
,
J.
Chang
, and
S. I.
Sandler
,
Mol. Phys.
86
,
1451
(
1995
).
101.
W. G.
Chapman
,
G.
Jackson
, and
K. E.
Gubbins
,
Mol. Phys.
65
,
1057
(
1988
).
102.
G.
Jackson
,
W. G.
Chapman
, and
K. E.
Gubbins
,
Mol. Phys.
65
,
1
(
1988
).
103.
N. M. P.
Kakalis
,
A. I.
Kakhu
, and
C. C.
Pantelides
,
Ind. Eng. Chem. Res.
45
,
6056
(
2006
).
104.
Y.
Tang
and
B. C.-Y.
Lu
,
J. Chem. Phys.
99
,
9828
(
1993
).
105.
Y.
Tang
and
B. C.-Y.
Lu
,
J. Chem. Phys.
100
,
6665
(
1994
).
106.
E. A.
Müller
and
K. E.
Gubbins
,
Ind. Eng. Chem. Res.
34
,
3662
(
1995
).
107.
Y.
Tang
and
B. C.-Y.
Lu
,
Fluid Phase Equilib.
171
,
27
(
2000
).
108.
J. K.
Johnson
and
K. E.
Gubbins
,
Mol. Phys.
77
,
1033
(
1992
).
109.
C.
Domb
,
The Critical Point: A Historical Introduction to the Modern Theory of Critical Phenomena
(
Taylor and Francis
,
London
,
1996
).
110.
J. J.
Binney
,
N. J.
Dowrick
,
A. J.
Fisher
, and
M. E. J.
Newman
,
The Theory of Critical Phenomena: An Introduction to the Renormalization Group
(
Clarendon
,
Oxford
,
1992
).
111.
M. E.
Fisher
,
Rev. Mod. Phys.
70
,
653
(
1998
).
112.
M. A.
Anisimov
and
J. V.
Sengers
, in
Equations of State for Fluids and Fluid Mixtures
, edited by
J. V.
Sengers
,
R. F.
Kayser
,
C. J.
Peters
, and
H. J.
White
(
Elsevier
,
Amsterdam
,
2000
), Chap. XI.
113.
A.
Lofti
,
J.
Vrabec
, and
J.
Fischer
,
Mol. Phys.
76
,
1319
(
1992
).
114.
D.
Ben-Amotz
and
G.
Stell
,
J. Chem. Phys.
120
,
4844
(
2004
).
115.
P.
Orea
,
Y.
Reyes-Mercado
, and
Y.
Duda
,
Phys. Lett. A
372
,
7024
(
2008
).
116.
L. A.
Girifalco
,
J. Phys. Chem.
95
,
5370
(
1991
).
117.
L. A.
Girifalco
,
J. Phys. Chem.
96
,
858
(
1992
).
118.
L. A.
Girifalco
,
M.
Hodak
, and
R. S.
Lee
,
Phys. Rev. B
62
,
13104
(
2000
).
119.
F. A.
Escobedo
and
J. J.
de Pablo
,
Mol. Phys.
87
,
347
(
1996
).
120.
C.
Vega
,
C.
McBride
,
E.
de Miguel
,
F. J.
Blas
, and
A.
Galindo
,
J. Chem. Phys.
118
,
10696
(
2003
).
121.
L. G.
MacDowell
and
F. J.
Blas
,
J. Chem. Phys.
131
,
074705
(
2009
).
122.
P.
Paricaud
,
S.
Varga
, and
G.
Jackson
,
J. Chem. Phys.
118
,
8525
(
2003
).
123.
F. J.
Blas
and
L. F.
Vega
,
J. Chem. Phys.
115
,
4355
(
2001
).
124.
D.
Ghonasgi
and
W. G.
Chapman
,
J. Chem. Phys.
100
,
6633
(
1994
).
125.
J.
Chang
and
S. I.
Sandler
,
Chem. Eng. Sci.
49
,
2777
(
1994
).
126.
C.
Avendaño
,
T.
Lafitte
,
A.
Galindo
,
C. S.
Adjiman
,
G.
Jackson
, and
E. A.
Müller
,
J. Phys. Chem. B
115
,
11154
(
2011
).
127.
C.
Avendaño
,
T.
Lafitte
,
C. S.
Adjiman
,
A.
Galindo
,
E. A.
Müller
, and
G.
Jackson
,
J. Phys. Chem. B
117
,
2717
(
2013
).
128.
T.
Lafitte
,
C.
Avendaño
,
V.
Papaioannou
,
A.
Galindo
,
C. S.
Adjiman
,
G.
Jackson
, and
E. A.
Müller
,
Mol. Phys.
110
,
1189
(
2012
).
129.
P.
Linstrom
and
W.
Mallard
, “
NIST Chemistry Webbook, NIST Standard Reference Database Number 69
,” http://webbook.nist.gov, retrieved 10 December 2012.
130.
B. D.
Smith
and
R.
Srivastava
,
Thermodynamic Data for Pure Compounds: Part A Hydrocarbons and Ketones
(
Elsevier
,
Amsterdam
,
1986
).
131.
B. D.
Smith
and
R.
Srivastava
,
Thermodynamic Data for Pure Compounds: Part B Halogenated Hydrocarbons and Alcohols
(
Elsevier
,
Amsterdam
,
1986
).
132.
L.
Yelash
,
M.
Müller
,
W.
Paul
, and
K.
Binder
,
Phys. Chem. Chem. Phys.
7
,
3728
(
2005
).
133.
L.
Yelash
,
M.
Müller
,
W.
Paul
, and
K.
Binder
,
J. Chem. Phys.
123
,
014908
(
2005
).
134.
R.
Privat
,
R.
Gani
, and
J.-N.
Jaubert
,
Fluid Phase Equilib.
295
,
76
(
2010
).
135.
G. G.
Yee
,
J. L.
Fulton
, and
R. D.
Smith
,
J. Phys. Chem.
96
,
6172
(
1992
).
136.
S. H.
Huang
and
M.
Radosz
,
Ind. Eng. Chem. Res.
29
,
2284
(
1990
).
137.
B. A.
Bufkin
,
R. L.
Robinson
 Jr.
,
S. S.
Estrera
, and
D. K.
Luks
,
J. Chem. Eng. Data
31
,
421
(
1986
).
138.
H. H.
Reamer
and
B. H.
Sage
,
J. Chem. Eng. Data
8
,
508
(
1963
).
139.
V.
Papaioannou
,
T.
Lafitte
,
C.
Avendaño
,
C. S.
Adjiman
,
G.
Jackson
,
E. A.
Müller
, and
A.
Galindo
, “
Group contribution methodology based on the statistical associating fluid theory for heteronuclear molecules formed from Mie segments
,”
J. Chem. Phys.
(submitted).
140.
A.
Lymperiadis
,
C. S.
Adjiman
,
A.
Galindo
, and
G.
Jackson
,
J. Chem. Phys.
127
,
234903
(
2007
).
141.
A.
Lymperiadis
,
C. S.
Adjiman
,
G.
Jackson
, and
A.
Galindo
,
Fluid Phase Equilib.
274
,
85
(
2008
).
142.
L. L.
Lee
,
Molecular Thermodynamics of Nonideal Fluids
(
Butterworth
,
Boston
,
1988
).
143.
C.
McCabe
,
A.
Galindo
,
A.
Gil-Villegas
, and
G.
Jackson
,
J. Phys. Chem. B
102
,
8060
(
1998
).
144.
T.
Boublík
,
J. Chem. Phys.
53
,
471
(
1970
).
145.
G. A.
Mansoori
,
N. F.
Carnahan
,
K. E.
Starling
, and
T. W.
Leland
,
J. Chem. Phys.
54
,
1523
(
1971
).
146.
J. S.
Rowlinson
and
F. L.
Swinton
,
Liquids and Liquid Mixtures
(
Butterworth
,
London
,
1982
).
147.
A. J.
Haslam
,
A.
Galindo
, and
G.
Jackson
,
Fluid Phase Equilib.
266
,
105
(
2008
).
148.
Y.-J.
Sheng
,
A. Z.
Panagiotopoulos
,
S. K.
Kumar
, and
I.
Szleifer
,
Macromolecules
27
,
400
(
1994
).