Statistical fluid theory for systems of variable range interacting via triangular-well pair potential

Journal of Molecular Liquids 265 (2018) 337–346

Contents lists available at ScienceDirect

Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq

Statistical fluid theory for systems of variable range interacting via triangular-well pair potential Víctor M. Trejos a , Alejandro Martínez b , Néstor E. Valadez-Pérez c, * a b c

Instituto de Química, Universidad Nacional Autónoma de México, Ciudad de México Apdo. Postal 70213, Coyoacán, 04510, Mexico División de Ciencias e Ingenierías, Campus León, Universidad de Guanajuato, Loma del Bosque 103, Lomas del Campestre, León 37150, Guanajuato, Mexico Department of Chemical and Biological Engineering, Illinois Institute of Technology, 3440 S. Dearborn St., Chicago, IL 60616, USA

A R T I C L E

I N F O

Article history: Received 2 April 2018 Accepted 25 May 2018 Available online 31 May 2018 Keywords: Barker-Henderson Equation of state Triangular-well fluids

A B S T R A C T The interest for describing the thermodynamic properties of fluids is growing in many fields from pure theoretical to those with technological applications. In this work, we have developed an analytical expression for the Helmholtz free energy of particles interacting via the triangular-well pair potential (TW), in the Barker-Henderson framework. This theoretical equation of state (EoS) uses the analytical radial distribution function of the hard-sphere (HS) fluid coupled with the Barker-Henderson perturbation theory. Thermodynamic properties as vapor-liquid coexistence, vapor pressures, internal energies and critical properties are obtained for attractive particles whose interaction potential range lies between 1.2 ≤ k ≤ 2.6. The performance of the TW equation of state is assessed by comparing with Monte Carlo (MC) simulation data. Finally, this theoretical approach is used to model real fluids as methane, oxygen, fluoromethane and hydrogen sulfide, in a wide range of temperature and pressures. © 2018 Elsevier B.V. All rights reserved.

1. Introduction In the context of the liquid theory, the hard sphere system is the simplest representation of particles in a real fluid at high temperature. The existence of an attraction among particles lets the HS system to present other phases, additional to the fluid and crystal ones. The attractive interaction makes possible the vapor-liquid coexistence or even a crystal-crystal coexistence [1]. The square-well (SW) pair potential is the simplest attractive potential, the thermodynamic properties and phase diagram of SW fluids have been extensively studied in the literature, through theoretical approaches [1–19] and simulations [20–24]. Even though the SW potential is an unphysical interaction, it has proved to be very useful due to its simple form and because other kind of potentials, with a more complex functional form, as the Yukawa, Lennard-Jones, AsakuraOosawa or Morse potentials, can be represented by a suitable SW according to the extended law of corresponding states [25]. Besides, the SW potential could be a crude representation of the effective interaction in complex fluids, as colloids and protein suspensions. In those systems the attraction range is smaller that the molecular diameter [19,26]. But if the attraction range is longer than the

* Corresponding author. E-mail address: [email protected] (N.E. Valadez-Pérez).

https://doi.org/10.1016/j.molliq.2018.05.116 0167-7322/© 2018 Elsevier B.V. All rights reserved.

molecular diameter, the system can be a representation of a simple fluid composed by neutral molecules. Another attractive potential, with a slightly more complicated mathematical form than the SW, is the TW potential. A schematic representation of the TW potential with an attraction range of k = 1.5 is represented in Fig. 1. Despite its simple form, it has not received the same attention as the SW potential. Some of earlier work about TW fluids has been focused on developing equations of state to predict the thermodynamic properties as well as the phase behavior, for different attraction ranges, and to compare the predictions with simulation data. Largo and Solana [27] obtained an analytical expression for the first and second order perturbation terms of the Barker-Henderson (BH) perturbation theory (PT) [2,28] based on the Chang-Sandler expression for the radial distribution function (RDF) [29]. The BH [2,28] theory has been successfully applied to model the thermodynamic properties of other kind of fluids, as for example square well (SW) and Lennard-Jones (LJ) fluids. In a more recently work, they analyzed the performance of several approximations for the second-order perturbation term [30]. BetancourtCárdenas et al. [31] developed an EoS based on the BH PT to calculate the properties of TW fluids and solids with a short-ranged potentials, k = 1.1, to intermediate-range ones, k = 2.0 and in the Ref. [32] they calculated thermodynamic properties as well as the phase diagram of TW fluids through perturbation theory, self-consistent Orstein-Zernike approach and Monte Carlo simulation. They focused

338

V. Trejos et al. / Journal of Molecular Liquids 265 (2018) 337–346

scheme and some details about our implementation to calculate the free energy and the vapor-liquid equilibrium (VLE). In Section 4, the results of the TW-EoS against MC molecular simulation results are presented for different thermodynamic properties as VLE, vapor pressure, critical points and internal energy. Additionally, we have included the VLE and vapor pressure calculation of fluoromethane (R41), hydrogen sulfide (H2 S), methane (CH4 ) and oxygen (O2 ) by using the TW-EoS, the molecular parameters were obtained by fitting our results to the data from Ref. [45]. Finally, in Section 5, some conclusions and perspectives are provided. 2. Theoretical framework 2.1. Perturbation theory for TW fluids In this work, the description of the fluids is given in the framework of the BH PT [2,28]. We consider a system of molecules interacting through the potential, u(r), which can be split in two contributions Fig. 1. a) Schematic representation of the triangular-well pair potential model for k = 1.5. b) Particles separated by a distance smaller than s exert an infinite repulsion, but if the distance is less than k they attract each other, otherwise they have not interaction. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

u(r) = uHS (r) + uatt (r),

