Scaled Particle Theory for Multicomponent Hard Sphere Fluids

Article pubs.acs.org/JPCB

Scaled Particle Theory for Multicomponent Hard Sphere Fluids Confined in Random Porous Media W. Chen,†,‡ S. L. Zhao,§ M. Holovko,∥ X. S. Chen,⊥ and W. Dong*,†,⊥ †

Université de Lyon, CNRS, Ecole Normale Supérieure de Lyon, Université Lyon 1, Laboratoire de Chimie, UMR 5182, 46 Allée d’Italie, 69364 Lyon Cedex 07, France ‡ Computer Network Information Center, Chinese Academy of Sciences, P.O. Box 349, 100190 Beijing, China § State Key Laboratory of Chemical Engineering, East China University of Science and Technology, 130 Meilong Road, 200237 Shanghai, China ∥ Institute for Condensed Matter Physics, National Academy of Sciences, 1 Svientsitskii Street, 79011 Lviv, Ukraine ⊥ State Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, P.O. Box 2735, 100190 Beijing, China S Supporting Information *

ABSTRACT: The formulation of scaled particle theory (SPT) is presented for a quite general model of fluids confined in a random porous media, i.e., a multicomponent hard sphere (HS) fluid in a multicomponent hard sphere or a multicomponent overlapping hard sphere (OHS) matrix. The analytical expressions for pressure, Helmholtz free energy, and chemical potential are derived. The thermodynamic consistency of the proposed theory is established. Moreover, we show that there is an isomorphism between the SPT for a multicomponent system and that for a one-component system. Results from grand canonical ensemble Monte Carlo simulations are also presented for a binary HS mixture in a one-component HS or a one-component OHS matrix. The accuracy of various variants derived from the basic SPT formulation is appraised against the simulation results. Scaled particle theory, initially formulated for a bulk HS fluid, has not only provided an analytical tool for calculating thermodynamic properties of HS fluid but also helped to gain very useful insight for elaborating other theoretical approaches such as the fundamental measure theory (FMT). We expect that the general SPT for multicomponent systems developed in this work can contribute to the study of confined fluids in a similar way.

1. INTRODUCTION

rule for the porosity has a very appealing form with the ratio of packing fraction as weight rather than the species concentration. To appraise the accuracy of the theory, grand canonical ensemble Monte Carlo (GCEMC) simulations are carried out for a binary HS mixture confined in either a hard sphere matrix or an overlapping hard sphere (OHS) matrix under a variety of conditions. The comparison between theory and simulation shows that the SPT results for the chemical potential of the larger species are systematically less accurate than those of the smaller species. The article is arranged as follows. The theory and the information about the GCEMC simulations are presented in section 2. In section 3, the simulation results and the comparison between theory and simulation are presented and discussed in details. We conclude with a few remarks in the last section.

In 2009, Holovko and Dong initiated the extension of the very successful scaled particle theory (SPT)1−3 to a hard sphere (HS) fluid confined in random porous media.4 Subsequent refinement and reformulation have been made by Chen et al.5 and Patsahan et al.6 Despite of the much effort devoted to the study of fluids confined in random porous media until now,7−14 much less attention has been paid to confined mixtures. In this article, we present the formulation of SPT for a multicomponent HS fluid confined in a multicomponent HS matrix. A remarkable isomorphism is found between some thermodynamic quantities, e.g., pressure and Helmholtz free energy, for the general multicomponent system (multicomponent fluid and multicomponent matrix) and those for a one-component fluid in a one-component matrix. This isomorphism shows that the SPT naturally has the structure of a one-fluid theory. The mixing rules for some properties resulting from SPT take the usual form with the species concentrations as the weights, whereas the mixing © 2016 American Chemical Society

Received: March 22, 2016 Revised: May 10, 2016 Published: May 11, 2016 5491

DOI: 10.1021/acs.jpcb.6b02957 J. Phys. Chem. B 2016, 120, 5491−5504

Article

The Journal of Physical Chemistry B n

2. METHOD 2.1. Theory. We consider here a quite general model of an ncomponent hard sphere (HS) fluid in an m-component porous matrix. Negative integers and zero will be used as indexes to denote matrix species, i.e., α = −m + 1, ..., −1, 0, whereas positive integers are used as indexes to denote fluid species, i.e., α = 1, 2, ..., n. There are Nα particles for the species α. With such a notation, the Hamiltonian of the considered fluid-matrix system can be written as n

0

Nκ

ΔH =

The chemical potential of the scaled particle is given by μsα = −kT ln

Nλ

1 2

n

n

Nκ

Nλ

∑∑∑

∑

κ=1 λ=1 i=1

j=1 j ≠ i (if κ = λ)

uκλ(|r iκ − r λj |) (1)

where the indices and superscripts in Greek denote species, the indices in English denote different particles of a given species, and the various interaction potentials, i.e., uκλ, will be given more explicitly later. When a scaled particle of a given fluid species, e.g., α whose position is denoted by rs, is inserted into the above system, the Hamiltonian becomes

H Nα + 1 = H + ΔH

Z N1,..., Nα + 1,...Nn(r N −m+1, ..., r N0) Z N1,..., Nα ,...Nn(r N −m+1, ..., r N0)

{m}

(4)

where k is the Boltzmann constant, T is the temperature, ZN1,...,Nα+1,...Nn(rN−m+1, ..., rN0) and ZN1,...,Nα,...Nn(rN−m+1, ..., rN0) are the partition functions of the system with or without, respectively, the inserted scaled particle for a particular realization of matrix configurations, and ⟨...⟩{m} denotes the average over matrix configurations. rNα (α = −m + 1, ..., −1, 0) is a shorthand notation of all the position vectors of the matrix particles of species α. A few straightforward algebraic arrangements lead to

κ =−m + 1 λ = 1 i = 1 j = 1

+

(3)

κ =−m + 1 i = 1

∑ ∑ ∑ ∑ uκλ(|r iκ − r λj |)

H=

Nκ

∑ ∑ uακ(|rs − r iκ|)

Z N1,..., Nα + 1,...Nn(r N −m+1, ..., r N0) Z N1,..., Nα ,...Nn(r N −m+1, ..., r N0) =

(2)

1 (Nα + 1)Λα 3

∫ drs pα (rs, R s; r N

, ..., r N0)

−m + 1

(5)

where Λα is the thermal wavelength of species α and

where pα (rs, R s ; r N −m+1, ..., r N0) =

1 n N −m + 1 Z N1,..., Nα ,...Nn(r , ..., r N0) ∏γ = 1 (Nγ !Λ γ 3Nγ)

⎤

γ=1

(6)

(10)

where Rκ and Rλ are respectively the radius of the hard sphere of species κ and that of species λ. The interactions potentials for a system with an n-component HS fluid confined in an mcomponent OHS matrix are characterized by eqs 11 and 12.

(8)

In deriving eq 7 from eq 4, we have changed the order of the average over matrix configurations and the logarithm. As for a one-component fluid,4 this change of order is rigorously valid for a small scaled particle. To derive a consistent SPT for a confined fluid, it has been shown that one needs to factorize the insertion probability into two parts.6 For the multicomponent systems considered here, we make the similar factorization pα (R s) = pIα (R s) pIIα (R s)

⎡

κ = −m + 1, ..., n; λ = −m + 1, ..., n

(7)

where ρα = Nα/V is the density of species α and pα (R s) = ⟨pα (rs, R s ; r N −m+1, ..., r N0)⟩{m}

γ

⎧∞ |r κ − r λ| < R + R i j κ λ ⎪ uκλ(|r iκ − r λj |) = ⎨ ⎪ 0 |r iκ − r λj | ≥ R κ + Rλ ⎩

is the insertion probability before the average over matrix configurations. Substituting eq 5 into eq 4, we obtain the following expression for the excess chemical potential μsex = μsα − kT ln(Λα 3ρα ) = −kT ln pα (R s) α

n

H ∫ ∏ dr N e−H /kT exp⎢⎣ −Δ ⎥ kT ⎦

uκλ(|r iκ − r λj |) = 0 κ = −m + 1, ..., 0; λ = −m + 1, ..., 0

uκλ(|r iκ

(9)

−

r λj |)

(11)

⎧∞ |r κ − r λ| < R + R i j κ λ ⎪ =⎨ ⎪ 0 |r iκ − r λj | ≥ R κ + Rλ ⎩

κ = −m + 1, ..., n; λ = −m + 1, ..., n

where pαI (Rs) is the probability of inserting a scaled particle into an empty matrix and pαII(Rs) is the conditional probability to insert a scaled particle into a matrix imbibed of fluid under the condition that the insertion is made in a matrix pore that is at least as large as the scaled particle. It is obvious that pαI (Rs) depends on the model for describing the matrix. In the present work, we consider two matrix models, i.e., an m-component HS matrix and an m-component overlapping HS matrix (OHS). The interactions potentials for a system with an n-component HS fluid confined in an m-component HS matrix are characterized by

(if κ ≤ 0, λ > 0; if κ > 0, λ ≤ 0; if κ > 0, λ > 0)

(12)

For the OHS matrix, pαI (Rs) can be calculated exactly, and we found the following analytic expression pIα (R s)

0 ⎡ ⎤ ⎢ = exp − ∑ ηκ (1 + R s/R κ )3 ⎥ ⎢⎣ κ =−m + 1 ⎥⎦

(13)

where 5492


Article

The Journal of Physical Chemistry B ηκ =

4πR κ 3ρκ

results, we use an ansatz similar to that we proposed previously for a one-component confined fluid, i.e., matching βμαII(Rs) = −ln pαII(Rs) with the following thermodynamics expression

(14)

3

For the HS matrix, it is straightforward to show that pαI (Rs) is the probability for inserting a HS of radius equal to Rs into an equilibrium m-component HS fluid if the HS matrix is made by quenching such an equilibrium fluid; thus, we have ⎡ μex (R s) ⎤ ⎥ pIα (R s) = exp⎢ − α I kT ⎦ ⎣

βW α(R s) = w0α + w1αR s +

the two functions and their first two derivatives at Rs = Rα). The matching yields

(15)

⎡ ⎤ ξ ⎥ w0α = −ln⎢1 − α F3 ⎢⎣ pI (R s = 0) ⎥⎦

where is the excess chemical potential of a scaled particle with a radius of Rs in an equilibrium m-component HS fluid, but it cannot be calculated exactly. SPT allows for obtaining the following approximate but analytic expression

⎛ξ R ⎞ bα I = 18⎜ M2 s ⎟ ⎝ ξM3 ⎠ ξMl =

π 6

⎡ α ⎢ d lnpI (R s) ⎢ξF3 dR s ⎣

(17)

2

−

(l = 1, 2, 3) (19)

λ =−m + 1

ξFl = (20)

ρλ

λ =−m + 1

(21)

36ξM1ξM2 πρM ξM3

(22)

πρM ξM32

dR s

R s= 0

⎫ ⎪ − 24ξF1⎬ ⎪ ⎭

(28)

π 6

n

∑ ρλ (2Rλ)l

(l = 1, 2, 3) (29)

λ=1