in particles interacting through a long-ranged potential, 2 ≤ k ≤ 4. Adhikari and Kofke [33] obtained the vapor-liquid, fluid-solid and solid-solid coexistence boundaries using MC simulations for several attraction ranges, 1.05 ≤ k ≤ 2.5. Zhou [34] used his own thermodynamic perturbation theory to obtain a fifth-order EoS valid for ranges, 1.5 ≤ k ≤ 3. Rivera et al. [35] obtained an EoS based on perturbation theory using the macroscopic compressibility approximation, they obtained accurate results for ranges around k = 2.0. Bárcenas et al. [36] obtained the binodal and interfacial properties for TW with 1.5 ≤ k ≤ 3.0 and in Ref. [37] they used MC simulations with the replica exchange technique to obtain the binodal of short-ranged TW fluids, they considered some values in the range: 1.15 ≤ k ≤ 1.4. Sengupta et al. [38,39] studied the properties of TW fluids in confinement. Guérin [40,41] developed a framework to predict the thermodynamic behavior of TW fluids, also useful to study other kind of fluids as the SW and Sutherland ones. The TW potential has been used to model experimental systems in a pair of works, those are the following. In Ref. [42], it is shown that particles in a microemulsion can be modeled with a short-ranged attractive TW with k = 1.47. In Ref. [31], the EoS developed is used to model methane and ethane, the range of the well of these fluids are 2.4925 and 2.1834, respectively. The previous works serve as our motivation to develop an EoS valid in the colloidal regime, for short-ranged potentials, as well as in the simple fluids regime, which correspond to long-ranged attractions. The main purpose of this work is to present an analytic EoS for the description of thermodynamic properties of fluids in the framework of the Barker-Henderson perturbation theory [2,28] and compare its predictions with recent MC simulation data reported the literature [31,32,37]. In this framework, the Helmholtz free energy, A, is expressed as the sum of a reference term, obtained from a reference potential, and a perturbation. The perturbation contribution is obtained by following a standard procedure, in particular we follow the one developed for SW fluids [43,44]. This paper is organized as follows. In Section 2, the fundamentals of the BH PT [2,28] for TW pair potentials are revised. Additionally, an analytic parameterization for the first and the second order perturbation theory terms of the so-obtained EoS for TW potentials, is presented. Section 3 presents the Monte Carlo (MC) simulation

uHS (r) =

(1)

where uHS (r) and uatt (r) are the repulsive HS potential and the attractive pair potential (which plays the role of the perturbation), respectively. The HS contribution is given by 

∞, if r < s, 0, if r ≥ s,

(2)

where r is the distance between the centers of two particles and s is the particle diameter. For a system of particles interacting via TW pair potential, the attractive contribution is given by

u

attr

⎧ if r < s, ⎨ 0, (r) = −e(k − r)/(k − s), if s ≤ r < k, ⎩ 0, if r ≥ k,

(3)

where e and k are the energy-depth and range of the potential, respectively. Once the pair potential, u(r), is split in two contributions the Helmholtz free energy, A, for a system of interacting particles can be written in the BH [2,28] approach as A Aid Aexc = + , NkT NkT NkT

(4)

where k is the Boltzmann’s constant, N is the number of particles in the system, T is the absolute temperature and Aid and Aexc are the ideal and excess Helmholtz free energy. The ideal gas contribution to the free energy is given by [46]   Aid = ln qk3B − 1, NkT

(5)

where kB the de Broglie wavelength and q the density of the system. The excess energy is composed by different parts to take into account the different interactions among particles, then there is a contribution for the excluded volume, a second contribution due to the dispersion forces and another one that comes from the hydrogen bonds, etc. The excess free energy could be written as Aexc AHS Aattr Aassoc = + + + ··· NkT NkT NkT NkT

(6)

to take into account the contribution of the reference, the perturbation and also the associating potentials, respectively. The free

V. Trejos et al. / Journal of Molecular Liquids 265 (2018) 337–346

339

energy of a HS fluid is well described by the Carnahan-Starling expression [47], given by

where Z is the compressibility factor. The internal energy U is obtained from the Helmholtz free energy as follows

AHS 4g − 3g2 = , NkT (1 − g)2

U ∂ = −T NkT ∂T

(7)

where g is the packing fraction, g = pqs 3 /6, being q the number density. The attractive Helmholtz free energy, Aattr , in hightemperature approach could be represented as a power expansion in the inverse of the temperature, A1 Aattr A2 =b + b2 + ··· NkT NkT NkT

(8)

where b = 1/kT and the coefficients Ai are the i = 1, 2, . . . , norder perturbation terms of the BH PT [2,28]. Here, we consider a second order expansion and neglect terms of O(b3 ). The A1 and A2 terms are the main contributions to the free energy in the range of temperatures we are interested and also because there are difficulties to precisely calculate higher order terms, as was pointed out by Espíndola et al. [18]. The first-order contribution to the Helmholtz free energy, A1 , can be written as A1 = 2pq NkT



∞

gHS (r)uatt (r)r2 dr,

(9)

s

where gHS (r) is the radial distribution function (RDF) of the HS fluid. In the so-called “local-compressibility approximation” (LCA) the second-order perturbation term is given by,  ∞

A2 ∂ q = −pK HS q gHS (r)[uatt (r)]2 r2 dr , NkT ∂q s



A NkT



.

(14)

V,N

The thermodynamic equilibrium between liquid (l) and vapor (v) phases at a given temperature for any interaction range is obtained by solving the chemical (l l = l v ), mechanical (P l = Pv ) and thermal equilibrium (T l = Tv ) for each phase. All the points that satisfy those equations form the so-called binodal line in the q vs. T plane. 2.2. Parameterization and effective packing fraction Here, we detail the calculation of the perturbation terms, A1 and A2 , by following the methodology employed for the case of SW fluids [43,44]. They together with the reference term, AHS , will be refereed as the TW-EoS. The first-order perturbation term, A1 , is expressed as A1 = 12g NkT



k

gHS (x)uatt (x)x2 dx,

1



k

= 12ggHS (n; g)

uatt (x)x2 dx,

1

  = −4ggHS (n; g) k3 + k2 + k − 3 ,

where g is the packing fraction, x = r/s is the reduced distance and n is a distance that satisfied the mean-value theorem, details can be founded in Refs. [43,44]. As it is discussed in Ref. [43], the RDF of the HS system is replaced by the term gHS (n; g), which comes from obtaining g(n) at a given effective packing fraction, geff , such that gHS (n; g) = gHS (1; geff ),

(10)

(15)

(16)

where gHS (1; geff ) is the RDF at the contact value, which is approximated by the Carnahan-Starling expression [47],

HS

where KHS = kT(∂ q/∂ P )T is the isothermal compressibility of the HS fluid, which is given by

K HS =

(1 − g)4 . 1 + 4g + 4g2

(11)

according with the Perkus-Yevick expression [48]. Eqs. (9)–(10) require the knowledge of a suitable expression for the RDF of the HS fluid for any value of g and k. Therefore, to have an analytical approximation for the gHS (r), for arbitrary values of g and k, is needed in order to determine the Helmholtz free energy for A1 and A2 . Here, we use the approximation for gHS (r) proposed in Ref. [29] to obtain A1 and A2 . From the numerical solution of Eqs. (9)– (10), we proposed a packing fraction parameterization in order to show an analytic EoS for TW fluids in a wide range of values of g and k. Once the analytic expression for the Helmholtz free energy is known as function of the packing fraction, thermodynamic variables as the pressure (P) and chemical potential (l), can be obtained by the standard thermodynamic relations  P=−

∂A ∂V



and

l=

T,N

∂A ∂N

,

(12)

T,V

and Z=

1 PV = NkT kT



∂A ∂N

− T,V

A , NkT

(13)

gHS (1; geff ) =

(1 − geff /2) . (1 − geff )3

(17)

The second order perturbation term A2 in the LCA is written as

  k  2 A2 1 ∂ 12g = − K HS g gHS (x) uatt (x) x2 dx NkT 2 ∂g 1

∂ 1 HS 4 A1 g − Ir , = K k 2 k − 1 ∂g NkT

(18)

where the second term in the right hand side is given by 

k

Ir = 12g

gHS (x)uatt (x)x3 dx,

(19)

1

and the function Ir can be approximated by using the first-order perturbation term A1 , as follows Ir = f (k)A1 /NkT,

(20)

being f(k) an arbitrary function of k. Here we use f(k) = 7k/11, this function gives us a good approximation to Ir in a wide range of values for the density and k. The comparison between this approximation and the numerical integration is presented in Appendix A. The A2 term is given in a compact form as 2 k4 A2 ∂ = K HS g NkT 11 (k − 1) ∂g



A1 NkT



.

(21)

340

V. Trejos et al. / Journal of Molecular Liquids 265 (2018) 337–346

The parameterization for geff is the missing part to determine the TW-EoS. Here, we follow the approach used by Gil-Villegas et al. [43], in which the effective packing fraction geff is a function of g and k. We have obtained values for A1 /NkT solving numerically Eq. (9) by using RDF for HS given by Chang et al. [29], along the range of 1.1 ≤ k ≤ 2.6 and 0.0001 ≤ g ≤ 0.42. The obtained parameterization for geff is given by geff = c1 g + c2 g2 + c3 g3 ,

4.1. Free energy

⎛

(23)

The parameterization given by Eqs. (22)–(23) coupled with Eqs. (15)– (18) provides analytical expressions for A1 and A2 in order to predict thermodynamic properties of the TW fluids. 3. Monte Carlo simulation We considered HS particles of diameter s interacting via a TW pair-potential, characterized by the variable potential range k and energy depth 4 (see Fig. 1(b)). In that figure, the HS diameter is represented in blue color and the pair-potential interactions are represented in red color. These kinds of models have been used to predict coexistence and interfacial properties of noble gases [49], methane and ethane [32]. Monte Carlo (MC) simulations are performed to calculate the free energy and the phase equilibrium of the TW fluid for different attraction ranges, k. The A1 and A2 terms in Eqs. (9)–(10) are computed directly from MC simulations of HS particles in the canonical ensemble. The perturbation terms are given by [2,18,31]   U0 A1 (q∗ , k) = NkT N     (U − U 0 )2 0 A2 (q∗ , k) =− NkT 2N

For the sake of simplicity, the following notation for reduced units is used throughout the paper: packing fraction g = (p/6)qs 3 , temperature T∗ = kT/4, reduced density q∗ = qs 3 and pressure P∗ = Ps 3 /4. The simulation data results and experimental data were displayed with filled and empty symbols whereas the theoretical TW-EoS results were represented by solid and dashed lines.

(22)

where ci , with i = 1, 2, 3 coefficients are given by ⎞ ⎛ ⎞ ⎞⎛ c1 1.94785 −1.03659 0.14141 1 ⎝ c2 ⎠ = ⎝ 2.65578 −3.37315 0.85580 ⎠ ⎝ k ⎠ . −2.83219 2.92387 −0.66341 c3 k2

4. Results and discussion

In Fig. 2(a), the A1 curves obtained by TW-EoS and MC simulations for several values of the attraction ranges 1.2 ≤ k ≤ 2.5 are shown. From short, k = 1.2, to long-range attractions k = 2.5, the A1 terms of TW-EoS results are in concordance with the simulation data. In all the cases the error bar of the simulation results was smaller than the size of the symbols. As is expected, A1 is a monotonous decreasing function of q∗ and depends strongly on the value of k. Fig. 2(b) shows the results of the term A2 at different values of k. The A2 follows the same trend as the simulation results but the prediction of the theoretical expression overestimates its value in the whole range of densities. This difference is close to 200%, which seems to be huge but the prediction is closer to the simulation than the reported in previous theoretical studies for the same system, see for example Ref. [31]. Our approach seems to be more adequate than the normal perturbation theory using the macroscopic compressibility expression [31] up to second order, however our expression for the A2 term is not that accurate as in the case of the fifth-order EoS presented by Zhou [34].