⎛ ξF3/φ ⎞2 4πR α 3βP 0 ⎟⎟ + + bα⎜⎜ 3ϕα ⎝ 1 − ξF3/φ0 ⎠

Now, for pαII(Rs), we obtain the following exact result n

∑ ηλ(1 + R s/Rλ)3 ,

φ0 = pIα (R s = 0)

(31)

ϕα =

pIα (R s

(32)

aα =

α ⎛ 6ξF2R α 12ξF1R α 2 6ξ R ⎞ d ln pI (R s) + − ⎜1 + F2 α ⎟R α ξF3 ξF3 ξF3 ⎠ dR s ⎝

λ=1

−min{R −m + 1 , ..., R n} ≤ R s ≤ 0

(30)

where

(23)

1 pIα (R s)

d lnpIα (R s)

⎞2 ⎤ ⎟⎥ ⎟⎥ R s= 0 ⎠ ⎥ ⎦

⎛ ξF3/φ0 ξ ⎞ = −ln pIα (R s = R α) − ln⎜⎜1 − F3 ⎟⎟ + aα φ0 ⎠ 1 − ξF3/φ0 ⎝

27ξM2 3

pIIα (R s) = 1 −

R s= 0

⎛ d lnpIα (R s) − ⎜⎜ dR s ⎝

βμαex = −ln pIα (R s = R α) + βW α(R s = R α)

0

BM =

⎤2 ⎥ − 6ξF2 ⎥ ⎦

The excess chemical potential of species α is given by

where

AM =

(27)

where

ξM3 ξM32 βP0 A 2BM 1 = + M + ρM 1 − ξM3 2 (1 − ξM3)2 3 (1 − ξM3)3

∑

⎤ ⎥ ⎥ R s= 0 ⎦

1 pIα (R s = 0)[1 − ξF3/pIα (R s = 0)]

+ 12ξF2

P0 in eq 16 is the pressure given by SPT for the equilibrium mcomponent HS fluid from which the matrix is obtained by quenching.

ρM =

R s= 0

⎧ ⎡ ⎪ ⎢ d 2 lnpIα (R s) ⎨ξF3⎢ dR s 2 ⎪ ⎢ ⎩ ⎣

(18)

ρλ (2Rλ)l

1 = 0) 1 − ξF3/pIα (R s = 0)

1 {pIα (R s = 0)[1 − ξF3/pIα (R s = 0)]}2

w2α =

0

∑

(26)

(16)

where β = 1/(kT) and 6ξM2R s + 12ξM1R s 2 ξM3

pIα (R s

⎡ α ⎢6ξ − ξ d ln pI (R s) F3 ⎢ F2 dR s ⎣

2

aα I =

1

w1α =

ξM3 1 − ξM3

⎛ ξ ⎞ 4πR s 3βP0 + bα I⎜ M3 ⎟ + 3 ⎝ 1 − ξM3 ⎠

(25)

where ϕα = pαI (matching here means requiring the equal values of

μex αI(Rs)

βμαexI (R s) = −ln(1 − ξM3) + aα I

3 1 α 2 4πR s βP w2 R s + 2 3ϕα

(24)

where ηλ = 4πR3λρλ/3 (packing fraction of species λ). For any reader who feels uncomfortable with the negative value of Rs as shown in eq 24, we recommend reading the clear description about the meaning of this given in the original paper of Reiss, Frisch, and Lebowitz (page 370 of ref 1). To obtain the SPT

= R α)

⎧⎡ α 1 ⎪⎢ d ln pI (R s) + ⎨⎢R α 2⎪ dR s ⎩⎣

R s= 0

⎫ ⎤ d 2 ln pIα (R s) ⎪ ⎥ 2 ⎬ 2 ⎥ − Rα dR s ⎪ R s= 0 ⎦ R s= 0 ⎭ 2

(33) 5493


Article

The Journal of Physical Chemistry B ⎡ d ln pIα (R s) 1 bα = ⎢⎢R α 2 dR s ⎣

R s= 0

⎤2 6ξF2R α ⎥ − ξF3 ⎥ ⎦

(1 + A)φ βP φ =− ln(1 − ξF3/φ) + ρ ξF3 (φ − φ0) ⎡φ ⎤ φ ×⎢ ln(1 − ξF3/φ) − 0 ln(1 − ξF3/φ0)⎥ ξF3 ⎣ ξF3 ⎦

(34)

It is noteworthy that pαI (Rs = 0) is in fact a quantity independent of the fluid species because it is simply the geometric porosity of the matrix, and this is why in the new symbol for representing this quantity in eq 30 and afterward, i.e., φ0, the superscript referring to a fluid species is removed. To obtain an equation containing only one unknown instead of eq 30, which contains two unknowns (chemical potential and pressure), we take the derivative with respect to total fluid density n ( ρ = ∑α = 1 ρα ) on the both sides of eq 30 by keeping the fluid composition unchanged (i.e., xα = ρα/ρ, α = 1, ..., n, constant) ∂(βμα ) ∂ρ

⎤⎫ φ φ ⎡φ ⎢ ln(1 − ξF3/φ) − 0 ln(1 − ξF3/φ0)⎥⎬ φ − φ0 ⎣ ξF3 ξF3 ⎦⎭

+

ξF3/φ0 2Bφ ⎧ 1 ⎨ − φ − φ0 ⎩ 2(1 − ξF3/φ0)2 1 − ξF3/φ0

⎛ ξF3/φ ⎞2 ⎛ ξF3/φ ⎞3⎤ 0 0 ⎟⎟ + 2bα⎜⎜ ⎟⎟ ⎥ × ⎜⎜ ⎝ 1 − ξF3/φ0 ⎠ ⎝ 1 − ξF3/φ0 ⎠ ⎥⎦

+

n

∂(βP) = ∂ρ

− (35)

∑ ρα α=1

∂(βμα ) ∂ρ

(36)

ξF3/φ0 ∂(βP) 1 1+A = + ∂ρ 1 − ξF3/φ 1 − ξF3/φ 1 − ξF3/φ0 A + 2B ⎡ ξF3/φ0 ⎤ ⎢ ⎥ 1 − ξF3/φ ⎢⎣ 1 − ξF3/φ0 ⎥⎦

+

∑ xαaα (38)

n

∑ xαbα 1 ξF3

(39) n

∑ α=1

4πR α 3ρα 3

ϕα −1

⎪

⎪

⎪

φ0 ξF3 φ0 ξF3 φ0 ξF3

ln(1 − ξF3/φ0) −

φ ⎡ 1 ⎢ φ − φ0 ⎢⎣ 1 − ξF3/φ0

⎤ ⎛ φ ⎞2 ⎡ φ ⎟⎢ ln(1 − ξF3/φ0)⎥+ ⎜⎜ ln(1 − ξF3/φ) ⎥⎦ ⎝ φ − φ0 ⎟⎠ ⎣ ξF3 ⎫ ⎤⎪ ln(1 − ξF3/φ0)⎥⎬ ⎦⎪ ⎭

(41)

(42)

A technical point for taking properly the limit φ → φ0 of the RHS of eq 41: To remove the singularity (φ − φ0)−1 (l = 1, 2, 3), one has to expand φln(1 − ξF3/φ)/ξF3 as a Taylor series around φ0 to the corresponding order. As pointed out above, it is impossible to obtain a differential equation containing just the chemical potential of a given species from eqs 35 and 36; we have to calculate the chemical potential in another way for a multicomponent system. For this, we first calculate the Helmholtz free energy from the pressure with the help of the following expression

(37)

n

α=1

⎪

2 2B (ξF3/φ0) + 3 (1 − ξF3/φ0)3

where

α=1

⎪

ξF3/φ0 βP SPT2a 1 A = + 1 − ξF3/φ0 2 (1 − ξF3/φ0)2 ρ

2

⎡ ξF3/φ ⎤3 2B 0 ⎢ ⎥ + 1 − ξF3/ϕ ⎢⎣ 1 − ξF3/φ0 ⎥⎦

⎪

We point out that the result given in the above equation for a system with an n-component fluid in an m-component matrix is isomorphic to that for a one-component fluid in a onecomponent matrix6 (see Appendix A for more details). In the case of a one-component fluid in a one-component matrix, it has been shown that the original version of SPT does not give the most accurate numerical results and some variants were derived.6 Here, we will examine also similar variants. When φ is replaced by φ0 in every term on the right-hand side (RHS) of eq 41 (more precisely, the limit φ → φ0 is taken), we obtain the following expression for pressure of the variant named SPT2a (we keep here the same designation of the variants as in ref 6).

and we see immediately that it does not allow for obtaining an equation for the chemical potential of one species from eq 35. Nevertheless, the combination of eqs 35 and 36 leads readily to

φ −1 =

−

−

In the case of a one-component fluid, we can eliminate one unknown, either μα or P, from eq 35 by using Gibbs−Duhem equation. For a multicomponent fluid, Gibbs−Duhem equation becomes

B=

φ (A + 2B)φ ⎧ 1 ⎨ + 0 ln(1 − ξF3/φ0) φ − φ0 ⎩ 1 − ξF3/φ0 ξF3

⎡ ξF3/φ0 1⎢ 1 + (1 + aα) = + (aα + 2bα) ⎢ 1 − ξF3/φ0 ρ⎣

4πR α 3 ∂(βP) + 3ϕα ∂ρ

A=

+

βF =ρ V

(40)

It is noteworthy that aα and bα as well as A and B are only functions of the concentration, i.e., xα, α = 1, ..., n, but independent of the total density, ρ. Integrating eq 37 with respect to ρ by keeping the concentration constant, we find finally

⎛

⎞

∫ dρ′ ρ1′ ⎜⎝ βρP′ ⎟⎠

(43)

We carry out this integration by choosing a particular integration path that consists of keeping the concentration (xα, α = 1, ..., n) constant and obtain 5494


Article

The Journal of Physical Chemistry B βF SPT2a βF id = − V V

n

∑ ργ ln

n

pIγ (R s

C=

= R γ)

γ=1

2⎤ B ⎛ ξF3/φ0 ⎞ ⎥ ⎟⎟ + ⎜⎜ 3 ⎝ 1 − ξF3/φ0 ⎠ ⎥⎦

n

E=

⎡ n ⎤ βF id = ρ⎢∑ xγ ln(Λ γ 3ργ ) − 1⎥ ⎢⎣ γ = 1 ⎥⎦ V

(45)

⎡ ∂(β FSPT2a /V ) ⎤ ⎥ βμαSPT2a = ⎢ = ln(Λα 3ρα ) ⎢⎣ ⎥⎦ ∂ρα {ρ , κ ≠ α}

n

κ

μ̃ =

⎛ ξF3/φ0 ξ ⎞ = R α) − ln⎜⎜1 − F3 ⎟⎟ + aα 1 − ξF3/φ0 φ0 ⎠ ⎝

2

⎡ ⎡ (2R )3 ξF3/φ0 (2R α)2 ⎤⎢ α ⎢ ⎥ − ⎢ ⎢⎣ (ξF3/φ0) (ξF2/φ0) ⎥⎦⎣ 2(1 − ξF3/φ0)