(24)

(25)

where 0 indicates the average over the reference system and U is the virtual total energy of interaction, i.e., we consider that particles interact through the potential in Eq. (3), even though the actual interaction is only the HS potential. To calculate the vapor-liquid equilibrium we performed the usual Gibbs ensemble MC simulation (GEMC) [50] for a system of N = 1000 particles initially placed in two boxes with N1 = N2 = 500 particles. Details of the implementation can be found in Ref. [51]. The simulation ran for 105 cycles to reach the equilibrium and for an equal number of cycles to compute the densities in both coexisting phases, vapor and liquid reduced densities q∗v and q∗l . Since we are not dealing with very short-ranged attractions and the interaction is very simple, the Gibbs ensemble technique still a good technique to develop our study [19]. For the estimation of the critical point, we used the scaling law and the rectilinear diameters recipe as a first estimate for the Wegner expansion [52], i.e., the reduced critical parameters q∗ c and T∗ c are given by q± = q∗ c + C2 |t| ±

1 Bo |t|0.325 . 2

(26)

where q− and q+ denote the vapor and liquid densities, respectively, and t = 1 − T∗ /T∗ c . This expression allowed us to fit the coexistence data and obtain estimates for q∗ c , T∗ c and the amplitude terms C2 and Bo .

Fig. 2. Free energy perturbation terms for the triangular-well potential. a) First and b) second-order terms, A1 /NkT and A2 /NkT. Our simulation results are represented by symbols and the solid lines correspond to the TW-EoS.

V. Trejos et al. / Journal of Molecular Liquids 265 (2018) 337–346

341

Besides, the A2 shows interesting features. The range of values for A2 is very narrow, compared with the range of values that A1 takes, except for low densities when both A1 and A2 → 0. For very shortranged attractions, A2 is comparable to A1 , then it is expected that the contribution of A2 , and higher order terms, to the free energy would become more important for systems where k → 1. Another important feature of A2 is that it has only one minimum which is different of the SW case, which has two minimum as was reported by Espíndola et al. [18]. Finally, for the TW fluid the minimum of A2 occurs at a density close to the critical value, as we will see below. 4.2. Phase equlibria We classified the attraction ranges in two categories: long-ranges (1.75 ≤ k ≤ 2.5) and short-ranges (1.2 ≤ k < 1.75). Although in the literature all ranges longer than 2.0 are considered long-ranges, we use this classification just to make comparison with similar SW potentials where short-ranged attraction correspond to kSW ≤ 1.25, and given that the effective attraction range of a TW is almost one half of the one in an SW. In Fig. 3, we present the binodal line for particles interacting through a) long and b) short ranged attractions. We included in this plot other simulation results, previously reported, those are: squares [31], diamonds [32] and triangles [37]. For all cases, our simulation results are consistent with the other simulation results. Here, the performance of the TW-EoS is compared with the numerical results in order to show the reliability of the equation. For longranged attractions, we obtain almost the same results using both approaches but for short-ranged attractions they differ for about 10%. For the long-ranged attractions the TW-EoS reproduces the shape of the binodal obtained with simulations. Close to the critical point the TW-EoS overestimates the binodal as is expected for theories bases on a perturbation approaches, the critical temperature is overestimated by less than 10% however the localization of the critical density is well estimated. At high densities, we also observe small deviations, such deviations become small as the attraction range decreases. For short ranges, we see that the TW-EoS nicely reproduce the binodal even at high densities. However, for values of the attraction range k < 1.5, the localization of the critical density shifts to higher densities and the binodal get biased. The critical temperature in TW fluids follows a same trend as the one of SW fluids [19], Tc ∗ is almost a linear function of the attraction range for the values we considered in this work, as can be seen in Fig. 4(a). In that figure, we plot the simulation data, ours and from

Fig. 4. a) Critical temperature and b) critical density as a function of the triangularwell attraction range. Solid lines correspond to the theoretical predictions using the TW-EoS and symbols are MC simulation data obtained in this work and taken from other references, as is indicated. In panel a) we also included the temperature obtained from the relation: B∗2 (T) = −1.5, which is a good predictor for the critical temperature, see Ref. [25].

Refs. [31] (inverted triangles), [32] (squares), [37] (triangles) and [36] (diamonds) along with the TW-EoS prediction and the critical temperatures estimated by solving the equation B∗2 (Tc∗ ) = −1.5. At small attraction ranges, k < 1.5, second virial coefficient evaluated at critical value, takes a value of −1.5, this is expected from the extended law of corresponding states [25]. Besides, our prediction for Tc∗ is

Fig. 3. Vapor-liquid coexistence for triangular-well fluids. Cases presented: a) long ranges k = 1.75, 2.0 2.2 and 2.5, and b) short ranges k = 1.3, 1.35 1.4 and 1.5. Solid lines correspond to the theoretical predictions with the TW-EoS. The empty symbols are our Gibbs ensemble simulation results and solid symbols correspond to MC simulations from Ref. [32] (squares), Ref. [31] (diamonds) and Ref. [37] (triangles).

342

V. Trejos et al. / Journal of Molecular Liquids 265 (2018) 337–346

very close to the values obtained in MC simulation in the full range from 1.2 to 2.6. Our theoretical prediction give us an overestimation almost constant in all the range of k values. In Fig. 4(b), we plot the q∗c obtained from simulation and TW-EoS. At long attraction ranges the deviation is almost constant the TW-EoS underestimated the critical value for about 10%. For small attraction ranges the trend inverts and now the TW-EoS overestimated the critical value. This combined effect could be due to the inaccuracy of the A2 term, which is biased to higher densities for long attraction ranges and biased to small densities for short attraction ranges. In Fig. 5, we compare the theoretical and simulation results for vapor pressures of TW fluids with ranges k = 1.414, 1.5, 1.736, 1.95, and 2.5. Here, for longest attraction range we observe a good agreement with the simulation data, however as k becomes small our prediction only gives us a qualitative agreement. This inconsistence could be explained in terms of the quality of the vapor pressures simulation data reported in the Ref. [33]. The results presented in Fig. 5 are frequently known in the literature as the ClausiusClapeyron representation of the vapor pressure and this results are useful to estimate the enthalpy of vaporization. 4.3. Excess energy and pressure The derived EoS for TW potentials agrees reasonably well with the MC simulation data for all the interaction ranges studied. In this section, we analyze the performance of the EoS at supercritical temperatures where we expect even a better agreement. Fig. 6(a) presents the excess energy vs. q∗ for particles at T∗ = 2.0 with different attraction ranges, obtained from MC simulations reported in Ref. [32] (symbols) and our EoS (lines). For all attraction ranges, the EoS presents a small deviation from the simulation results however it reproduces well the behavior in each case. Then at this temperature, the contribution of the higher order perturbation terms becomes very small compared with the one of the first term. In Fig. 6(b), we plot the energy for a long-ranged potential, k = 2.5, at different temperatures. Again, the MC data was taken from Ref. [32] (symbols). We observe a nice agreement for all temperatures. Surprisingly at T∗ = 2.0 the agreement is good but the systems must be already phase separated, according to Ref. [32] Tc∗ = 2.34 and the coexisting densities are around 0.1 and 0.5 for the gas and liquid, respectively. In Fig. 7, we show the P∗ corresponding to some representative values of k, at different temperatures. Symbols are the MC results from Ref. [31]. For k = 1.2 and 1.5, the EoS predicts very well

Fig. 6. Excess energy for triangular-well fluids. a) Fluids at T∗ = 2.0 with different ranges k = 1.2, 1.5, 1.75, and 2.0 and b) fluids with the same attraction range k∗ = 2.5 but different temperatures T∗ = 2.0, 2.5, 3.0, and 4.0. Solid lines correspond to the theoretical predictions using the TW-EoS. The MC simulation results, were taken from Ref. [32] (k = 2.5) and Ref. [31] (k = 1.2, 1.5, 1.75 and 2.0), are symbols.

the MC results for T∗ > 1, which are temperatures well above the corresponding critical values: Tc∗ = 0.475 and 0.75 respectively. For the isotherm T∗ = 0.5, the EoS follows a similar trend as the simulation however the agreement is no that good because the simulation data were obtained at densities at which we did not consider for the EoS fitting. For k = 1.75 and 2.0, the critical temperatures are around 1.0 and 1.125, respectively (see Fig. 7), for those states outside of the binodal the EoS performs very well even for those temperatures close to the Tc∗ , this is a very nice feature for our EoS. For the largest temperatures, the EoS is indistinguishable from the simulation. Also, we point out that our prediction deviate from MC data when the densities are close to the fluid-solid coexistence. This occurs typically in the range 0.9 < q∗ < 1.1, as was reported in the Ref. [34]. Besides, the calculation for A2 was done by using a fitted expression valid in the range of densities q∗ ≤ 0.9. Then we did not expected our TW-EoS to be valid for metastable fluids at high densities. 4.4. Modeling of real substances

Fig. 5. Vapor pressure for triangular-well fluids of range: k = 1.414, 1.5, 1.736, 1.95, and 2.5. The solid lines are the theoretical predictions with the TW-EoS and the symbols are the simulation data of Adhikari and Kofke [33].

In this section, we use our EoS to find the parameters that model some simple fluids as well as predict their corresponding phase equilibrium and vapor pressure. For this comparison, we choose the following substances: methane (CH4 ), oxygen (O2 ), fluoromethane (R41) and hydrogen sulfide (H2 S). For all fluids, the parameters have been obtained from a best fit to experimental vapor pressures and saturated liquid densities by using

V. Trejos et al. / Journal of Molecular Liquids 265 (2018) 337–346

343

Fig. 7. Pressure for triangular-well fluids with different attraction ranges at some given temperatures. Cases presented: (a) k = 1.2 and T∗ = 0.5, 1.0, 1.5, 2.0, (b) k = 1.5 and T∗ = 0.5, 1.0, 1.5, 2.0, (c) k = 1.75 and T∗ = 0.5, 1.0, 1.5, 2.0, (d) k = 2.0 and T∗ = 0.5, 1.0, 1.5, 2.0, Symbols are MC simulation data from Ref. [31] and solid lines correspond to the theoretical predictions with the TW-EoS.

the optimization method presented in Ref. [53]. The experimental data that is used as the fitting target was reported in Ref. [45] and those values were obtained trough a minimization process of the function,

Fob =

Nq 



exp

− qi qcal i exp

i=1

qi

2 +

NP 



exp 2

Pjcal − Pj exp

j=1

Pj

(27)

where Fob is the objective function, NP and Nq are the number of experimental reference data points for saturation pressure and satuexp exp rated liquid-density, respectively, Pi and qj are the experimental reference value of the vapor pressure and saturated liquid of each point, respectively, and Pical and qcal are the corresponding vapor j pressure and saturated liquid densities calculated with the EoS for a give selection of the molecular parameters. In the minimization of the objective function Eq. (27), we use the Levenberg Marquardt method. The optimized values for the HS diameter (s), the range of the TW attractive interaction (k) and the energy (4) are listed in Table 1. The set of parameter values for CH4 agrees with that reported by Betancourt et al. in Ref. [32]. In Figs. 8 and 9, we show the binodal line and vapor pressures, respectively, for the substances listed in Table 1 obtained though the TW-EoS. For R41 and H2 S, the binodal match very well the experimental data, except in the vicinity of the critical point. Such feature is expected since the same happens with the simulation data. For CH4 and O2 , our fitting describes reasonable well the experimental. In those cases, the experimental binodal has a shape slightly different from the previous ones, that could indicate that the interaction could