ξ ⎞ βP SPT2b βP SPT2a φ ⎛ = − ln⎜1 − F3 ⎟ ρ ρ ξF3 ⎝ φ⎠ + (46)

where

eα = R αD1 +

D1 =

(R D )2 ξF2 + α 1 2 ξF3

(47)

R α 2D2 2

(48)

D2 =

R s= 0

R s= 0

ξF3

⎛ ξ ⎞ ln⎜⎜1 − F3 ⎟⎟ φ0 ⎠ ⎝

(54)

⎡ ⎛ ξ ⎞ βF SPT2b βF SPT2a φ⎞ ⎛ = + ρ ⎢ − ⎜1 − ⎟ ln⎜1 − F3 ⎟ ⎢⎣ ⎝ ξF3 ⎠ ⎝ φ⎠ V V

(49)

⎛ φ ⎞ ⎛ ξ ⎞⎤ + ⎜1 − 0 ⎟ ln⎜⎜1 − F3 ⎟⎟⎥ ξF3 ⎠ ⎝ φ0 ⎠⎥⎦ ⎝

d 2 ln pIα (R s) dR s 2

φ0

In eq 54, we express PSPT2b in terms of PSPT2a because this allows us to have a more compact expressions for different thermodynamics functions, i.e., pressure, free energy, and chemical potential. In the following, we use also similar expressions for the other variants of SPT2. Integration pressure allows for obtaining the Helmholtz free energy.

d ln pIα (R s) dR s

(53)

is isomorphic to that of the chemical potential for a onecomponent fluid, and this holds for SPT2a and all the variants to be given in the following sections. One can check also that the pressure given in eq 42 and the chemical potential given in eq 46 satisfy Gibbs−Duhem equation (see Appendix B for more details). Keeping φ only in the first term on the RHS of eq 41 and taking the limit φ → φ0 in the other terms yields the following expression for the pressure of the variant SPT2b.

⎤ ⎡ 4πR 3ρC 3D1(ξF2/φ0)2 α × ⎢ − cα ⎥ + 2(ξF3/φ0) ⎦ ⎣ 3ξF3

cα = −6R α 2D1

∑ xαμα α=1

⎛ ξF3/φ ⎞2 ξF3/φ0 4πR α 3βP SPT2a 0 ⎟⎟ + + bα⎜⎜ + 3φ0 2(1 − ξF3/φ0) ⎝ 1 − ξF3/φ0 ⎠

2⎤ 1 ⎛ ξF3/φ0 ⎞ ⎥ ⎟⎟ + ⎜⎜ 3 ⎝ 1 − ξF3/φ0 ⎠ ⎥⎦

(52)

It is straightforward to check that all the terms after the sixth on the RHS of eq 46 vanish in the case of a one-component fluid. In this case, ξF3 = 4πR31ρ1/3 = η1, a1 = A, and b1 = B. and eq 46 reduces readily to the results we obtained previously for a onecomponent fluid (see eq 42 in ref 6) when the result for βPSPT2a given in eq 42 is substituted into eq 46. Unlike pressure or free energy that is a thermodynamic function for the whole system, the chemical potentials in a multicomponent fluid are specific to each species. This is why there is no isomorphism between the chemical potential of a specific species in the multicomponent fluid with that of a one-component fluid. Nevertheless, it is interesting to note that the expression of the average chemical potential defined by

(44)

Recall that the result given in eq 44 is valid in general although we have chosen a particular integration path to obtain it because dF is a total differential and its integration does not depend on the integration path. Now, we can obtain readily the chemical potential from the Helmholtz free energy

⎤ ⎡ 4πR 3ρ(C − E) 2 ⎛ ξF3/φ0 ⎞ α ⎟⎟ ×⎢ − cα + eα ⎥ + ⎜⎜ 3ξF3 3 ⎝ 1 − ξF3/φ0 ⎠ ⎦ ⎣

∑ x γ eγ γ=1

where Fid is the Helmholtz free energy of a bulk ideal gas mixture.

− ln

(51)

γ=1

⎡ ⎛ ξ ⎞ A ξF3/φ0 + ρ⎢ −ln⎜⎜1 − F3 ⎟⎟ + ⎢ ⎝ φ0 ⎠ 2 1 − ξF3/φ0 ⎣