no be well represented by a solely radial interaction. However, for all the analyzed substances, we observe a good agreement between vapor pressure obtained via the EoS and the experimental data, even though such good agreement is not clear in the binodal line for CH4 and O2 . Finally, a VLE comparison between molecular simulation data, experimental data and TW-EoS predictions for argon (Ar) and xenon (Xe) is shown in Fig. 10. Solid lines are the TW-EoS, filled symbols correspond to the experimental data taken Chemistry WebBook [45] and empty symbols are simulation from Ref. [49]. Here, the binodal line predictions for Ar and Xe were obtained by using the molecular parameters (s, 4 and k) reported by Bárcenas et al. [49]. As can be seen in this figure, the TW-EoS match very well with the experimental data, in a broad range of values of densities and temperatures. In both cases (k = 2.045 and 2.030), noticeable deviations in the liquid densities are observed at low temperatures when the results are

Table 1 Optimized TW-EoS parameters for fluoromethane (R41), hydrogen sulfide (H2 S), methane (CH4 ), and oxygen (O2 ) fluids obtained by fitting to experimental data liquid density and vapor pressure. Molecular parameters are: the HS diameter (s), the attraction range of the TW interaction (k), and the energy-depth in the TW potential (4). Also we included the critical values for the reduced temperature and pressure. Fluid

k

s [nm]

4/kB [K]

Tc∗

Pc∗

CH4 O2 R41 H2 S

1.8837 1.8512 1.6761 1.7574

3.7371 3.3890 3.8824 3.7441

168.2307 141.7071 376.8620 385.5400

1.1691 1.1263 0.9126 1.0084

0.1188 0.1162 0.1041 0.1092

344

V. Trejos et al. / Journal of Molecular Liquids 265 (2018) 337–346

Fig. 8. Vapor-liquid coexistence for the following substances are: (a) fluoromethane (R41), hydrogen sulfide (H2 S), and (b) methane (CH4 ), oxygen (O2 ). The lines correspond to the TW-EoS for fluids characterized by the parameters of Table 1, and symbols are the experimental data taken from the NIST Chemistry WebBook [45].

Fig. 9. Vapor pressure curves for the following substances are: (a) fluoromethane (R41), hydrogen sulfide (H2 S), and (b) methane (CH4 ), oxygen (O2 ). The lines correspond to the TW-EoS for fluids characterized by the parameters of Table 1, and symbols are the experimental data taken from the NIST Chemistry WebBook [45].

compared with experimental data. As was discussed in Section 4.2, the critical temperature is overestimated as expected in this kind of PT (Fig. 10). 5. Conclusions and perspectives In this paper, we present an accurate EoS to predict the thermodynamic properties of TW fluids based on the BH PT [2,28]. This TW-EoS performs very well against Monte Carlo simulation for systems of particles whose attraction range lies between 1.2 and 2.6, then this equation could be used to characterize the properties of model systems in colloidal and simple fluids regimes. In particular, we used the EoS to model the molecular interaction of a few simple fluids, see Table 1. In all cases, our calculation for the saturation pressure shows a good agreement when it is compared with the experimental data [45]. However, our prediction for the binodal curve agrees with the experimental data for fluoromethane and hydrogen sulfide but in the of methane and oxygen there are some discrepancies liquid branch predictions. Here, we limited the use of the TW-EoS to study particles with a radial interactions but it can be also used to model chains, associating

Fig. 10. Vapor-liquid coexistence for Argon and Xenon modeled as TW fluids. The molecular parameters s, 4 and k are reported by Bárcenas et al. [49]. Here, we plot our results using the TW-EoS along with the simulation of Ref. [49] and the experimental data, taken from the NIST Chemistry WebBook [45].

V. Trejos et al. / Journal of Molecular Liquids 265 (2018) 337–346

and mixture of fluids, as it has been done for other pair potentials, as the LJ [54-56] and the SW [43,44]. On the other hand, the TW-EoS can be also used for the prediction of adsorption isotherms of monomeric fluids with associating sites [57,58]. Additionally, the procedure presented here can be used to obtain an EoS for very short-ranged potentials to be applied to the study of colloidal systems [19]. Acknowledgments N. E. V.-P. acknowledges the financial support provided by CONACYT through a postdoctoral scholarship (Convocatoria de Estancias Posdoctorales en el Extranjero 2017). We also thank S. F.Gerstenmaier (UG) for high performance computer facilities (Fondos Mixtos CONACYT-CONCYTEG 2011, PROMEP 2010) that have contributed to the research results reported within this paper. Appendix A. Analytical evaluation of the integral I(r)

Fig. A.11. Numerical and calculated function given by Eq. (19) as a function of the reduced density q∗ in the range of values 1.20 ≤ k ≤ 2.6.

There are several expressions for f(k) that satisfied Eq. (19). In this work, we have used the most simple approximation in order to obtain a useful analytic solution for the second term A2 . In this way, the following steps were performed: i) we solve numerically the integral equations A1 by using the HS radial distribution function (RDF) proposed by Chang et al. [29], ii) the obtained results for A1 were compared with MC simulation results, iii) we solve numerically the integral term Ir by using the RDF by Chang et al. [29] for different k values and finally, iv) the analytical and numerical results for Ir and A1 at different values of k were compared in order to obtain an expression for f(k). As can be seen in Fig. A.11 the simple function f(k) = 7k/11 satisfied Eq. (19) in the range of values 1.2 ≤ k ≤ 2.6. Here, we can observe that some deviations are observed at low values of k and high values of q∗ . Therefore, more complex functions for f(k) would be proposed for widely ranges of k and q∗ in order to improve the concordance of Eq. (19). We used the more simple expression in order to show the most simple solution of the term A2 .

References [1] P. Bolhuis, M. Hagen, D. Frenkel, Isostructural solid-solid transition in crystalline systems with short-ranged interaction, Phys. Rev. E 50 (1994) 4880–4890.

345

[2] J.A. Barker, D. Henderson, Perturbation theory and equation of state for fluids: the square well potential, J. Chem. Phys. 47 (8) (1967) 2856–2861. [3] W.R. Smith, D. Henderson, J.A. Barker, Approximate evaluation of the second-order term in the perturbation theory of fluids, J. Chem. Phys. 53 (1970) 508–515. [4] W.R. Smith, D. Henderson, J.A. Barker, Perturbation theory and the radial distribution function of the squarewell fluid, J. Chem. Phys. 55 (8) (1971) 4027–4033. [5] D. Henderson, J.A. Barker, W.R. Smith, Calculation of the contact value of the first-and second-order terms in the perturbation expansion of the radial distribution function for the square-well potential, J. Chem. Phys. 64 (10) (1976) 4244–4245. [6] D. Henderson, O.H. Scalise, W.R. Smith, Monte Carlo calculations of the equation of state of the square well fluid as a function of well width, J. Chem. Phys. 72 (4) (1980) 2431–2438. [7] D.D. Carley, Equations of state for a square well gas from a parametric integral equation, J. Chem. Phys. 67 (3) (1977) 1267–1272. [8] D.D. Carley, A.C. Dotson, Integral equation and perturbation method of calculating thermodynamic functions for a square-well fluid, Phys. Rev. A 23 (3) (1981) 1411–1418. [9] D.D. Carley, Thermodynamic properties of a squarewell fluid in the liquid and vapor regions, J. Chem. Phys. 78 (9) (1983) 5776–5781. [10] A. de Lonngi, F. de1 Rio, Square well perturbation theory of fluids, Mol. Phys. 48 (2) (1983) 293–313. [11] A. de Lonngi, F. de1 Rio, Square-well perturbation theory for the structure of simple fluids, Mol. Phys. 56 (3) (1985) 691–700. [12] F. de1 Rio, L. Lira, Properties of the square-well fluid of variable width I. Short-range expansion, Mol. Phys. 61 (2) (1987) 275–292. [13] F. de1 Rio, L. Lira, Properties of the square-well fluid of variable width. II. The mean field term, J. Chem. Phys. 87 (12) (1987) 7179–7183. [14] A.L. Benavides, F. de1 Rio, Properties of the square-well fluid of variable width III. Long-range expansion, Mol. Phys. 68 (5) (1989) 983–1000. [15] A.L. Benavides, J. Alejandre, F. de1 Rio, Properties of the square-well fluid of variable width: IV. Molecular dynamics test of the van der Waals and long-range approximations, Mol. Phys. 74 (2) (1991) 321–331. [16] A. Gil-Villegas, F. de1 Rio, A.L. Benavides, Deviations from correspondingstates behavior in the vapor-liquid equilibrium of the square-well fluid, Fluid Phase Equilib. 119 (2) (1996) 97–112. [17] E. Schöll-Paschinger, A.L. Benavides, R. Castañeda-Priego, Vapor-liquid equilibrium and critical behavior of the square-well fluid of variable range: a theoretical study, J. Chem. Phys. 123 (23) (2005) 234513. [18] R. Espíndola-Heredia, F. del Río, A. Malijevsky, Optimized equation of the state of the square-well fluid of variable range based on a fourth-order free-energy expansion, J. Chem. Phys. 130 (2) (2009) 024509. [19] N.E. Valadez-Pérez, A.L. Benavides, E. Schöll-Paschinger, R. Castañeda-Priego, Phase behavior of colloids and proteins in aqueous suspensions: theory and computer simulations, J. Chem. Phys. 137 (8) (2012) 084905. [20] Rotenberg, Monte Carlo equation of state for hard spheres in an attractive square well, J. Chem. Phys. 43 (4) (1965) 1198–1201. [21] F. Lado, W.W. Wood, N dependence in Monte Carlo studies of the square-well system, J. Chem. Phys. 49 (9) (1968) 4244–4245. [22] Y. Rosenfeld, R. Thieberger, Monte Carlo and perturbation calculations for the square well fluid: dependence on the square well range, J. Chem. Phys. 63 (5) (1975) 1875–1877. [23] K.D. Scarfe, I.L. McLaughlin, A.F. Collings, The transport coefficients for a fluid of square-well rough spheres: comparison with methane, Mol. Phys. 65 (8) (1976) 2991–2994. [24] L. Vega, E. de Miguel, L.F. Rull, G. Jackson, I.A. McLure, Phase equilibria and critical behavior of squarewell fluids of variable width by Gibbs ensemble Monte Carlo simulation, J. Chem. Phys. 96 (3) (1992) 2296–2305. [25] M.G. Noro, D. Frenkel, Extended corresponding-states behavior for particles with variable range attractions, J. Chem. Phys. 113 (8) (2000) 2941–2944. [26] F. Platten, N.E. Valadez-Pérez, R. Castañeda-Priego, S.U. Egelhaaf, Extended law of corresponding states for protein solutions, J. Chem. Phys. 142 (17) (2015) 174905. [27] J. Largo, J. Solana, A simplified perturbation theory for equilibrium properties of triangular-well fluids, Physica A 284 (1) (2000) 68–78. [28] J.A. Barker, D. Henderson, Perturbation theory and equation of state for fluids: II. A successful theory of liquids, J. Chem. Phys. 47 (1967) 4714–4721. [29] J. Chang, S.I. Sandler, A real function representation for the structure of the hard-sphere fluid, Mol. Phys. 81 (1994) 735–744. [30] J. Largo, J. Solana, Theory and computer simulation of the first- and second-order perturbative contributions to the free energy of square-well fluids, Mol. Simul. 29 (6–7) (2003) 363–371. [31] F.F. Betancourt-Cárdenas, L.A. Galicia-Luna, S.I. Sandler, Thermodynamic properties for the triangular-well fluid, Mol. Phys. 105 (23–24) (2007) 2987–2998. [32] F. Betancourt-Cárdenas, L. Galicia-Luna, A. Benavides, J. Ramírez, E. Schöll-Paschinger, Thermodynamics of a long-range triangle-well fluid, Mol. Phys. 106 (1) (2008) 113–126. [33] J. Adhikari, D.A. Kofke, Monte Carlo and cell model calculations for the solid-fluid phase behaviour of the triangle-well model, Mol. Phys. 100 (10) (2002) 1543–1550. [34] S. Zhou, Thermodynamics and phase behavior of a triangle-well model and density-dependent variety, J. Chem. Phys. 130 (1) (2009) 014502.