pIα (R s

∑ xγ c γ

(50)

(55)

The chemical potential of the variant SPT2b is given by 5495


Article

The Journal of Physical Chemistry B ⎛ ⎛ ξ ⎞ ξ ⎞ βμαSPT2b = βμαSPT2a − ln⎜1 − F3 ⎟ + ln⎜⎜1 − F3 ⎟⎟ φ⎠ φ0 ⎠ ⎝ ⎝

⎧ ⎪⎛ βP SPT2d βP SPT2c φ ⎞ ⎟⎟ = + (A + 2B)⎨⎜⎜1 + ⎪ ρ ρ φ − φ0 ⎠ ⎩⎝ ⎤ ⎡φ ⎛ ξF3/φ0 ξ ⎞ 1 ⎥− × ⎢ 0 ln⎜⎜1 − F3 ⎟⎟ + ⎢⎣ ξF3 ⎝ φ0 ⎠ 1 − ξF3/φ0 ⎥⎦ 2(1 − ξF3/φ0)2

⎞ β(P SPT2b − P SPT2a) ⎛ 4πR 3ρ α +⎜ − 1⎟ ρ ⎠ ⎝ 3ξF3 −

⎫ ⎛ φ ⎞2 ⎡ φ ⎛ φ ⎛ ξ ⎞ ξ ⎞⎤ ⎪ ⎟⎟ ⎢ ln⎜1 − F3 ⎟ − 0 ln⎜⎜1 − F3 ⎟⎟⎥⎬ − ⎜⎜ ⎪ φ⎠ ξF3 ⎝ φ0 ⎠⎥⎦⎭ ⎝ φ − φ0 ⎠ ⎢⎣ ξF3 ⎝

⎞⎡ φ ⎛ ⎤ ξ ⎞ 4πR α 3ρ ⎛ φ ⎜⎜ − 1⎟⎟⎢ ln⎜1 − F3 ⎟ + 1⎥ φ⎠ 3ξF3 ⎝ ϕα ⎦ ⎠⎣ ξF3 ⎝

(60)

(56)

Integrating the pressure, we obtain the following expression for

Now, we obtain the pressure for the variant SPT2c by keeping φ in the two first terms on the RHS of eq 41 and taking the limit φ → φ0 in the other terms.

the Helmholtz free energy ⎧ ⎪⎛ βF SPT2d βF SPT2c φ ⎞ ⎟⎟ = − ρ(A + 2B)⎨⎜⎜1 + ⎪ V V φ − φ0 ⎠ ⎩⎝ ⎤ ⎡φ ⎛ ξF3/φ0 ξ ⎞ × ⎢ 0 ln⎜⎜1 − F3 ⎟⎟ + 1⎥ + ⎥⎦ ⎢⎣ ξF3 ⎝ 2(1 − ξF3/φ0) φ0 ⎠

⎧ φ ⎛ ⎪ ξ ⎞ βP SPT2c βP SPT2b = + (1 + A)⎨ − 0 ln⎜⎜1 − F3 ⎟⎟ ⎪ ρ ρ φ0 ⎠ ⎩ ξF3 ⎝ ⎡ ξ ⎞ φ ⎢φ ⎛ 1 ln⎜1 − F3 ⎟ − + 1 − ξF3/φ0 φ − φ0 ⎢⎣ ξF3 ⎝ φ⎠ ⎫ φ ⎛ ξ ⎞⎤ ⎪ − 0 ln⎜⎜1 − F3 ⎟⎟⎥⎬ ⎪ ξF3 ⎝ φ0 ⎠⎥⎦⎭

⎛ φ ⎞ 2 ⎡⎛ ξ ⎞ φ⎞ ⎛ ⎟⎟ ⎢⎜1 − + ⎜⎜ ⎟ ln⎜1 − F3 ⎟ ξF3 ⎠ ⎝ φ⎠ ⎝ φ − φ0 ⎠ ⎢⎣⎝ ⎫ ⎛ φ ⎞ ⎛ ξ ⎞⎤⎪ − ⎜1 − 0 ⎟ ln⎜⎜1 − F3 ⎟⎟⎥⎬ ⎪ ξF3 ⎠ ⎝ φ0 ⎠⎥⎦⎭ ⎝

(57)

The Helmholtz free energy is given by β FSPT2c β FSPT2b = V V

⎧ φ ⎛ ⎪ ξ ⎞ 1 + 0 ln⎜⎜1 − F3 ⎟⎟ + ρ(1 + A)⎨ ⎪ ξF3 ⎝ φ0 ⎠ ⎩

The chemical potential is given by ⎧ (2R α)2 3ξ ξ ⎡ 2R βμαSPT2d = βμαSPT2c − ⎨aα + 2bα + F1 F2 ⎢ α + ξF3 ⎣ ξF1 ξF2 ⎩ 3⎤ 3 2 2⎛ (2R α) (2R α) ⎤⎫ 9ξ 4ξ ⎞⎡ (2R α) ⎥⎬ ⎥ + F2 ⎜D1 − F2 ⎟⎢ −2 − ξF3 ⎦ 2ξF3 ⎝ ξF3 ⎠⎣ ξF3 ξF2 ⎦⎭ ⎧ ⎤ ⎡ ⎪⎛ ξF3/φ0 ξ ⎞ φ ⎞⎢ φ0 ⎛ ⎟⎟ ln⎜⎜1 − F3 ⎟⎟ + 1⎥ + × ⎨⎜⎜1 + ⎪ ⎥ ⎢ 2(1 φ − φ ξ φ − ξF3/φ0) ⎠ ⎠ ⎝ ⎝ ⎦ 0 ⎣ F3 0 ⎩

⎡ φ ⎞ ξ ⎞ ⎛ φ ⎢⎛ φ⎞ ⎛ + ⎜1 − ⎟ ln⎜1 − F3 ⎟ − ⎜1 − 0 ⎟ φ − φ0 ⎢⎣⎝ ξF3 ⎠ ⎝ φ⎠ ⎝ ξF3 ⎠ ⎫ ⎛ ξ ⎞⎤⎪ × ln⎜⎜1 − F3 ⎟⎟⎥⎬ ⎪ φ0 ⎠⎥⎦⎭ ⎝

⎪

⎪

⎪

⎪

(58)

The chemical potential has the following expression

⎛ φ ⎞2 ⎡⎛ φ ⎞ ξ ⎞ ⎛ φ⎞ ⎛ ⎟⎟ ⎢⎜1 − + ⎜⎜ ⎟ ln⎜1 − F3 ⎟ − ⎜1 − 0 ⎟ φ⎠ ⎝ ξF3 ⎠ ξF3 ⎠ ⎝ ⎝ φ − φ0 ⎠ ⎢⎣⎝

⎧ (2R α)2 3ξ ξ ⎡ 2R βμαSPT2c = βμαSPT2b + ⎨1 + aα + F1 F2 ⎢ α + ξF3 ⎣ ξF1 ξF2 ⎩ (2R α)3 ⎤ 3D1ξF2 2 ⎡ (2R α)3 (2R α)2 ⎤⎫ ⎥+ ⎢ ⎥⎬ −2 − 2ξF3 ⎣ ξF3 ξF3 ⎦ ξF2 ⎦⎭ ⎧ ⎡ ⎪ φ ⎛ ξ ⎞ ξ ⎞ φ ⎢⎛ φ⎞ ⎛ × ⎨1 + 0 ln⎜⎜1 − F3 ⎟⎟ + ⎜1 − ⎟ ln⎜1 − F3 ⎟ ⎪ ⎢ ξ φ φ φ ξ φ⎠ − ⎝ ⎝ F3 F3 ⎠ 0⎠ 0 ⎣⎝ ⎩ ⎪

⎪

⎫ ⎛ ξ ⎞⎤⎪ 4πR α 3ρ β(P SPT2d − P SPT2c) × ln⎜⎜1 − F3 ⎟⎟⎥⎬ + ⎪ 3ξF3 φ0 ⎠⎥⎦⎭ ρ ⎝

⎪

⎪

⎧ φ ⎞ φ ⎪ 4πR α 3ρ ⎛ φ 0 ⎜⎜ − 1⎟⎟ ⎨ 3ξF3 ⎝ ϕα ⎠ φ − φ0 ⎪ ⎩ φ − φ0 ⎤ ⎡ϕ ⎛ ⎤ ξ ⎞ ξ ⎞ φ ⎡φ ⎛ ⎢ ln⎜1 − F3 ⎟ + 1⎥ × ⎢ 0 ln⎜⎜1 − F3 ⎟⎟ + 1⎥ + ⎥⎦ ⎢⎣ ξF3 ⎝ φ0 ⎠ φ − φ0 ⎣ ξF3 ⎝ φ⎠ ⎦

+ (A + 2B)

⎫ ⎛ φ ⎞ ⎛ ξ ⎞⎤⎪ 4πR α 3ρ β(P SPT2c − P SPT2b) − ⎜1 − 0 ⎟ ln⎜⎜1 − F3 ⎟⎟⎥⎬ + ⎪ 3ξF3 ξF3 ⎠ ⎝ ϕ0 ⎠⎥⎦⎭ ρ ⎝ ⎞ φ ⎧ ⎪ ξ ⎞ 4πR α 3ρ ⎛ φ ϕ ⎛ ⎜⎜ − 1⎟⎟ ⎨1 + ln⎜1 − F3 ⎟ ⎪ 3ξF3 ⎝ ϕα ξF3 ⎝ φ⎠ ⎠ φ − φ0 ⎩ ⎞ ⎛ φ0 ⎡⎛ φ ⎞ ξ ⎞ ⎛ ⎢⎜1 − φ ⎟ ln⎜1 − F3 ⎟ − ⎜1 − 0 ⎟ + φ − φ0 ⎢⎣⎝ ξF3 ⎠ ⎝ φ⎠ ⎝ ξF3 ⎠ + (1 + A)

⎫ ⎛ ξ ⎞⎤⎪ × ln⎜⎜1 − F3 ⎟⎟⎥⎬ ⎪ φ0 ⎠⎥⎦⎭ ⎝

(61)

+

⎡⎛ ⎞ ⎛ φ ⎞ ξ ⎞ ⎛ ⎢⎜1 − φ ⎟ ln⎜1 − F3 ⎟ − ⎜1 − 0 ⎟ ξF3 ⎠ ⎝ φ⎠ ⎝ ξF3 ⎠ (φ − φ0) ⎢⎣⎝ 2φ0φ

2

⎫ ⎛ ξ ⎞⎤⎪ × ln⎜⎜1 − F3 ⎟⎟⎥⎬ ⎪ φ0 ⎠⎥⎦⎭ ⎝

(62)

Now, we proceed to calculating the Helmholtz free energy of the (59)

original version of SPT2, i.e., integrating the pressure given in eq

To obtain SPT2d, we take the limit φ → φ0 in the fourth term on the RHS of eq 41, and this yields

41 without any modification on it. To facilitate this calculation, we first rewrite eq 41 as 5496


Article

The Journal of Physical Chemistry B ⎧ ⎡ ⎛ φ ⎞2 ⎤ ⎪ βP SPT2 βP SPT2d φ = + 2B⎨ − ⎢1 + +⎜ ⎟⎥ ⎪ ⎢ ρ ρ φ − φ0 ⎝ φ − φ 0 ⎠ ⎥⎦ ⎩ ⎣ ⎤ ⎛ ⎡φ ⎛ ξ ⎞ φ ⎞ 1 ⎥ + ⎜⎜1 + ⎟⎟ × ⎢ 0 ln⎜⎜1 − F3 ⎟⎟ + ⎢⎣ ξF3 ⎝ φ0 ⎠ 1 − ξF3/φ0 ⎥⎦ ⎝ φ − φ0 ⎠ ⎛ φ ⎞3 ⎟⎟ ⎜⎜ + 2(1 − ξF3/φ0)2 3(1 − ξF3/φ0)3 ⎝ φ − φ0 ⎠ ⎫ ⎡ φ ⎛ ξ ⎞ ξ ⎞⎤⎪ φ ⎛ ln⎜1 − F3 ⎟ − 0 ln⎜⎜1 − F3 ⎟⎟⎥⎬ ×⎢ ⎪ ⎢⎣ ξF3 ⎝ φ⎠ ξF3 ⎝ φ0 ⎠⎥⎦⎭ ξF3/φ0

×

2.2. Simulation. To make comparisons between the results given by different SPT variants and simulation ones, we carried out some GCEMC simulations of a binary HS mixture in a onecomponent HS matrix or a one-component OHS matrix. Species 0 denotes the matrix species, whereas species 1 and species 2 are used to designate the two fluid species. We will use the diameter of the first fluid species as unit of length (i.e., σ1 = 2R1 = σ). Matrices with different particle sizes (R0 = 0.5σ and 0.9σ) as well as different porosities were considered. A cubic simulation box with the edge length equal to 10σ is used with periodic boundary conditions. Each simulation run is performed with a particular matrix realization, and the fluid particles can move in the volume that is not occupied by matrix particles. The matrix configurations are generated by using Monte Carlo simulations in a canonical ensemble. Because matrices of a finite size are used in the simulation, any observable quantity fluctuates with matrix realizations; hence, an average over matrix realizations should be also taken. In the present work, we use 25 realizations to make such averages. For each matrix realization, about 9 × 105 trial moves per particle are performed during which the average fluid density is calculated. Here, “trial move” has the general meaning of the displacement of particle position, particle insertion, or removal. In our simulations, an insertion or removal trial is attempted after each Monte Carlo cycle, i.e., one attempted displacement for every particle. The primary purpose of the simulations performed in the present work is to provide numerical experimental data in order to assess the accuracy of our theory. For this purpose, the different parameters (matrix particle size and fluid and matrix densities) were chosen to cover a range of situations. Moreover, the considered matrices are theoretically tractable models. Usually, real porous materials have more complex morphological and geometrical structures. However, porosity and specific pore size are common characteristics of porous materials. It is noteworthy that even the simple models considered in the present work allow one to adjust the model parameters (matrix particle size and density) to generate matrix samples with specified characteristics equivalent to those of the real porous materials.

(ξF3/φ0)2

−

(63)

Integrating the pressure given in the above expression leads to the following expression for the Helmholtz free energy ⎧⎡ ⎛ φ ⎞2 ⎤ ⎪⎢ βF SPT2 βF SPT2d φ ⎟⎟ ⎥ = + 2Bρ⎨ 1 + + ⎜⎜ V φ − φ0 φ − φ0 ⎠ ⎥⎦ V ⎪⎢⎣ ⎝ ⎩ ⎤ ⎛ ⎡φ ⎛ ⎞ ξF3/φ0 ξ φ ⎞ ⎟⎟ × ⎢ 0 ln⎜⎜1 − F3 ⎟⎟ + 1⎥ + ⎜⎜1 + ⎥⎦ ⎝ ⎢⎣ ξF3 ⎝ φ0 ⎠ φ − φ0 ⎠ 2(1 − ξF3/φ0) −

2 ⎛ φ ⎞3⎡⎛ ξ ⎞ φ⎞ ⎛ 1 ⎛ ξF3/φ0 ⎞ ⎟⎟ ⎢⎜1 − ⎟⎟ + ⎜⎜ ⎜⎜ ⎟ ln⎜1 − F3 ⎟ ξF3 ⎠ ⎝ φ⎠ 6 ⎝ 1 − ξF3/φ0 ⎠ ⎝ φ − φ0 ⎠ ⎢⎣⎝

⎫ ⎛ φ ⎞ ⎛ ξ ⎞⎤⎪ − ⎜1 − 0 ⎟ ln⎜⎜1 − F3 ⎟⎟⎥⎬ ⎪ φ0 ⎠⎥⎦⎭ ξF3 ⎠ ⎝ ⎝

(64)

We obtain finally the following expression for chemical potential ⎧ 3ξ 2 ⎛ 6ξ ⎞ βμαSPT2 = βμαSPT2d + ⎨2bα + F2 ⎜D1 − F2 ⎟ ξF3 ⎝ ξF3 ⎠ ⎩ ⎧⎡ ⎛ φ ⎞2 ⎤ ⎡ (2R )3 (2R α)2 ⎤⎫⎪⎢ φ α ⎟⎟ ⎥ ⎥⎬⎨ 1 + ×⎢ − + ⎜⎜ ξF2 ⎦⎭⎪⎢⎣ φ − φ0 φ − φ0 ⎠ ⎥⎦ ⎣ ξF3 ⎝ ⎩ ⎪

⎪

⎪

⎪

⎤ ⎛ ⎡φ ⎛ ξF3/φ0 ξ ⎞ φ ⎞ ⎟⎟ ×⎢ 0 ln⎜⎜1 − F3 ⎟⎟ + 1⎥ + ⎜⎜1 + ⎢⎣ ξF3 ⎝ φ0 ⎠ φ − φ0 ⎠ 2(1 − ξF3/φ0) ⎦⎥ ⎝

3. RESULTS AND DISCUSSION 3.1. Porosity. From the presentation of the formulation of SPT2 given in section 2, we see that different porosities appear in SPT2, i.e., φ0 (geometric one) and ϕ1 and ϕ2 (probe-specific porosities for the two fluid species considered here). It is clear that the accuracy of the results for the chemical potentials and pressure depends on the accuracy with which these porosities can be calculated. First, it is noteworthy that in the case of an OHS matrix probed by a HS particle both φ0 and ϕα (α = 1, 2) can be determined exactly. Thus, we have

2 ⎛ φ ⎞3⎡⎛ ξ ⎞ φ⎞ ⎛ 1 ⎛ ξF3/φ0 ⎞ ⎟⎟ ⎢⎜1 − ⎟⎟ + ⎜⎜ − ⎜⎜ ⎟ ln⎜1 − F3 ⎟ 6 ⎝ 1 − ξF3/φ0 ⎠ ξF3 ⎠ ⎝ φ⎠ ⎝ φ − φ0 ⎠ ⎢⎣⎝

⎫ ⎛ φ ⎞ ⎛ ξ ⎞⎤⎪ 4πR α 3ρ β(P SPT2 − P SPT2d) − ⎜1 − 0 ⎟ ln⎜⎜1 − F3 ⎟⎟⎥⎬ + ⎪ 3ξF3 ξF3 ⎠ ⎝ φ0 ⎠⎥⎦⎭ ρ ⎝ + 2B

⎞ φ ⎧ ⎪ φ0 ⎛ 4πR α 3ρ ⎛ φ 2φ ⎞ ⎟⎟ ⎜⎜1 + ⎜⎜ − 1⎟⎟ ⎨ ⎪ 3ξF3 ⎝ ϕα φ − φ0 ⎠ ⎠ φ − φ0 ⎩ φ − φ0 ⎝

2 ⎤ ⎛ ⎡φ ⎛ ξF3 ⎞ ξ ⎞ φ ⎞⎡φ ⎛ 0 ⎥ ⎢ ⎟ ⎟ ⎜ ⎜ × ln⎜1 − ln⎜1 − F3 ⎟ ⎟⎢ ⎟+1 +⎜ ⎥⎦ ⎝ φ − φ0 ⎠ ⎣ ξF3 ⎝ ⎢⎣ ξF3 ⎝ φ0 ⎠ φ⎠

3

ϕαOHS = e−4π(R 0+ R α) ρ0 /3

(66)

and φ0 can be obtained by substituting Rα = 0 into eq 66. In the case of a HS matrix, only φ0 can be calculated exactly as

⎤ 3φ0φ2 φ0 ξF3/φ0 + 1⎥ + + × φ − φ0 2(1 − ξF3/φ0) (φ − φ0)3 ⎦ ⎫ ⎡⎛ ⎞ ⎛ φ ⎞ ⎛ ξ ⎞ ⎛ ξ ⎞⎤⎪ ⎢⎜1 − φ ⎟ ln⎜1 − F3 ⎟ − ⎜1 − 0 ⎟ ln⎜⎜1 − F3 ⎟⎟⎥⎬ ⎪ ⎢⎣⎝ ξF3 ⎠ ⎝ φ⎠ ⎝ ξF3 ⎠ ⎝ φ0 ⎠⎥⎦⎭

α = 1, 2

φ0HS = 1 − 4πR 0 3ρ0 /3

(67)

For probing particle of finite size, we can use SPT to calculate ϕHS α

(65)

ϕαHS = pIα (R s = R α)

The consistency of the results for all the variants given above has been checked, and a more detailed description is provided in Appendix B. The numerical performance of all these variants will be accessed against Monte Carlo simulation results, as discussed in section 3.

α = 1, 2

(68)

pαI (Rs

with = Rα) being calculated by applying eqs 15−23 to the particular case of m = 1 because we consider here a onecomponent HS matrix. Thus, it is natural to assess first the accuracy of SPT for ϕHS α (α = 1, 2). For this, we also carried out 5497


Article

The Journal of Physical Chemistry B

ex Figure 1. Relative errors (Error α = [βμex α (SPT) − βμα (GCEMC)]/ βμex α (GCEMC), α = 1, 2) of different SPT2 variants (SPT2a: black curve; SPT2b: red curve; SPT2c: green curve; SPT2d: blue curve; and SPT2: pink curve) with respect to the simulation results given in Table S1 for a binary HS mixture (R1 = 0.5σ, R2 = 0.4σ) in a one-component HS matrix under different conditions: (a) ρ0σ3 = 0.1, R0 = R1; (b) ρ0σ3 = 0.3, R0 = R1; (c) ρ0σ3 = 0.017, R0 = 1.8R1; and (d) ρ0σ3 = 0.051, R0 = 1.8R1. Results in the left column are for the chemical potential of species 1 (big fluid sphere), and those on the right column for the chemical potential of species 2 (small fluid sphere).

ex Figure 2. Relative errors (Error α = [βμex α (SPT) − βμα (GCEMC)]/ (GCEMC), α = 1, 2) of different SPT2 variants (SPT2a: black βμex α curve; SPT2b: red curve; SPT2c: green curve; SPT2d: blue curve; and SPT2: pink curve) with respect to the simulation results given in Table S2 for a binary HS mixture (R1 = 0.5σ, R2 = 0.4σ) in a one-component HS matrix under different conditions: (a) ρ0σ3 = 0.1, R0 = R1; (b) ρ0σ3 = 0.3, R0 = R1; (c) ρ0σ3 = 0.017, R0 = 1.8R1; and (d) ρ0σ3 = 0.051, R0 = 1.8R1. Results in the left column are for the chemical potential of species 1 (big fluid sphere), and those on the right column for the chemical potential of species 2 (small fluid sphere).

Monte Carlo simulations to determine ϕHS α . A probing HS particle of radius Rα (α = 1, 2) is inserted randomly into a matrix realization (with about 107 attempted insertions), and the successful insertion probability (i.e., rate of successful insertions) can be easily calculated, which is equal to ϕHS α (α = 1, 2). The simulation results along with the SPT ones are presented in Table 1 for HS matrices. The statistical error of the simulation results given in Table 1 is less than 0.5%. The results in Table 1 show that SPT gives very accurate results for ϕHS α (α = 1, 2) under the considered conditions with an error less than 0.2% (smaller than the statistical errors of simulations under the considered conditions). The matrix density and the size of matrix particles are chosen in such a way that we have two groups of matrices. The matrices in the same group have the same geometric porosity, i.e., in Table 1 the matrices in the second and the third columns with φ0 = 0.95 and the matrices in the fourth and fifth columns with φ0 = 0.84. From Table 1, we see that the porosity determined by using a probing particle of finite size can be significantly lower than the geometric porosity, in particular when the matrix density is high and the probing particle is large. Comparing the results in the second and third columns (or those in the fourth and fifth columns) in Table 1, we see that at the same geometric porosity the porosity determined by using a

particle of finite size is increased when the size of matrix particle is larger, and this implies that more fluid particles can be accommodated. 3.2. Chemical Potential. To assess the accuracy of our different SPT variants, we carried out several series of GCEMC simulations for a binary HS fluid in a one-component HS matrix or in a one-component OHS matrix. In each series, the chemical potential of one species is kept at a fixed value, whereas the chemical potential of the other species is varied systematically. The simulation results are presented in the Supporting Information along with the numerical results obtained from the different variants of SPT2. The statistical errors of our simulation results do not exceed 1.5%. From the simulation results, one can easily see the following general trend. For given chemical potentials, the total fluid density decreases with the geometric porosity. For the same geometric porosity, the fluid density is always higher in the matrices with larger matrix particles. This conforms perfectly to the prediction we can make from the results for the porosity determined by using a probing particle of finite size (see the discussions given at the end of section 3.1). It is also noteworthy that for given geometric porosity and chemical potentials that the 5498


Article


Table 1. Porosities of HS Matrices Calculated by Using Respectively Species 1 (R1 = 0.5σ) or Species 2 (R2 = 0.4σ) as Probe φ0 HS MC: ϕHS 1 ; ϕ2 HS SPT: ϕHS ; ϕ 1 2

R0 = R1 ρ0σ3 = 0.1

R0 = 1.8R1 ρ0σ3 = 0.017

R0 = R1 ρ0σ3 = 0.3

R0 = 1.8R1 ρ0σ3 = 0.051

0.95 0.6279; 0.7153 0.6287; 0.7156

0.95 0.8104; 0.8462 0.8102; 0.8461

0.84 0.1747; 0.2925 0.1728; 0.2909

0.84 0.4737; 0.5592 0.4737; 0.5590

ex Figure 3. Relative errors (Error α = [βμex α (SPT) − βμα (GCEMC)]/ ex βμα (GCEMC), α = 1, 2) of different SPT2 variants (SPT2a: black curve; SPT2b: red curve; SPT2c: green curve; SPT2d: blue curve; and SPT2: pink curve) with respect to the simulation results given in Table S3 for a binary HS mixture (R1 = 0.5σ, R2 = 0.4σ) in a one-component OHS matrix under different conditions: (a) ρ0σ3 = 0.1, R0 = R1; (b) ρ0σ3 = 0.3, R0 = R1; (c) ρ0σ3 = 0.017, R0 = 1.8R1; and (d) ρ0σ3 = 0.051, R0 = 1.8R1. Results in the left column are for the chemical potential of species 1 (big fluid sphere), and those on the right column for the chemical potential of species 2 (small fluid sphere).

ex Figure 4. Relative errors (Error α = [βμex α (SPT) − βμα (GCEMC)]/ βμex α (GCEMC), α = 1, 2) of different SPT2 variants (SPT2a: black curve; SPT2b: red curve; SPT2c: green curve; SPT2d: blue curve; and SPT2: pink curve) with respect to the simulation results given in Table S4 for a binary HS mixture (R1 = 0.5σ, R2 = 0.4σ) in a one-component OHS matrix under different conditions: (a) ρ0σ3 = 0.1, R0 = R1; (b) ρ0σ3 = 0.3, R0 = R1; (c) ρ0σ3 = 0.017, R0 = 1.8R1; and (d) ρ0σ3 = 0.051, R0 = 1.8R1. Results in the left column are for the chemical potential of species 1 (big fluid sphere), and those on the right column for the chemical potential of species 2 (small fluid sphere).

fluid concentration is nearly the same in matrices with different particle sizes. To assess the accuracy of SPT2 and its different variants described in section 2.1, we present on one hand the numerical results of these theoretical approximations with the simulation ones in form of tables given in Supporting Information and show on the other hand the relative errors of SPT with respect to simulation results in the form of the figures given below. Figures 1 and 2 show, respectively, the comparison of different SPT approximations with the simulation results for a binary HS mixture in a one-component HS matrix presented in Tables S1 and S2. Figures 3 and 4 show, respectively, the comparison for a one-component OHS matrix with the simulation results given in Tables S3 and S4. The relative errors of different SPT variants for 3 excess chemical potentials, defined as βμex α = βμα − ln(ρασ ), are plotted in Figures 1−4 as a function of fluid density. As for a one-

component fluid,6 the values of chemical potential given by different approximations become larger and larger following the order SPT2a, SPT2b, SPT2c, SPT2d, and SPT2. For the onecomponent fluid, it was found that SPT2b gives the most accurate results.6 This holds in most cases for a binary mixture. Nevertheless, we see from Figures 1(b-1) and 1(b-2) that in the case of a HS matrix with ρ0σ3 = 0.3 and R0 = R1 SPT2a gives more accurate results than SPT2b in particular at high fluid densities. In most cases, the accuracy of SPT is similar for both species. We observe also some situations where the chemical potential of the large species (Figures 1(b-1) and 2(b-1)) is described less accurately by SPT than that of the small species (Figures 1(b-2) and 2(b-2)). As for the one-component fluid,6 the accuracy of different SPT2 approximations is higher for OHS matrices than for HS matrices. 5499


Article


ϕ ⎛ 1 − η1/ϕ ⎞ ϕ ⎟⎟ = ln⎜⎜ ln(1 − η1/ϕ) η1 ⎝ 1 − η1/ϕ0 ⎠ η1 ϕ ϕ − ϕ0 ϕ − 0 ln(1 − η1/ϕ0) − ln(1 − η1/ϕ0) η1 ϕ η1

4. CONCLUDING REMARKS In the present work, we extended the SPT to a quite general model for a multicomponent fluid in a multicomponent matrix. We find that there is a remarkable isomorphism between such a multicomponent system and the system of a one-component fluid in a one-component matrix for some thermodynamics properties such as pressure and Helmholtz free energy. SPT for multicomponent systems naturally takes the structure of a onefluid theory. It is noteworthy that the mixing rule for porosity has a quite appealing form with the ratio of packing fractions as weight rather than the species concentrations. The comparison between the theory and simulation shows that in most cases the approximation SPT2b has the best accuracy. In contrast to the case of a one-component fluid in a one-component matrix, it was found that under some conditions (in particular in matrices of low porosity and at high fluid densities) SPT2a gives more accurate results. In almost all cases, the accuracy of the different approximations for the chemical potential of the small species is better than that for the large species.

■

(A3)

ϕ ϕ ⎛ 1 − η1/ϕ ⎞ ϕ ⎟⎟ = ln⎜⎜ ln(1 − η1/ϕ) − 0 ln(1 − η1/ϕ0) η1 ⎝ 1 − η1/ϕ0 ⎠ η1 η1 ⎡⎛ ϕ − ϕ ⎞ 2 ϕ − ϕ0 ϕ0 ⎤ ϕ 0 ⎥ ln(1 − η /ϕ ) − ⎢⎜ + ⎟ 1 0 ⎢⎣⎝ ϕ ⎠ η1 ϕ η1 ⎥⎦ (A4)

Equations A3 and A4 can be obtained from eq A2 by using the following identity, η1/ϕ0 ϕ = ϕ − ϕ0 η1/ϕ0 − η1/ϕ

Substituting eqs A2−A4 into the second, third, and fourth terms on the RHS of eq A1 respectively, we obtain

APPENDIX A

Isomorphism between Multicomponent and One-Component Systems

(1 + a)ϕ ⎡ ϕ βP ϕ ⎢ ln(1 − η1/ϕ) = − ln(1 − η1/ϕ) + (ϕ − ϕ0) ⎢⎣ η1 ρ η1

In this appendix, we discuss, in more details, the isomorphism between the SPT results for a multicomponent system and those for a one-component system. In ref 6, we obtained the following expression for the pressure (eq 39 in ref 6) in the case of a one-component fluid: βP ϕ ⎡ 1 − η1/ϕ ⎤ ϕ ϕ ⎥ + (1 + a) = − ln⎢ ρ η1 ⎢⎣ 1 − η1/ϕ0 ⎥⎦ η1 (ϕ − ϕ0)

⎧ ⎫ ⎡ 1 − η / ϕ ⎤⎪ ϕ ϕ 1 1 ⎢ ⎥⎬ ⎨ − ln ⎪ η1 (ϕ − ϕ0) ⎢⎣ 1 − η1/ϕ0 ⎥⎦⎪ ⎩ 1 − η1/ϕ0 ⎭ ⎧ η1/ϕ0 2ϕ − ϕ0 ϕ ⎪ 1 ⎨ − 2 ⎪ ϕ − ϕ0 ⎩ 2(1 − η1/ϕ0) ϕ − ϕ0 1 − η1/ϕ0

(A1)

ln(1 − η1/ϕ0)

η1

ln(1 − η1/ϕ0) −

ϕ ⎡ϕ ⎢ ln(1 − η1/ϕ) ϕ − ϕ0 ⎢⎣ η1

ϕ0

⎪

⎪

(A6)

Now, it is clear that eq A6 is isomorphic with eq 41 since one obtains eq 41 by substituting the various quantities of the onecomponent fluid in eq A6 by the corresponding averaged quantities of the multicomponent fluid, i.e., a → A, b → B, η1 → ξF3, and ϕ → φ. It is noteworthy that ϕ0 = φ0. This isomorphism shows that SPT can be considered in fact as a one-fluid theory. The mixing rule for parameters, A and B has the usual form [see eqs(38 and (39)], i.e., with the concentrations (xα, α = 1,...,n) as weights. ξF3 can be considered as the mixing rule for the particle volumes of different species with the number density as weight.

ϕ ⎛ 1 − η1/ϕ ⎞ ϕ ⎟⎟ = ln⎜⎜ ln(1 − η1/ϕ) η1 ⎝ 1 − η1/ϕ0 ⎠ η1 ϕ ϕ − 0 ln(1 − η1/ϕ0) − ln(1 − η1/ϕ0) η1 η1 η1

ϕ0

⎫ ⎡ϕ ⎤⎪ ϕ0 ⎥ ⎢ × ln(1 − η1/ϕ) − ln(1 − η1/ϕ0) ⎬ ⎢⎣ η1 ⎥⎦⎪ η1 ⎭

In eq A1, a slight modification of notation, with respect to the original eq 39 in ref 6, is made and we write now the two parameters with lowercase letters, i.e., a and b, because we reserve the capital letters for the averaged parameters in the case of a multicomponent fluid. The direct comparison of eq A1 with eq 41 does not show the isomorphism straightforwardly. Some algebraic rearrangements are necessary to show this by using the following identities,

ϕ0

+

⎡ ⎤ ⎛ ϕ ⎞2 ϕ0 1 ⎟ ×⎢ + ln(1 − η1/ϕ0)⎥ +⎜⎜ ⎢⎣ 1 − η1/ϕ0 ⎥⎦ ⎝ ϕ − ϕ0 ⎟⎠ η1

2

+

ϕ0

−

⎪

⎫ ⎡ 1 − η / ϕ ⎤⎪ ϕ ϕ 1 ⎥⎬ + ln⎢ 2 η1 (ϕ − ϕ0) ⎢⎣ 1 − η1/ϕ0 ⎥⎦⎪ ⎭

⎧ ⎤ (a + 2b)ϕ ⎪ 1 ⎨ ln(1 − η1/ϕ0)⎥ + ⎥⎦ − η1 ϕ − ϕ0 ⎪ η1/ϕ0 1 ⎩

−

⎫ ⎤⎪ η1/ϕ0 2bϕ ⎧ ⎬+ ⎨ ln(1 − η1/ϕ0)⎥⎪ ⎥⎦⎭ η1 ϕ − ϕ0 ⎩ 2(1 − η1/ϕ0)2 ϕ ϕ 1 − − 0 ln(1 − η1/ϕ0) − η1 ϕ − ϕ0 1 − η1/ϕ0

⎡ 1 − η /ϕ ⎤ ϕ 1 ⎥ + (a + 2b) × ln⎢ × ⎢⎣ 1 − η1/ϕ0 ⎥⎦ ϕ − ϕ0

+ 2b

(A5)

n

ξF3 =

∑ ρλ λ=1

4πRλ 3 3

(A7)

It is also quite appealing that the weight in the mixing rule for porosity resulting from SPT is the ratio of packing fraction, i.e., ηα/ξF3, rather than the species concentration [see eq 40]. It is noteworthy also that the mixing rule concerning the porosity is

(A2) 5500


Article

The Journal of Physical Chemistry B expressed in terms of φ−1 and ϕ−1 α (α = 1, ..., n) rather than φ and ϕα (α = 1, ..., n). We would also like to point out that the one-fluid theory structure of SPT holds for both confined and bulk multicomponent fluids. It does not seem that this has been noted and received any attention until now. The expression for pressure we gave in eq 20 for a bulk multicomponent HS fluid is different from that in the original paper of Lebowitz, Helfand and Praestgaard.3 The purpose of recasting it in such a form is to highlight the striking similarity of the general structure of SPT for all the cases (bulk one-component and multicomponent fluids as well as confined one-component and multicomponent fluids). In the following sections, we will show the isomorphism between one-component and multicomponent fluids for all the SPT2 variants. Although the expressions for Helmholtz free energy were not worked out in ref 6 for one-component systems, we checked that isomorphism holds for free energy as well. The self-contained expressions for Helmholtz free energy will be given below (the expressions of free energy given in section 2.1 are more compact but expressed in term of another variant except for SPT2a). SPT2a.

ξ ⎞ (A + 1)φ ⎡ φ ⎛ ξ ⎞ βP SPT2c φ ⎛ ⎢ ln⎜1 − F3 ⎟ + ln⎜1 − F3 ⎟ =− ρ ξF3 ⎝ φ⎠ φ − φ0 ⎢⎣ ξF3 ⎝ φ⎠ −

−

Helmholtz free energy is βF SPT2c βF id = − V V

n

∑ ργ ln pIγ (R s = R γ) γ=1

2 ⎧ ⎛ ⎪ ξ ⎞ A ξF3/φ0 B ⎛ ξF3/φ0 ⎞ ⎟⎟ + ρ⎨ − ln⎜⎜1 − F3 ⎟⎟ + + ⎜⎜ ⎪ 2 1 − ξF3/φ0 3 ⎝ 1 − ξF3/φ0 ⎠ φ0 ⎠ ⎩ ⎝ ⎡ φ ⎛ ξ ⎞⎤ + (1 + A)⎢1 + 0 ln⎜⎜1 − F3 ⎟⎟⎥ ⎢⎣ ξF3 ⎝ φ0 ⎠⎥⎦

⎡ ⎡ ξ ⎞ φ ⎤⎢⎛ φ⎞ ⎛ ⎥ ⎜1 − + ⎢ − 1 + (1 + A) ⎟ ln⎜1 − F3 ⎟ ⎢⎣ ⎥ φ − φ0 ⎦⎢⎣⎝ ξF3 ⎠ ⎝ φ⎠ ⎫ ⎛ φ ⎞ ⎛ ξ ⎞⎤⎪ − ⎜1 − 0 ⎟ ln⎜⎜1 − F3 ⎟⎟⎥⎬ ⎪ ξF3 ⎠ ⎝ φ0 ⎠⎥⎦⎭ ⎝

(A12)

SPT2d. Substituting eq A11 into eq 60, we obtain the following more explicit expression for pressure given by SPT2d,

2 ξF3/φ0 1 A 2B (ξF3/φ0) + + 1 − ξF3/φ0 2 (1 − ξF3/φ0)2 3 (1 − ξF3/φ0)3

ξ ⎞ (A + 1)φ ⎡ φ ⎛ ξ ⎞ βP SPT2d φ ⎛ ⎢ ln⎜1 − F3 ⎟ + ln⎜1 − F3 ⎟ =− ρ ξF3 ⎝ φ⎠ φ − φ0 ⎢⎣ ξF3 ⎝ φ⎠

(A9)

⎧ ⎛ ξ ⎞⎤ φ ⎪ 1 ⎨ ln⎜⎜1 − F3 ⎟⎟⎥ + (A + 2B) ⎪ 1 − ξ /φ ξF3 ⎝ φ0 ⎠⎥⎦ φ − φ0 ⎩ F3 0 ⎡ φ ⎛ ξ ⎞ ξ ⎞ φ ⎢φ ⎛ ln⎜1 − F3 ⎟ + 0 ln⎜⎜1 − F3 ⎟⎟ − ξF3 ⎝ φ0 ⎠ φ − φ0 ⎢⎣ ξF3 ⎝ φ⎠

One sees straightforwardly that eq A9 is isomorphic to eq 46 in ref 6. The self-contained expression for Helmholtz free energy is n

∑ ργ ln pIγ (R s = R γ)

−

φ0

−

φ0

γ=1

⎡ ⎛ ξ ⎞ A ξF3/φ0 + ρ⎢ −ln⎜⎜1 − F3 ⎟⎟ + ⎢ ⎝ φ0 ⎠ 2 1 − ξF3/φ0 ⎣

+

2 ⎛ ξ ⎞ φ⎞ ⎛ B ⎛ ξF3/φ0 ⎞ ⎜⎜ ⎟⎟ − ⎜1 − ⎟ ln⎜1 − F3 ⎟ 3 ⎝ 1 − ξF3/φ0 ⎠ ξF3 ⎠ ⎝ φ⎠ ⎝

⎤ ⎛ φ0 ⎞ ⎛ ξF3 ⎞⎥ ⎟⎟ + ⎜1 − ⎟ ln⎜⎜1 − ξF3 ⎠ ⎝ φ0 ⎠⎥⎦ ⎝

(A11)

corrected in an erratum). The self-contained expression for

φ ⎛ ξ ⎞ ξ ⎞ βP SPT2b φ ⎛ =− ln⎜1 − F3 ⎟ + 0 ln⎜⎜1 − F3 ⎟⎟ ρ ξF3 ⎝ φ⎠ ξF3 ⎝ φ0 ⎠

+

⎤ (ξF3/φ0)2 ⎥ + 2B 3 (1 − ξF3/φ0)3 2(1 − ξF3/φ0)2 ⎥⎦ ξF3/φ0

seen immediately (note that a misprinting in ref 6, eq 48, was

The comparison of eq A8 to eq 60 in ref 6 shows immediately their isomorphism. The expression for Helmholtz free energy given in eq 44 is already a self-contained one. SPT2b. Substituting eq 42 into eq 54, we obtain the following more explicit expression for the pressure given by SPT2b.

βF SPT2b βF id = − V V

ξF3

⎡φ ⎛ ⎛ ξ ⎞⎤ ξ ⎞ 1 ln⎜⎜1 − F3 ⎟⎟⎥ − A⎢ 0 ln⎜⎜1 − F3 ⎟⎟ + ⎢⎣ ξF3 ⎝ φ0 ⎠⎥⎦ φ0 ⎠ 1 − ξF3/φ0 ⎝

Again, the isomorphism between eq A11 and eq 48 in ref 6 can be

2 ξF3/φ0 βP SPT2a 1 A 2B (ξF3/φ0) = + + ρ 1 − ξF3/φ0 2 (1 − ξF3/φ0)2 3 (1 − ξF3/φ0)3 (A8)

+

φ0

ξF3

⎫ ⎡ ⎛ ξF3/φ0 ξ ⎞⎤⎪ 1 ln⎜⎜1 − F3 ⎟⎟⎥⎬ + 2B⎢ − ⎪ ⎥ ⎢ 1 / − φ ξ φ 2(1 − ξF3/φ0)2 ⎝ ⎣ F3 0 0 ⎠⎦⎭ (ξF3/φ0)2

3(1 − ξF3/φ0)3

+

φ0 ξF3

⎛ ξ ⎞⎤ ln⎜⎜1 − F3 ⎟⎟⎥ φ0 ⎠⎥⎦ ⎝

(A13)

After some straightforward algebraic rearrangement of the term with (A+2B) as coefficient of eq 50 in ref 6, its isomorphism to eq

(A10)

A13 becomes evident. The self-contained expression for

SPT2c. As for SPT2b, we first obtain a more explicit expression for the pressure given by SPT2c by substituting eq A9 into eq 57.

Helmholtz free energy is 5501


Article

The Journal of Physical Chemistry B βF SPT2d βF id = − V V

n

n

∂ [ ∑ ρ βμ − βP] = βμκ ∂ρκ α = 1 α α

∑ ργ ln pIγ (R s = R γ) γ=1

2 ⎧ ⎛ ⎪ ξF3/φ0 ξF3 ⎞ B ⎛ ξF3/φ0 ⎞ ⎟⎟ ⎟⎟ + B + ρ⎨ − ln⎜⎜1 − + ⎜⎜ ⎪ φ0 ⎠ 1 − ξF3/φ0 3 ⎝ 1 − ξF3/φ0 ⎠ ⎩ ⎝

⎡ ⎛ φ ⎞2 ⎤ φ ⎟⎟ ⎥ + ⎢ − 1 + (1 + A) − (A + 2B)⎜⎜ ⎥ ⎢ − − φ φ φ φ ⎝ 0 0⎠ ⎦ ⎣

n

∑ ρα βμα SPT2a − βP SPT2a

⎫ ⎡⎛ φ ⎞ ⎛ ξ ⎞ ⎛ ξ ⎞⎤⎪ φ⎞ ⎛ × ⎢⎜1 − ⎟ ln⎜1 − F3 ⎟ − ⎜1 − 0 ⎟ ln⎜⎜1 − F3 ⎟⎟⎥⎬ ⎢⎣⎝ ξF3 ⎠ ⎝ φ⎠ ⎝ ξF3 ⎠ ⎝ φ0 ⎠⎥⎦⎪ ⎭

α=1 n

=

⎡ ξF3/φ0 − ρ ln(1 − ξF3/φ0) + ρ⎢A ⎢ 1 − ξF3/φ 0 ⎣

SPT2. The expression for pressure given in eq 41 for SPT2 (complete) is already a self-contained one. The self-contained expression for Helmholtz free energy is,

⎛ ξF3/φ ⎞2 ⎤ 0 ⎟⎟ ⎥ − (1 − ξF3/φ0)βP SPT2a + B⎜⎜ ⎝ 1 − ξF3/φ0 ⎠ ⎥⎦

n

∑ ργ ln

= R γ)

γ=1

⎧ ⎛ ξF3/φ0 ⎪ ξ ⎞ Bφ + ρ⎨ − ln⎜⎜1 − F3 ⎟⎟ + ⎪ φ0 ⎠ φ − φ0 1 − ξF3/φ0 ⎩ ⎝

n

∑ ρα βμα SPT2a − βP SPT2a α=1

⎡ ⎛ φ ⎞2 φ ⎢ ⎟⎟ + − 1 + (1 + A) − (A + 2B)⎜⎜ ⎢ φ − φ0 ⎝ φ − φ0 ⎠ ⎣

n

=

⎡ A ξF3/φ0 − ρ ln(1 − ξF3/φ0) + ρ⎢ ⎢ 2 1 − ξF3/φ 0 ⎣ +

(A15)

APPENDIX B

In this appendix, we show in details that the different SPT variants for multicomponent systems obtained in section 2.1 satisfy Gibbs−Duhem equation. For a multicomponent system, Gibbs−Duhem equation is given by

∑ ρα d(βμα )

n α=1

n

∑ ρα α=1

∂(βμα ) ∂ρκ

(B6)

Finally, we have n ⎤ β FSPT2a /V ∂ ⎡⎢ = βμκ SPT2a ∑ ρα βμα SPT2a − βP SPT2a ⎥ = ⎥⎦ ∂ρκ ⎢⎣ α = 1 ∂ρκ

(B1)

when temperature is kept constant. It is to be emphasized that eq B1 is more general than the form given in eq 36 since the latter is used under the condition that the concentrations are kept constant. With eq B1, we can vary the density of any species, e.g., ρκ. In this case, we have ∂(βP) = ∂ρκ

(B5)

∑ ρα βμα SPT2a − βP SPT2a = β FSPT2a /V

n α=1

2⎤ B ⎛ ξF3/φ0 ⎞ ⎥ ⎜⎜ ⎟⎟ 3 ⎝ 1 − ξF3/φ0 ⎠ ⎥⎦

Comparing the RHS of eq B5 to that of eq 44, we see immediately it is nothing else but βFSPT2a/V.

Consistency Check with Gibbs−Duhem Equation

d (β P ) =

∑ ρα [ln(Λα 3ρα ) − 1 − ln(pIα (R s = R α))] α=1

⎛ φ ⎞3⎤⎡⎛ φ ⎞ ξ ⎞ ⎛ φ⎞ ⎛ ⎟⎟ ⎥⎢⎜1 − + 2B⎜⎜ ⎟ ln⎜1 − F3 ⎟ − ⎜1 − 0 ⎟ ⎥ ξF3 ⎠ ⎝ φ⎠ ⎝ ξF3 ⎠ ⎝ φ − φ0 ⎠ ⎦⎢⎣⎝

■

(B4)

Substituting the expression of βPSPT2a given in eq 42 into the RHS of eq B4, we obtain

⎡ ⎛ φ ⎞2 ⎤⎡ φ ⎛ ξ ⎞⎤ Aφ ⎢ ⎟⎟ ⎥⎢1 + 0 ln⎜⎜1 − F3 ⎟⎟⎥ + 1− + 2B⎜⎜ ⎢ φ − φ0 ξF3 ⎝ φ0 ⎠⎥⎦ ⎝ φ − φ0 ⎠ ⎥⎦⎢⎣ ⎣

⎫ ⎛ ξ ⎞⎤ ⎪ ln⎜⎜1 − F3 ⎟⎟⎥⎬ φ0 ⎠⎥⎦⎪ ⎝ ⎭

∑ ρα [ln(Λα 3ρα ) − ln(pIα (R s = R α))] α=1

(A14)

pIγ (R s

(B3)

Consistency Check of SPT2a. From the expression for chemical potential given in eq 46, it can be readily checked that the last three terms on the RHS of eq 46 vanish when the weighted summation over species is taken with ρα as weight. Thus, we have

⎡ ⎡ φ ⎛ ξ ⎞⎤ φ ⎤⎢ ⎥ 1 + 0 ln⎜⎜1 − F3 ⎟⎟⎥ + ⎢1 − 2B − (A + 2B) ⎢⎣ φ − φ0 ⎥⎦⎢⎣ ξF3 ⎝ φ0 ⎠⎥⎦

βF SPT2 βF id = − V V

(κ = 1, 2, ..., n)

(B7)

Thus, it is established that Gibbs−Duhem equation holds for SPT2a. Consistency Check of SPT2b. Proceeding in a similar way as for SPT2a, we start with the expression for chemical potential given in eq 56 and take the weighted summation over species. The last two terms on the RHS of eq 56 vanish after this summation and we have

(κ = 1, 2, ..., n) (B2)

For the convenience of subsequent application, we first rewrite eq B2 in the following equivalent form, 5502


Article

The Journal of Physical Chemistry B n

n

Consistency Check of SPT2d. Taking weighted summation over species with the expression for chemical potential given in eq 62, we obtain,

∑ ρα βμα SPT2b − βP SPT2b = ∑ ρα βμα SPT2a − βP SPT2a α=1

α=1

+ ρ[ln(1 − ξF3/φ0) − ln(1 − ξF3/φ)] − β(P SPT2b − P SPT2a) n

=

∑ ρα βμα

SPT2a

− βP

SPT2a

α=1

⎛ φ ⎞ ⎛ ξ ⎞⎤ + ⎜1 − 0 ⎟ ln⎜⎜1 − F3 ⎟⎟⎥ = ξF3 ⎠ ⎝ φ0 ⎠⎥⎦ ⎝ +

n

⎡ ⎛ ξ ⎞ φ⎞ ⎛ + ρ⎢ − ⎜1 − ⎟ ln⎜1 − F3 ⎟ ⎢⎣ ⎝ ξF3 ⎠ ⎝ φ⎠

∑ ρα βμα SPT2d − βP SPT2d α=1 n

n

∑ ρα βμα

SPT2a

− βP

=

SPT2a

α=1

α=1

β(FSPT2b − FSPT2a ) V

⎧ ⎡ ⎪⎛ φ ⎛ ξ ⎞⎤ φ ⎞⎢ ⎟⎟ 1 + 0 ln⎜⎜1 − F3 ⎟⎟⎥ × ⎨⎜⎜1 + ⎪ φ − φ0 ⎠⎢⎣ ξF3 ⎝ φ0 ⎠⎥⎦ ⎩⎝

(B8)

Equation 54 was used when going to the second equality on the RHS of eq B8, whereas eq 55 was used for obtaining the last equality on the RHS of eq B8. Now, we take the derivative with respect to ρκ on the both sides of eq B8 and obtain

+

n

n ⎤ ∂ ⎡⎢ = ∑ ρα βμα SPT2a − βP SPT2a ⎥ − βμκ SPT2a + βμκ SPT2b ⎥⎦ ∂ρκ ⎢⎣ α = 1

= βμκ SPT2b

=

Equation 61 was used when going from the first equality to the second equality on the RHS of eq B12. Taking the derivative with respect to ρκ on the both sides of eq B12, we obtain n ⎤ ∂ ⎡⎢ ∑ ρα βμα SPT2d − βP SPT2d ⎥ ⎥⎦ ∂ρκ ⎢⎣ α = 1

n

=

SPT2c

α=1 n

=

α=1

⎫ ⎛ φ ⎞ ⎛ ξ ⎞ ⎛ ξ ⎞ ⎤⎪ × ln⎜1 − F3 ⎟ − ⎜1 − 0 ⎟ ln⎜⎜1 − F3 ⎟⎟⎥⎬ φ⎠ ⎝ ξF3 ⎠ ⎝ φ0 ⎠⎥⎦⎪ ⎝ ⎭ =

∑ ρα βμα SPT2b − βP SPT2b + α=1

β (F

−F V

n α=1

α=1

⎧⎡ ⎛ φ ⎞2 ⎤⎡ ⎪⎢ φ ⎛ ξ ⎞⎤ φ ⎟⎟ ⎥⎢1 + 0 ln⎜⎜1 − F3 ⎟⎟⎥ + 2ρB⎨ 1 + + ⎜⎜ φ − φ0 ξF3 ⎝ φ0 ⎠⎥⎦ ⎪⎢⎣ ⎝ φ − φ0 ⎠ ⎥⎦⎢⎣ ⎩

SPT2b

)

2 ⎛ ξF3/φ0 φ ⎞ 1 ⎛ ξF3/φ0 ⎞ ⎟⎟ ⎟⎟ + ⎜⎜1 + − ⎜⎜ 6 ⎝ 1 − ξF3/φ0 ⎠ φ − φ0 ⎠ 2(1 − ξF3/φ0) ⎝

Equation 58 was used when going from the first equality to the second equality on the RHS of eq B10. Taking the derivative with respect to ρκ leads to

⎛ φ ⎞3⎡⎛ φ ⎞ ξ ⎞ ⎛ φ⎞ ⎛ ⎟⎟ ⎢⎜1 − + ⎜⎜ ⎟ ln⎜1 − F3 ⎟ − ⎜1 − 0 ⎟ ξF3 ⎠ ⎝ φ⎠ ⎝ ξF3 ⎠ ⎝ φ − φ0 ⎠ ⎢⎣⎝

n ⎤ ∂ ⎡⎢ ∑ ρα βμα SPT2c − βP SPT2c⎥ ⎥⎦ ∂ρκ ⎢⎣ α = 1

⎫ ⎛ ξ ⎞⎤⎪ × ln⎜⎜1 − F3 ⎟⎟⎥⎬ = ⎪ φ0 ⎠⎥⎦⎭ ⎝

n ⎤ ∂ ⎡⎢ ∑ ρα βμα SPT2b − βP SPT2b⎥ − βμκ SPT2b + βμκ SPT2c ⎥⎦ ∂ρκ ⎢⎣ α = 1

= βμκ SPT2c

n

∑ ρα βμα SPT2 − βP SPT2 = ∑ ρα βμα SPT2d − βP SPT2d

(B10)

=

(B13)

Gibbs−Duhem equations holds for this variant also. Consistency Check of SPT2 (Complete). The summation over species, weighted with their density, of chemical potential given in eq 65 leads to

⎡ ⎧ φ0 ⎛ ⎪ ξ ⎞ φ ⎢⎛ φ⎞ ⎜⎜1 − F3 ⎟⎟ + ×⎨ + − 1 ln 1 ⎟ ⎜ ⎪ ξF3 ⎝ φ0 ⎠ φ − φ0 ⎢⎣⎝ ξF3 ⎠ ⎩

SPT2c

n ⎤ ∂ ⎡⎢ ∑ ρα βμα SPT2c − βP SPT2c⎥ − βμκ SPT2c + βμκ SPT2d ⎥⎦ ∂ρκ ⎢⎣ α = 1

= βμκ SPT2d

∑ ρα βμα SPT2b − βP SPT2b + ρ(1 + A)

n

β(F SPT2d − F SPT2c) V (B12)

(B9)

− βP

∑ ρα βμα SPT2c − βP SPT2c + α=1

Equation B7 was used for going to the second equality of eq B9 which shows that Gibbs−Duhem equation holds. Consistency Check of SPT2c. Taking weighted summation over species with the expression for chemical potential given in eq 59, we obtain,

∑ ρα βμα

⎛ φ ⎞ 2 ⎡⎛ φ⎞ ⎜⎜ ⎟⎟ ⎢⎜1 − + ⎟ 2 ξF3 ⎠ 2(1 − ξF3/φ0) ⎝ φ − φ0 ⎠ ⎢⎣⎝ ξF3/φ0

⎫ ⎛ φ ⎞ ⎛ ξ ⎞ ⎛ ξ ⎞ ⎤⎪ × ln⎜1 − F3 ⎟ − ⎜1 − 0 ⎟ ln⎜⎜1 − F3 ⎟⎟⎥⎬ φ⎠ ⎝ ξF3 ⎠ ⎝ φ0 ⎠⎥⎦⎪ ⎝ ⎭

n ⎤ ∂ ⎡⎢ ∑ ρα βμα SPT2b − βP SPT2b⎥ ⎥⎦ ∂ρκ ⎢⎣ α = 1

SPT2c

∑ ρα βμα SPT2c − βP SPT2c − ρ(A + 2B)

+

(B11)

β(F

SPT2

−F V

SPT2d

n

∑ ρα βμα SPT2d − βP SPT2d α=1

) (B14)

Equation 64 was used when going from the first equality to the second equality on the RHS of eq B14. Taking the derivative with respect to ρκ on the both sides of eq B14, we obtain

Equation B9 was used to go to the second equality on the RHS of eq B11 which is finally nothing else but Gibbs−Duhem equation. 5503


Article

The Journal of Physical Chemistry B n ⎤ ∂ ⎡⎢ ∑ ρα βμα SPT2 − βP SPT2⎥ ⎥⎦ ∂ρκ ⎢⎣ α = 1

=

(13) Sarkisov, L.; Van Tassel, P. R. J. Phys.: Condens. Matter 2008, 20, 333101. (14) Holovko, M.; Patsahan, T.; Dong, W. Pure Appl. Chem. 2012, 85, 115.

n ⎤ ∂ ⎡⎢ ∑ ρα βμα SPT2d − βP SPT2d ⎥ − βμκ SPT2d + βμκ SPT2 ⎥⎦ ∂ρκ ⎢⎣ α = 1

= βμκ SPT2

(B15)

Equation B15 shows that Gibbs−Duhem equations holds also for SPT2 (the complete version). This closes the present appendix which establishes the thermodynamics consistency of all the SPT variants for a multicomponent fluid in a multicomponent matrix.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b02957. Tables containing numerical results for chemical potentials of the two fluid species from GCEMC and the different variants of SPT. (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: 0033-472728844. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Dr. J. X. Hou for some discussions on this work at its initial stage and C.Z. Qiao for check reading of the manuscript. The support from CNRS-NASU cooperation project (N° 23999) is acknowledged. W. Chen is grateful for a French government scholarship (thèse en cotutelle). S.L. Zhao is grateful for an invited professor position at Ecole Normale Supérieure de Lyon. Parts of the present work were accomplished during the visits of W. Dong at East China University of Science and Technology, supported by the 111 project of the Ministry of Education of China (N° B08021), and at Institute of Theoretical Physics of Chinese Academy of Sciences, supported by an invited research professor position.

■

REFERENCES

(1) Reiss, H.; Frisch, H. L.; Lebowitz, J. L. J. Chem. Phys. 1959, 31, 369. (2) Reiss, H.; Frisch, H. L.; Helfand, E.; Lebowitz, J. L. J. Chem. Phys. 1960, 32, 119. (3) Lebowitz, J. L.; Helfand, E.; Praestgaard, E. J. Chem. Phys. 1965, 43, 774. (4) Holovko, M.; Dong, W. J. Phys. Chem. B 2009, 113 (6360), 16091. (5) Chen, W.; Dong, W.; Holovko, M.; Chen, X. S. J. Phys. Chem. B 2010, 114, 1225. (6) Patsahan, T.; Holovko, M.; Dong, W. J. Chem. Phys. 2011, 134, 074503; 2016, 144, 099903. (7) Madden, W. G.; Glandt, E. D. J. Stat. Phys. 1988, 51, 537. (8) Vega, C.; Kaminsky, R. D.; Monson, P. A. J. Chem. Phys. 1993, 99, 3003. (9) Meroni, A.; Levesque, D.; Weis, J. J. J. Chem. Phys. 1996, 105, 1101. (10) Rosinberg, M. L. In New Approaches to Problems in Liquid State Theory; NATO Science Series C, Vol. 529; Caccamo, C., Hansen, J.-P., Stell, G. , Eds.; Kluwer: New York, 1999; page 245. (11) Gelb, L. D.; Gubbins, K. E.; Radhakrishnan, R.; SliwinskaBartkowiak, M. Rep. Prog. Phys. 1999, 62, 1573. (12) Schmidt, M. J. Phys.: Condens. Matter 2005, 17, S3481. 5504


Scaled Particle Theory for Multicomponent Hard Sphere Fluids

Recommend Documents