346

V. Trejos et al. / Journal of Molecular Liquids 265 (2018) 337–346

[35] L. Rivera, M. Robles, M.L. de Haro, Equation of state and liquid-vapour equilibrium in a triangle-well fluid, Mol. Phys. 110 (11–12) (2012) 1317–1323. [36] M. Bárcenas, G. Odriozola, P. Orea, Coexistence and interfacial properties of triangle-well fluids, Mol. Phys. 112 (16) (2014) 2114–2121. [37] M. Bárcenas, V. Castellanos, Y. Reyes, G. Odriozola, P. Orea, Phase behaviour of short range triangle well fluids: a comparison with lysozyme suspensions, J. Mol. Liq. 225 (2017) 723–729. [38] A. Sengupta, P. Behera, J. Adhikari, Molecular simulation study of triangle-well fluids confined in slit pores, Mol. Phys. 112 (15) (2014) 1969–1978. [39] A. Sengupta, J. Adhikari, Fluid phase equilibria of triangle-well fluids confined inside slit pores: a transition matrix Monte Carlo simulation study, J. Mol. Liq. 221 (Supplement C) (2016) 1184–1196. [40] H. Guérin, Improved analytical thermodynamic properties of the triangular-well fluid from perturbation theory, J. Mol. Liq. 170 (2012) 37–40. [41] H. Guérin, Unified SAFT-VR theory for simple and chain fluids formed of square-well, triangular-well, Sutherland and Mie segments, J. Mol. Liq. 203 (2015) 187–197. [42] M.C. Somasekhara Reddy, A.K. Murthy, Perturbation theory for a microemulsion with triangular well potential, Pramana 20 (3) (1983) 217. [43] A. Gil-Villegas, A. Galindo, P.J. Whitehead, S.J. Mills, G. Jackson, A.N. Burgess, Statistical associating fluid theory for chain molecules with attractive potentials of variable range, J. Chem. Phys. 106 (1997) 4168–4186. [44] B.H. Patel, H. Docherty, S. Varga, A. Galindo, G.C. Maitland, Generalized equation of state for square-well potentials of variable range, Mol. Phys. 103 (1) (2005) 129–139. [45] E. Lemmon, M. McLinden, D. Friend, Thermophysical Properties of Fluid Systems, National Institute of Standards and Technology, Gaithersburg MD, 20899, 2017. [46] L.L. Lee, Molecular Thermodynamics of Nonideal Fluids, Butterworth Publishers, Boston, MA, 1988.

[47] N.F. Carnahan, K.E. Starling, Equation of state for nonattracting rigid spheres, J. Chem. Phys. 51 (2) (1969) 635–636. [48] J.A. Barker, D. Henderson, What is liquid? Understanding the states of matter, Rev. Mod. Phys. 48 (1976) 587–671. [49] M. Bárcenas, Y. Reyes, A. Romero-Martínez, G. Odriozola, P. Orea, Coexistence and interfacial properties of a triangle-well mimicking the Lennard-Jones fluid and a comparison with noble gases, J. Chem. Phys. 142 (2015)pp. 074706 1–5. [50] A.Z. Panagiotopoulos, Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble, Mol. Phys. 61 (4) (1987) 813–826. [51] D. Frenkel, B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Computational Science Series Elsevier Science. 2001. [52] F. Wegner, Corrections to scaling laws, Phys. Rev. B 5 (1972) 4529. [53] J.A. López, V.M. Trejos, C.A. Cardona, Objective functions analysis in the minimization of binary VLE data for asymmetric mixtures at high pressures, Fluid Phase Equilib. 248 (2006) 147–157. [54] W.G. Chapman, K.E. Gubbins, G. Jackson, M. Radosz, SAFT: Equation-of-state solution model for associating fluids, Fluid Phase Equilib. 52 (1989) 31–38. [55] W.G. Chapman, K.E. Gubbins, G. Jackson, M. Radosz, New reference equation of state for associating liquids, Ind. Eng. Chem. Res. 29 (1990) 1709–1721. [56] C. Avendaño, T. Lafitte, A. Galindo, C.S. Adjiman, G. Jackson, E. Müller, SAFT-c force field for the simulation of molecular fluids. 1. A single-site coarse grained model of carbon dioxide, J. Phys. Chem. B 115 (2011) 11154–11169. [57] A. Martínez, V.M. Trejos, A. Gil-Villegas, Predicting adsorption isotherms for methanol and water onto different surfaces using the SAFT-VR-2D approach and molecular simulation, Fluid Phase Equilib. 449 (2017) 207–216. [58] V.M. Trejos, J. Quintana-H, Thermodynamic properties of confined square-well fluids with multiple associating sites, J. Chem. Phys. 148 (7) (2018) 074703.