Phase Equilibrium with External Fields: Application to Confined Fluids

Article pubs.acs.org/jced

Phase Equilibrium with External Fields: Application to Confined Fluids Noura Dawass, Michelle L. D’Lima, Ioannis G. Economou, and Marcelo Castier* Chemical Engineering Program, Texas A&M University at Qatar, P.O. Box 23874, Doha, Qatar ABSTRACT: This work addresses the problem of finding the equilibrium spatial segregation of components in systems that have multiple regions, each of them subject to the effect of external fields. The specifications are the temperature, region volumes, component amounts, and additional variables that characterize the effect of such fields. The formulation leads to the mathematical problem of minimizing the Helmholtz energy of the system, subject to constraints that represent component mass balances and volume conservation equations applied to each region. Among other uses, the approach is suitable for determining the equilibrium conditions in batch adsorption. The formulation is general but the focus of this work is on the compositional segregation in isothermal reservoirs due to gravity and the spatial segregation of components close to pore walls. Calculations using the Steele and DRA potentials to model fluid−wall interactions demonstrate the formulation, and its solution procedure provides results that are generally in very good agreement with experimental data and simulations reported in the literature. The formulation enables the prediction of meaningful trends for local composition profiles for fluids inside pores, at a coarser level than those from molecular simulations, but with much smaller computational load.

■

INTRODUCTION The conventional formulation of equilibrium calculations for chemical process design often neglects the effect of external fields. As a consequence, each phase present is homogeneous, in the sense that all intensive properties have equal value anywhere within it. However, there are situations of practical interest in which the action of one or more external fields induces phase inhomogeneity. For example, the action of the gravitational field gives origin to compositional grading in oil reservoirs, that is, the oil composition depends on depth. Gravitational segregation studies follow two paths, namely, thermodynamic equilibrium assuming isothermal conditions1−3 and nonequilibrium calculations that take thermal diffusion into account.4−6 The presence of a gravitational field also induces the spatial segregation of light particles in sedimentation equilibrium.7 The separation of chemical components in ultracentrifuges8−13 exemplifies inhomogeneities induced by a centrifugal field. Another example is the action of an electrical field, which can create compositional grading and phase separation in electrolyte solutions.14 The effect of magnetic fields on wax precipitation15 and on the critical point of hydrocarbons mixtures16 has been studied recently. Adsorption is another phenomenon that falls within the scope of systems subject to the effect of external fields, because the walls of the solid adsorbent impose a field that causes spatial segregation of the fluid close to the solid’s wall. Therefore, the general problem of accounting for external fields on the conditions of thermodynamic equilibrium may span multiple size scales, from hundreds of meters in depth, as in compositional grading in oil reservoirs caused by gravity, to nanometric distances from solid walls, as in fluid adsorption. © XXXX American Chemical Society

This work presents a formulation to account for the action of one or more external fields on the conditions of thermodynamic equilibrium, with focus on fluid segregation caused by adsorption. Direct experimental probing of local structure and density of adsorbed fluids is difficult.17 Formal theoretical treatments to study this sort of segregation include Monte Carlo simulations, molecular dynamics, and classical density functional theory. Monte Carlo and molecular dynamics simulations18,19 provide much insight into the effect of the walls on the structure of a fluid confined within a pore, but generally involve long calculations, still too demanding for routine applications to chemical process design. Classical density functional theory20,21 demands less computational effort but models derived from it are seldom used for chemical process design, possibly due to their level of complexity. An interesting aspect is that studies of adsorption commonly adopt grand-canonical ensemble formulations, that is, the chemical potentials of the components in the bulk phase are specified and the adsorbed amounts are calculated. A technique of intermediate complexity and computational load is the development of equations of state (EOS) for confined fluids. For instance, under a suitable set of simplifying assumptions, it is possible to incorporate the effect of the confining wall into the statistical mechanical derivation of expressions for the Helmholtz energy, from which EOS can be Special Issue: In Honor of Kenneth R. Hall Received: March 8, 2016 Accepted: April 26, 2016

A

DOI: 10.1021/acs.jced.6b00209 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

Figure 1. Various representations of a system comprising a bulk region (white background) and confined regions: (a), bulk region and multiple replicas of a confined region (stripes); (b), bulk region and a confined region (background with color), with a single replica depicted. The arrows represent molecule−wall interactions; (c), bulk region and two confined regions (backgrounds with color), with a single replica depicted. The arrows represent molecule−wall interactions and the gravitational field; (d), schematic of the grid used to model a system with one bulk region and a confined region (background with color), with a single replica depicted. The arrows represent molecule−wall interactions and the gravitational field.

is important to account for when dealing with reservoir fluids trapped in nanopores.26 Using this formulation, it is possible to predict the composition of the reservoir fluid as a function of depth within the reservoir and distance to the confining walls, inside the rock pores. Unlike many papers on adsorption equilibrium, which adopt a grand-canonical ensemble formulation, the approach to adsorption and other problems with external fields adopted here uses a canonical ensemble formulation. In it, determining equilibrium conditions requires the minimization of the system’s Helmholtz energy. The specifications are the volume of each region (e.g., the volume of each pore type that is available in a heterogeneous adsorbent), and the system’s temperature and component amounts. These specifications are relevant, for example, to batch adsorption problems. The system’s Helmholtz energy is the summation of the Helmholtz energies of each of the regions. However, the thermodynamic properties change with position within each inhomogeneous region, and so does the molar Helmholtz energy. Therefore, it is necessary to integrate over the space of each region to find the corresponding contribution. For EOS typically used for chemical process design, the corresponding integrals do not have analytical expressions, so their numerical integration is necessary. This paper describes this formulation and the solution procedure in detail, and how it interfaces with known EOS and expressions that represent the effect of external fields on the fluid. It also shows how the model predicts the spatial segregation of the chemical components, as in results of molecular simulations, albeit at a coarser level. The examples deal with the adsorption properties of confined pure components and mixtures. One of the examples shows a prediction of the simultaneous effect of confinement and gravity, which is a situation of interest to the oil and gas industry.

obtained. Examples include the use of the generalized van der Waals theory as a starting point to develop models that extend the van der Waals,22 Redlich−Kwong, Soave−Redlich−Kwong, and Peng−Robinson23 EOS to correlate and predict the properties of fluids confined within cylindrical and spherical pores. The perturbed-chain statistical associating fluid theory (PC-SAFT EOS) has been recently coupled with the Young− Laplace equation to model phase equilibrium within cylindrical pores.24 A convenient feature of these extended models is that they revert to the originating EOS for bulk fluids if the confinement space is large and has nonattractive pore walls. Calculations with these models are in the same vein as those with their originating EOS, and they can be used to predict capillary condensation. Furthermore, it is possible to simulate heterogeneous solids by specifying the presence of pores with more than one set of characterizing parameters, such as pore radius and fluid−wall interactions. One shortcoming of these models is that they cannot be applied to evaluate the local fluid distribution inside the pore. An approach to adsorption that overcomes this shortcoming is to add terms to the internal energy of the fluid in order to account for the fields imposed by the confining solid walls, as in the multipotential theory of adsorption (MPTA).25 Further, in adsorbents with multiple pore types, the adsorbed fluid is submitted to a different field depending on the pore in which it is; thus, this approach allows the detailed modeling of adsorption on heterogeneous adsorbents. In addition, the modeling can also include the simultaneous effect of additional external fields, such as gravitational, electrical, and magnetic fields. Such a general formulation allows the modeling and simulation of a variety of phenomena in which multiple fields affect the conditions of thermodynamic equilibrium. An example is the modeling of compositional grading in oil reservoirs, accounting for the simultaneous effect of confinement and gravity. Confinement affects system properties, which B



■

Article

FORMULATION The formulation of this model adopts the following assumptions: 1. The temperature, T, is fixed and known. 2. The amount of each component i, ni, is fixed and known, and there is no chemical reaction. 3. The system has one or more regions, depending on the number of external fields present; each region m has a fixed and known total volume, Vm. 4. The total volume of a region may comprise a known number of replicas of smaller volume; for example, the total pore space may comprise a large number of identical pores, each of them of small volume. 5. The integral of the Helmholtz energy, A, over the space of each region is approximated by a summation of the Helmholtz energy of discrete grid elements, thereby accounting for inhomogeneity. 6. Each region may be subject to the effect of one or more external fields. For example, it may be necessary to account for the effect of gravitational field and interaction of the fluid with a confining wall simultaneously. 7. The effect of confining walls, when they exist in a problem, is fully accounted for by their interaction potentials with the molecules of the confined fluid. To account for the inhomogeneity of the system, the system is initially divided into regions. As an example, consider Figure 1a, which is a schematic of a system with one bulk region (white background) and a porous region of several identical slit pores (stripes), hereafter referred to as replicas. Therefore, the total pore volume and the total amount of each component in that region is the outcome of replicating the result of each individual pore. Figure 1b represents a system with just one pore that is magnified, and its contribution needs to be replicated. In Figure 1b, the arrows represent molecule−wall interactions. Figure 1c depicts a more complex situation, in which there is one bulk region and two different confined regions. For simplicity, only one magnified pore of each of these confined regions is shown. Molecule−wall interactions occur only in the confined regions, while both the bulk and confined regions are subject to the gravitational field. The space of each region may be split further into layers that define grid elements to account for inhomogeneity. Figure 1d shows the grids used to model a system with one bulk region and a confined region, where a single replica is depicted. In addition to the data that characterize the components present and the effect of the external fields, the temperature, component amounts, and region volumes are specified. Under this set of specifications, the minimization of the Helmholtz energy gives the equilibrium condition. Taking into account the multiple regions, replicas, and layers, the Helmholtz energy is given by r̂

A=

Equation 1 then takes the form: r̂

A=

The Helmholtz energy can be split in an internal contribution, associated with the EOS, and a contribution due to the effect of external fields; that is, A*jm = A*jm,in + A*jm,f

(4)

where the superscripts “in” and “f” denote the internal and field contributions, respectively. The total field contribution may result from the simultaneous effect of multiple individual fields, which give rise to layers in different directions within a given region. Thus, eq 4 can be rewritten as r̂

A=

∑

lm

∑ (A*jm,in +

m=1 j=1

fm

∑ A*fjm,f ) (5)

f =1

where f m represents the number of external fields acting on region m. As discussed later on, the procedure for minimizing the Helmholtz energy requires first and second derivatives Ajm * with respect to mole numbers. The first derivative is given by ⎛ ∂A* ⎞ * = ⎜ jm ⎟ μijm ⎜ ∂n * ⎟ ⎝ ijm ⎠T , V * , n * jm

≠ i ,jm

⎛ ∂Ajm ⎞ ⎟⎟ = ⎜⎜ ⎝ ∂nijm ⎠T , V

= μijm

jm , n≠ i ,jm

(6)

where μijm * is the chemical potential of component i, in layer j of region m, considering all replicas and the effect of all fields. The symbol n*ijm denotes the number of moles of component i at the same location. The notation ≠i indicates components other than i. To derive the two rightmost terms of eq 6, the relationship between n*ijm and nijm, which is the corresponding amount in a single replica, is used. It is * = rmnijm nijm

(7)

Using eq 7, the second derivative of A*jm with respect to mole numbers is given by ⎛ ⎞ ⎜ ∂ ⎛ ∂A*jm ⎞ ⎟ ⎜ * ⎜⎜ * ⎟⎟ ⎟ ⎜ ∂nkjm ⎝ ∂nijm ⎠ * * ⎟ T , V jm , n≠ i ,jm ⎠ ⎝ T ,V * ,n * jm

=

1 ⎛⎜ ∂μijm ⎞⎟ rm ⎜⎝ ∂nkjm ⎟⎠

T , Vjm , n≠ k ,jm

⎛ ∂μ * ⎞ ijm ⎟ = ⎜⎜ ⎟ * ∂ n ⎝ kjm ⎠T , V * , n *

≠ k ,jm

jm

≠ k ,jm

(8)

The first and second derivatives of Ajm * are the result of adding the derivatives of the internal and external contributions in that element. The following subsections discuss the evaluation of the internal and external field contributions. Internal Contribution. The reference state for the calculation of Helmholtz energies is chosen to be the chemical components at the system temperature and at a reference pressure, P0, in the hypothetical ideal gas state. To alleviate the notation, the indexes that denote the region and layer are temporarily dropped from the equations that follow. The path to deriving a working expression for the molar Helmholtz energy (a) as a function of temperature, molar volume (v), and component mole fractions (x), starts from an explicit expression for the molar residual Helmholtz energy,

lm

(1)

where r̂ is the number of regions, rm is the number of replicas of region m, lm is the number of layers in each replica of region m, and Ajm is the Helmholtz energy of layer j in each replica in region m. The definition of A*jm which is the joint contribution of all replicas of layer j in region m follows: A*jm = rmAjm

(3)

m=1 j=1

∑ ∑ rmAjm m=1 j=1

lm

∑ ∑ A*jm

(2) C



Article

ares(T,v,x), as is the case of most modern EOS based on statistical thermodynamics. The molar residual Helmholtz energy is defined here as the difference between the molar Helmholtz energies of the fluid and of the ideal gas at the same temperature, molar volume, and component mole fractions; that is, a res(T , v , x̲ ) = a(T , v , x̲ ) − a ig (T , v , x̲ )

a res(T , vjm , x̲ jm) RT

(9)

ĉ

am =

(16)

ĉ

∑ ∑ xixjaij (17)

i=1 j=1

aij =

i=1 ĉ

aiaj (1 − kij)

(18)

ĉ

bm =

= RT ∑ xi ln xi

v

2

Pv 1 dv = −RT ln 0 v RT

ĉ

cm =

(11)

∑ xici

(20)

i=1

The change in molar Helmholtz energy between states 3 and 4 is the residual molar Helmholtz energy, whose expression depends on the model adopted for a given calculation. Reintroducing the indexes for region and layer gives the expression for the Helmholtz energy:

Parameter ci is calculated as follows: ci = sibi

(21)

where si represents a parameter fitted to experimental liquid density data and its value is generally negative. Field Contribution. Different types of external fields may affect a given region, some of which are exemplified here. Gravitational Field. The effect of the gravitational field is relevant to compositional grading, in which the pressure and composition of a mixture change with the vertical position coordinate, like depth in an oil reservoir. Accounting for the effects of all replicas, the corresponding contribution of this field to the Helmholtz energy function is

⎡ ĉ ,in * * ln xijm Ajm = RT ⎢∑ nijm ⎢⎣ i = 1 ⎛ P0vjm a res(T , vjm , x̲ jm) ⎞⎤ * ⎟⎟⎥ − (∑ nijm)⎜⎜ln − RT RT ⎝ ⎠⎥⎦ i=1

(19)

where kij is a binary interaction parameter. The expressions for ai and bi as a function of critical temperature, critical pressure, and acentric factor are readily available in the literature.27 The volume translation parameter cm is given by

where ĉ denotes the number of components and g is the molar Gibbs energy. The change between states 2 and 3 is given by

∫v

∑ xibi i=1

(10)

i=1

2

(15)

Parameters am and bm are evaluated using the classical mixing rules:

= RT ∑ xi ln xi − R ΔT

⎛ ∂a ⎞ig ⎜ ⎟ dv = −RT ⎝ ∂v ⎠T , x̲

am (v + cm)(v + cm + bm) + bm(v + cm − bm)

−

ĉ

v

∫∞

⎛ P in ⎞ ⎜ jm − 1 ⎟ dvjm ⎜ RT vjm ⎟⎠ ⎝

RT v + cm − bm

P=

Δa12 = Δg12ig − Δ(Pv ig)12

∫v

vjm

The first term of the integrand corresponds to the contribution of the real fluid and the second term is the ideal gas contribution, both lie between infinite molar volume and the molar volume of the fluid in layer j of region m. Peng−Robinson EOS with Volume Translation. Volume translation improves density predictions of cubic EOS, and this refinement is essential for accurate evaluations of the amounts adsorbed within porous materials. The expression for the Peng−Robinson EOS27 with volume translation28 is as follows:

A four-state path is adopted, as follows: State 1: each pure component in the hypothetical ideal gas state at T and P0; State 2: ideal gas mixture with mole fractions x at T and P0, at which v2 = v(T,P0) = RT/P0; State 3: ideal gas mixture at T, with mole fractions x and molar volume v, equal to the molar volume of the system; State 4: real fluid at T, with mole fractions x, and molar volume v. The effect on the molar Helmholtz energy of forming an ideal gas mixture isothermally (State 2) from the pure components (State 1) can be computed as follows:

Δa 23 =

=−

ĉ

(12)

where xjm denotes the vector of mole fractions of all components in layer j of region m. For mole fractions and molar volumes, it holds that

ĉ

* Mi A*jm,f = gĥ jm ∑ nijm i=1

* = xijm xijm

(13)

v*jm = vjm

(14)

(22)

where ĝ is the gravitational constant, Mi is the molar mass of component i, and hjm is the height of layer j in region m. Field Imposed by a Flat Confining Wall: Steele Potential. There are several expressions available to describe the effect of a confining wall, depending on the pore geometry and type of fluid−wall interactions. The Steele potential29 is a common choice for modeling the confining effect of the slit pores of activated carbon. In a slit pore, there are two parallel confining walls, assumed to be graphite-like surfaces,30,31 that contribute

For this reason, the superscripted asterisk (∗) does not appear in these variables in eq 12. An alternative to determining the Helmholtz energy is to start from an EOS and substitute the following relation in eq 12, D



■

to the potential. Therefore, the overall field contribution of the Helmholtz energy is computed as follows:

SOLUTION PROCEDURE This section presents the discretization approach taken to consider inhomogeneity within each region. It also discusses the selection of independent variables, the generation of initial estimates, and a procedure to handle independent variables whose values may cause numerical difficulties. Grid Setup. The thermodynamic equilibrium formulation requires minimizing the Helmholtz energy, which is a result of adding the contributions of the regions present. Calculating the Helmholtz energy of each inhomogeneous region requires integration over the region. Because of the mathematical complexity of most thermodynamic models, it is necessary to replace such integration with a numerical approximation. This is accomplished by establishing layers within each region that generate a three-dimensional grid for each inhomogeneous region. There are several ways to implement this. In the cases considered in this paper, which are for slit pores, the direction of the fields present and the axes adopted for the problem are Cartesian and aligned. Discretization occurs in the direction perpendicular to each field, as illustrated in Figure 1. For simplicity, grid elements of equal volume are used within each region. However, grids in different regions may have different sizes, based on problem specifications. The center of each grid element is taken as its position for field intensity calculations. Alternatives to this approach, not adopted in this work, are the use of Gaussian quadrature or successive grid refinements using Romberg integration, among others. Numerical Procedure. The Helmholtz energy is minimized using Murray’s34 method as implemented by Michelsen.35 Since the systems studied in this work do not undergo phase change, a phase stability test is not utilized. For the calculations, distribution factors are utilized, as discussed in the text that follows. Independent Variables. It is possible to choose several sets of independent variables to minimize the Helmholtz function, with the general goal of establishing the amount of each component in each layer of each region. A convenient choice is the distribution factor, θ*ijm, defined as follows:

ĉ

A*jm,f =

fw fw * [Aijm (zjm) + Aijm (Hm − zjm)] ∑ nijm i=1

(23)

fw where Afw ijm(zjm) and Aijm(Hm − zjm) represent the contribution by the near and distant walls, respectively. The symbol Hm represents the wall to wall distance, while the symbol zjm represents the distance from the center of layer j to a confining wall, all within region m. The Steele potential for each component, i, can be defined for any distance z (equal to zjm or (Hm − zjm)) as

⎡ 2 ⎛ σ ⎞10 ⎛ σ ⎞4 s , im s , im fw ⎟ ⎟ Aijm (z) = 2πρsm Δmεs , imσs2, im⎢ ⎜ −⎜ ⎝ z ⎠ ⎢⎣ 5 ⎝ z ⎠ −

⎤ ⎥ 3Δm(0.61Δm + z)3 ⎥⎦ σs4, im

(24)

where εs,im is the solid−fluid energy interaction parameter and σs,im is the solid−fluid molecular diameter, both in region m, which are calculated using the Lorentz−Berthelot combining rules. The symbols ρsm and Δm are the density and interlayer spacing of graphite, respectively, both in region m. The Dubinin−Radushkevich−Astakhov Potential. Another potential used to describe the effect of a confining wall is the Dubinin−Radushkevich−Astakhov (DRA) potential.32,33 Its original formulation predicts the adsorption of pure components, disregarding fluid compressibility. Shapiro and Stenby25 modified the potential to account for adsorbate compressibility, and also adapted it to perform multicomponent adsorption calculations. Therefore, the potential experienced by each component, εi, is given by ⎡ ⎛ ⎡ v ̂ ⎤⎞1/ β ⎤ εi = ⎢ −ε0i⎜ln⎢ 0 ⎥⎟ ⎥ ⎢⎣ ⎝ ⎣ v ̂ ⎦⎠ ⎥⎦

(25)

* = θijm

where ε0i is the characteristic energy of adsorption of component i. Parameter β is related to the solid heterogeneity and is equal to 2 for activated carbon.25 The total pore volume, v̂0, is fitted to pure adsorption experimental data along with ε0i parameters. In the MPTA implementation of the DRA potential, calculations are performed at each volume from v̂ = 0 to v̂ = v̂0. In this work, the volume ratio that appears in eq 25 is replaced by a distance ratio. The DRA contribution to the Helmholtz energy based on the distance ratio is given by

* nijm ni

(27)

Due to mass conservation, note that for each component i r̂

lm

* =1 ∑ ∑ θijm (28)

m=1 j=1

In this work, the largest θ* value of each component i is chosen to be the dependent variable. The symbol θiJM * represents its value and the capital subscripts J and M indicate the layer and region where it occurs. Therefore, the value of θ*iJM is obtained as follows:

1/ β ⎡ ⎛ ⎡ ĉ ⎤⎞ ⎤ ⎢ ⎜ ⎢ Hm ⎥⎟ ⎥ ,f * * Ajm (zjm) = −∑ nijm⎢ε0i⎜ln ⎟ ⎥ i=1 ⎢⎣ ⎝ ⎢⎣ zjm ⎥⎦⎠ ⎥⎦

⎛ ⎡ H ⎤⎞1/ β c ̂ * ε0i = −⎜⎜ln⎢ m ⎥⎟⎟ ∑ nijm ⎝ ⎢⎣ zjm ⎥⎦⎠ i = 1

Article

r̂

* = 1 − (∑ θiJM

lm

* )≠ {iJM} ∑ θijm

m=1 j=1

(29)

The subscript ≠{iJM} on the right-hand side of eq 29 indicates that all θ*-values of component i are added, except the largest one. At the beginning of each iteration, the largest θ* value of each component i is flagged as a dependent variable and the numerical procedure generates new values for all the independent θ* values. Then, the dependent value is computed

(26)

where Hm is the distance beyond which the adsorption potential effects are negligible. Calculations are performed at distance ratios zjm/Hm sampled in the open interval (0,1). E



Article

elements are equal to zero, and the equilibrium molar distribution is found. Ideal Gas Initial Estimates. To generate initial estimates for the distribution factors, ideal gas behavior is assumed. This simplifies the problem, and allows for the number of moles in each grid element of the system to be found analytically. At equilibrium, the chemical potential of each component i should be the same in all grids, that is,

using eq 29, allowing the evaluation of the objective function and its derivatives during the subsequent iteration. Also, it is possible to write the Helmholtz energy derivatives in terms of the independent distribution factors. The first derivative of the Helmholtz function is ⎛ ∂A ⎞ * − μ* ) ⎜⎜ ⎟ = ni(μijm iJM * ⎟⎠ ⎝ ∂θijm T , V , θ′*

(30)

μi*,ref

These derivatives are evaluated at constant volume of all the layers in all the regions, denoted by V , and at constant value of all independent θ* values, other than θijm * . The latter variable set is denoted by θ′* in eq 30. Omitting the variables kept constant during differentiation, the second derivative of the Helmholtz function is ⎛ ⎞ ∂ 2A ⎜ ⎟ ⎜ ∂θ * ∂θ * ⎟ ⎝ i , ji , mi l , jl , ml ⎠ ⎡⎡ ⎛ ∂μ * ⎞⎤ ⎤ ⎢⎢ ⎜ i , ji , mi ⎟⎥ ⎥ ⎢ ⎢(δmi , mlδ ji , jl − δmi , Mlδ ji , Jl)⎜ * ⎟⎥− ⎥ ⎢⎣ ⎝ ∂nl , ji , mi ⎠⎦ ⎥ ⎢ ⎥ = ninl ⎢⎡ ⎛ ∂μ * ⎞⎤ ⎥ i , Ji , Mi ⎟⎥ ⎥ ⎢ ⎢ (δ δ J , j − δ Mi , Mlδ J , J )⎜ M , m i l ⎢ ⎢ i l ⎜ i l * ⎟⎥ ⎥ ⎝ ∂nl , Ji , Mi ⎠⎦ ⎥⎦ ⎢⎣ ⎣

RT

=

* μijm (32)

RT

The symbol μi,ref * represents the chemical potential of component i in a certain grid element that is taken as a reference. The symbol μijm * denotes the chemical potential of component i in any other grid element, characterized by the indexes j and m for layer and region, respectively. Omitting the variables kept constant during differentiation, it follows that ln xi ,ref + ln Pref +

*,f 1 ∂A ref RT ∂ni*,ref

= ln xijm + ln Pjm +

,f 1 ∂A*jm * RT ∂nijm

(33)

Applying the ideal gas EOS to eliminate the pressures from eq 33, one obtains: (31)

* nijm

where δ represents the Kronecker delta function. The lowercase subscripts refer to the independent grid elements of component i and l, while the capital subscripts represent the dependent ones (those having the highest distribution factor) of component i and l. Finally, the algorithm in Figure 2 summarizes the steps followed during the numerical optimization for a given iteration. In this work, to begin the minimization of the Helmholtz energy, a set of initial estimates are calculated as discussed in the next section. Then the steps demonstrated in Figure 2 are implemented and repeated until the gradient

Vjm

=

ni*,ref Vref

⎛ ⎛ *,f ∂A*jm,f ⎞⎞ 1 ⎜ ∂A ref ⎟⎟ exp⎜⎜ − ⎜ ∂n * ⎟⎟ * RT n ∂ ijm ⎠⎠ ⎝ i ,ref ⎝

(34)

The grid element with the lowest field potential is taken as a reference for other grid elements. The equations of the form given in eq 34 and the mass balance equation compose a linear system of equations that is solved to obtain the number of moles in each grid element. This procedure is repeated for each component. Finally, the amounts of each component in each grid element are converted to distribution factors in order to serve as a good starting point for the optimization process. Inactive Variables. During the course of the iterative procedure, some of the independent variables may become very small. Because of size exclusion, the potential energy becomes very large and therefore, the probability of observing a molecule at this position is very low and equal to zero for practical purposes. During the iterative procedure, if a distribution coefficient (θ*) falls below the value of 1 × 10−20, it is set to zero and removed from the list of independent variables. Physical violations may also occur when the density in a specific grid element is very large. In other words, this happens when the molar volume, vjm, is lower than the corresponding value of bm calculated from the EOS, which is the packing volume for that grid element. To overcome this problem, the excess amounts of all the violating grids are calculated, and are then redistributed in the nearby grid elements that do not experience any violations.

■

RESULTS AND DISCUSSION This section presents results for adsorption on activated carbon, through modeling as slit pores, using the Steele and DRA potentials. For the components studied, the parameters that characterize adsorption are taken from literature. This allows for a direct comparison between this work and similar approaches. The model results are compared to experimental

Figure 2. Algorithm detailing the evaluation of Helmholtz energy and its derivatives. F



Article

the adsorbed fluid accumulated near the pore walls due to the maximum potential effect experienced there. Moving away from the wall and toward the center of the pore, the solid−fluid interactions become weaker. Less amounts are adsorbed, until a near-constant density is reached in the middle of the pore. Another important observation is related to the effect of the bulk phase pressure on the density profile. As the pressure increases, the amounts inside the pore are expected to increase, with higher pressures corresponding to higher density maxima. Finally, varying the wall-to-wall distance of the slit pore affects the width of the density peaks, as displayed in Figure 3. Once the density of each layer in the pore is calculated, the density of each component in the confined region m is obtained, which is defined as

data for pure components, binary, and ternary mixtures. These data are reported at a given temperature and bulk phase pressure and, for mixtures, also the bulk phase composition. The methodology proposed in this work requires specification of the temperature, bulk, and pore volumes, and total amount of each component (T,V,N). Minimization of the Helmholtz energy under these specifications yields the density and composition in each region and the bulk phase pressure (P). To allow comparisons to the experimental data and validate the procedure’s implementation, the numbers of moles in the system were varied until matching the bulk phase pressure of each experimental data point. In addition to adsorption cases, this section presents a conceptual example in which the fluid is confined in slit-like pores and subject to the gravitational field, as in a deep porous reservoir. Results with the Steele Potential. This example investigates the adsorption of methane on activated carbon, which has been presented in work by Li et al.31 that uses DFT with the volume-translated Peng−Robinson EOS. Activated carbon is modeled as a slit pore, with pore walls that interact with the fluid molecules according to the Steele potential. The same EOS and fluid-wall interaction model are used here. Adsorption on a Single Slit Pore. Figure 3 displays the local density profile of methane at T = 298 K and P = 0.1, 1, and 2

ρim =

1 Hin, m

∫0

Hm

ρim (zm) dzm

(35)

where Hin,m is the internal pore width (excluding the radius of the adsorbent molecule). Using numerical discretization, the adsorbed amount is calculated by adding the number of moles at each distance and dividing by the pore volume in region m(Vm), according to l

m * Hm ∑ j = 1 nijm ρim = Hin, m Vm

(36)

where the factor Hm/Hin,m is used to obtain the density within the effective pore volume. Figure 4 shows the adsorption of methane for wall-to-wall distances of 2, 3, and 6 nm, at pressures up to 30 MPa. The

Figure 4. Methane adsorption isotherm in activated carbon at 303.15 K, and H = 2, 3, and 6 nm. Solid lines represent predictions of this work, while the points are predictions by DFT.31

results are plotted along with the DFT results of Li et al. There is excellent overall agreement between the two results, especially in the range of 5−30 MPa, where the maximum deviation is 4%. Additionally, the plots agree with the DFT results of Li et al. and demonstrate the effect of pore size on the amounts adsorbed. As the pore width decreases, the density inside the pore at a given bulk phase pressure increases. The model is capable of predicting compositional profiles inside pores. Figure 5 shows results at 323.15 K and bulk phase pressure of 0.5 MPa for an equimolar bulk phase mixture of methane and propane, confined in slit pores whose width is 10 nm. The Steele potential parameters are shown in Table 1. The

Figure 3. Local density profiles for methane at 298 K in slit pores of different widths and bulk phase pressures of 0.1, 1, and 2 MPa. Dashed lines represent bulk densities at the respective pressures. (a) H = 2 nm, (b) H = 3 nm, and (c) H = 6 nm.

MPa. The three profiles are plotted for activated carbon at different pore widths, namely H = 2, 3, and 6 nm. The plots highlight a few features of adsorption in slit pores at the given pressure range. Since both walls are characterized using the same energy, the density distribution is symmetric, with most of G



Article

Table 2. Pore Size Distribution and Pore Volumes Hm

PSD (Hm)

Vm

nm

m3/(kg/nm)

m3/kg

0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

3.030 7.333 4.029 8.653 9.925 7.373 4.020 1.751 6.462 2.109 6.280

× × × × × × × × × × ×

10−6 10−5 10−4 10−4 10−4 10−4 10−4 10−4 10−5 10−5 10−6

3.818 2.381 6.341 9.289 8.649 5.697 2.886 1.199 4.286 1.369 4.010

× × × × × × × × × × ×

10−6 10−5 10−5 10−5 10−5 10−5 10−5 10−5 10−6 10−6 10−7

rconf ̂

Nads =

Figure 5. Compositional profile for an equimolar mixture of methane and propane in an activated carbon pore of H = 10 nm, at 323.15 K and bulk phase pressure of 0.5 MPa. DFT predictions are from Li et al.31

∑ ρm Vm(Hm) m=1

(38) 36

Figure 6 displays the experimental data of Qiao et al., the DFT results of Li et al.,37 and two curves that present the

Table 1. Steele Potential Parameters Obtained from the DFT Model (Li et al.31) solid parameters −3

ρs (nm ) σs (nm)

ϵs/kB(K) s

114 δ (nm) 0.340 ϵs/kB (K) adsorbed fluid parameters

0.335 20

CO2

CH4

C3H8

1760 −0.0623

1178 −0.1533

1866 −0.0869

composition profiles are symmetric and only the wall-to-center distance is shown. Next to the pore wall, the methane mole fraction peaks, followed by a peak in the propane mole fraction. As the distance from the wall increases, the adsorption effect becomes less pronounced and both mole fractions tend to 0.5, which is their value in the bulk phase. Figure 5 also includes the propane composition profile for the same mixture as predicted by the DFT model, demonstrating that the two approaches agree qualitatively in this example. Adsorption on Multiple Slit Pores. The work by Li et al.31 describes the use of pore size distribution (PSD) to characterize the structure of the pore (heterogeneous confinement) and predict adsorption of mixtures based on pure components. The model can be validated by comparison with experimental results by Qiao et al.36 The approach taken in this section is to use their PSD characterization to obtain the adsorption isotherm of methane within a pressure range of 0.01−0.11 MPa at 363 K. The resulting calculations are then compared to the experimental data and the DFT results. To simplify the calculations, pores with small volumes are neglected. Table 2 presents the PSD values for the pores considered. This PSD characterization is used to obtain the pore volume, Vm(Hm), of region m with a specific pore width of Hm using the following relation,

Figure 6. Adsorption isotherm of CH4 on activated carbon at 363.15 K, εCH4/kB = 865 and 1178 K.31,36

results of this work. One of these curves was evaluated with the ε/kB value of methane equal to 1178 K, fitted by Li et al. using the attractive part of the Lennard-Jones potential at a low temperature. Setting the ε/kB value of methane to 865 K, it is possible to obtain results that are in much better agreement with the experimental data than those of the DFT calculations of Li et al.37 Comparison with Molecular Simulations. To understand the capabilities and limitations of this model, a comparison is made with grand canonical Monte Carlo (GCMC) simulations by Mosher et al.17 While a direct comparison between this model and molecular simulations cannot be expected, this example serves as a way to qualitatively evaluate the results of the two methods. The chosen system is the adsorption of methane in an activated carbon slit pore shown in Figure 7. On a qualitative basis, the model is able to capture the overall trend predicted by GCMC; however, there is a clear deviation observed in the maximum peak density, which is approximately 515 kg/m3 for molecular simulations and around 343 kg/m3 for this model. It should be noted that the underlying physics of the two models are different. In particular, the highest

PSD(Hm) + PSD(Hm + 1) (37) 2 The amount adsorbed on each pore is determined at a specific bulk pressure. The following relation is used to find the total adsorbed amount of a pure component in r̂conf confined regions, Vm(Hm) = (Hm + 1 − Hm)

H



Article

ADDΓ (%) =

100 p̂

ADDNads (%) =

ADDx (%) =

p̂

∑ |(Γ exp, i − Γ calc,im)/Γ exp,i| p̂

100 p̂

100 p̂

(42)

i=1 exp , i exp , i calc, i − Nads )/Nads | ∑ |(Nads i=1

(43)

p̂ calc | ∑ |xex,expi − xex,im i=1

(44)

In Figure 8, the surface excess predictions for pure methane and nitrogen are reported for a temperature of 298 K and bulk

Figure 7. Methane local density profile in an activated carbon pore of H = 3 nm at 298 K, εCH4/kB = 148 K, bulk phase pressure of 1 MPa. GCMC calculations are taken from Mosher et al.17

achievable local density that can be predicted is limited by the packing density of the EOS, whose value includes void spaces between molecules. Such a limitation does not exist in the case of GCMC. Results with the DRA Potential. In this case, the DRA potential is used to perform adsorption calculations for CH4, CO2, N2, and their mixtures on activated carbon. The predictions are compared with experimental data from Dreisbach et al.,38 and the same data are used by Monsalvo and Shapiro39 to fit the DRA pure component solid−fluid interaction parameters for this system. Their MPTA model is based on the volume-translated Soave−Redlich−Kwong EOS with the DRA potential, and the energy interaction parameter values are 8143 J/mol, 7980 J/mol, and 6328 J/mol for CH4, CO2, and N2, respectively, for a solid with specific pore volume of 4.093 × 10−4 m3/kg. The volume-translated Peng−Robinson EOS was used here with binary interaction parameters set to zero for all mixtures. The adsorption data used in this example is available in terms of the excess surface amount at a given temperature and bulk phase pressure. The surface excess of each component i in region m, Γim, is computed by subtracting the bulk amount from the adsorbed amount, Γim =

Vm Hm

∫0

Hm

(ρ(z)xi(z) − ρB xBi) dz

Figure 8. Surface excess amounts for the adsorption of pure N2 and CH4 on activated carbon at 298 K. Solid lines are predictions by this work and points are experimental data by Dreisbach et al.38

phase pressures up to 6 MPa. The model predicts the adsorption isotherm of each pure component reasonably well, with larger deviations from experimental data at higher bulk phase pressures. The total and individual surface excess values are computed at the same temperature and bulk phase pressure range for binary and ternary mixtures. The experimental data for these mixtures were reported at different bulk compositions, and the calculations are performed at the reported compositions to obtain an accurate comparison. Figures 9, 10, and 11 show the results for selected binary mixtures, and Figure 12 shows the results for the ternary mixture of CH4−CO2−N2. In general, satisfactory results are obtained at low pressures, as in the case of pure components. The adsorption of CH4 is accurately predicted over the entire pressure range, while relatively larger deviations are observed in the cases of N2 and CO2. For instance, the highest average deviation value is reported for the CO2−N2 mixture. Table 3 lists the deviations between this work and experimental data for all the computed points at various compositions. It includes the ADDxex,i, and ADDNads values obtained by Monsalvo and Shapiro for comparison. The model performance is satisfactory, especially with respect to the adsorbed compositions values (maximum ADDxex,i is 8.7%). Results for Multiple Fields. Gravitational Field and Steele Potential Results. In this case, carbon slit-pores are used to emulate adsorption on a carbonaceous formation. The

(39)

In this work, eq 39 simplifies to lm

* ) − Vmρ xBi Γim = (∑ nijm B j=1

(40)

The individual excess amounts are added to find the total surface excess in region m, Γm. Subsequently, the excess mole fraction of component i is xex, im =

Γim Γ = c ̂ im Γm ∑k = 1 Γkm

(41)

The deviation between experimental data and calculations for each data point is computed, and the average deviation in Γ, Nads, and xex,im for each mixture (for p̂ data points) is computed according to the following relations, I



Article

Figure 12. Ternary adsorption isotherm of CH4, CO2 and N2 on activated carbon at 298 K, mean bulk composition of 0.7 and 0.1 for CH4 and CO2.38

Figure 9. Binary adsorption isotherm of CH4 and N2 on activated carbon at 298 K, mean CH4 bulk composition of 0.4. Empty symbols are predictions by this work and filled symbols are experimental data.38 Dashed lines are linear interpolations between the calculated points.

Table 3. Average Deviations for Binary and Ternary Mixtures Using DRA potential

mixture CH4−N2 CH4−CO2 CO2−N2 CH4−CO2− N2

Γ

xex,i

Nads

total surface excess (%)

mole fraction in the adsorbate (%)

absolute adsorbed amount (%)

this work

this work

MPTA

this work

MPTA

6.1 5.3 9.6 5.9

3.4 3.9 4.2 8.7

2.5 4.2 3.8 17.5

10.5 7.7 11.8 6.7

4.5 7.4 7.9 9.4

equilibrium molar distribution of a 5500 mol mixture of 80% mol CH4 and 20% mol CO2 is obtained for a reservoir depth of 4000 m at T = 305 K. The specified volumes of the bulk phase and pore region are equal to 1.0 m3 and 0.1 m3, respectively. The Steele potential is used to describe confinement in an activated carbon pore with H = 6 nm. The effect of the gravitational field, which is assumed to act on a system consisting of a bulk phase and a slit pore region, is also included in the equilibrium calculations. The slit pores are discretized in 50 layers in the horizontal direction. Additionally, the bulk and slit regions are discretized in 10 vertical layers. Therefore, the bulk and slit pore regions are modeled by 10 and 500 grid elements, respectively. Figure 13 displays the bulk pressure dependence with depth, which is an outcome of the calculation. The slight change in curvature indicates that the rate of pressure increase with depth becomes more pronounced at larger depths as a consequence of the increasing fluid densities. Figure 14 shows the effect of depth on the molar concentrations of CH4 and CO2 in the porous space. The molar concentration of CH4 is higher than that of CO2 in the pore region at all depths but the molar ratio CH4/CO2 is less than 4, which is the overall molar ratio in the simulated system. It is interesting to observe that the effect of gravity on the molar concentration of CO2 is more pronounced than on the molar concentration of CH4 because of its larger molar mass. Figure 15 shows the compositional profile of CO2 inside the slit pore region at h = −200 m and h = −3800 m. As the

Figure 10. Binary adsorption isotherm of CH4 and CO2 on activated carbon at 298 K, mean CH4 bulk composition of 0.9.38

Figure 11. Binary adsorption isotherm of CO2 and N2 on activated carbon at 298 K, mean CO2 bulk composition of 0.2.38

J



Article

profiles are symmetric, the plot spans the distance from the pore wall to its center. For all heights, CO2 has a high mole fraction near the solid wall, which then decreases until it reaches a near-constant value near the center of the pore. The CO2 profiles inside the pore are shifted to higher values at greater depths, which correlates well with the pronounced increase in CO2 molar concentration in the pore region displayed in Figure 14. The FORTRAN implementation of the procedure presented here runs this example in about 4 min of CPU time in a Lenovo ThinkPad Yoga 2015 model (equipped with an Intel processor model i7-4500U and 8 GB memory). The calculations for the other cases presented in this work converge in around 60 s. Refinements to the computer code to take advantage of symmetries could possibly decrease these times by 50% or more. Figure 13. Effect of depth on the bulk phase pressure at T = 330 K.

■

CONCLUSIONS

■

AUTHOR INFORMATION

A general framework for the computation of thermodynamic equilibrium in systems that have multiple regions subject to multiple external fields has been developed. The specifications are the temperature, component amounts, and volume of each region present. The underlying mathematical problem is the minimization of the Helmholtz energy of the system, subject to constraints that represent component mass balances and volume conservation equations applied to each region. Among other applications, the formulation and its solution procedure are particularly suitable for determining the equilibrium conditions in batch adsorption. The formulation is general but the focus of this work was on the compositional segregation in isothermal reservoirs due to gravity and the spatial segregation of components close to pore walls. The formulation was validated by comparing its results with experimental adsorption data of light gases on activated carbon and simulation results from the literature that use the Steele and DRA potentials to model such systems. The overall observation is that the formulation and solution procedure presented here provide results that are in very good agreement with literature results in most of the studied cases. The compositional predictions in a conceptual example that deals with the simultaneous effect of confinement and gravity follow expected trends for the simulated situation. Future work could include the application of various EOS, adsorption potentials, and external fields. In addition, the model can be extended to describe adsorption of pure components and mixtures confined in pores of spherical and cylindrical geometries. Furthermore, implementing a phase stability test in the equilibrium calculations will enable the prediction of phase transitions in both the bulk and confined regions.

Figure 14. Effect of depth on the molar concentration of CH4 and CO2, inside a slit pore, under the combined effect of gravity and Steele potential.

Corresponding Author

*E-mail: [email protected]. Funding

This publication was made possible by NPRP Grant No. 5-3442-129 and No. 7-042-2-021 from the Qatar National Research Fund (a member of the Qatar Foundation). The statements made herein are solely the responsibility of the authors.

Figure 15. Local mole fraction profiles of CO2 at h = −200 m and h = −3800 m, under the combined effect of gravity and Steele potential.

Notes

The authors declare no competing financial interest. K



■

Article

μijm = chemical potential of component i in layer j in region m ρB = molar density of the bulk phase ρsm = density of the solid (adsorbent) in region m σs,im = solid−fluid diameter

ACKNOWLEDGMENTS The authors gratefully acknowledge Prof. Abbas Firoozabadi and Dr. Zhehui Jin for providing additional details about their work. The senior authors (IGE and MC) of this work are grateful to Kenneth R. Hall for his collegiality, friendship and leadership at Texas A University at Qatar and the junior authors (ND and MLDL) for being a great teacher and inspiration for his students.

■

Subscripts

i,k,l = component j = layer m = region m = mixture (in eqs 16, 17, 19, and 20) ref = reference state

LIST OF SYMBOLS

Roman Letters

Superscripts

a = molar Helmholtz energy A = Helmholtz energy ares = molar residual Helmholtz energy Ajm *,in = internal Helmholtz energy of component i in layer j of region m A*jm,f = external Helmholtz energy of component i in layer j of region m Ajm = Helmholtz energy in layer j of region m a,b,c = Peng−Robinson EOS parameters ĉ = number of chemical components f m = number of external fields acting on region m g = molar gibbs energy ĝ = gravitational constant hjm = depth of layer j in region m Hm = total pore width in region m Hin,m = internal pore width in region m (excluding the adsorbent radius, σs) lm = number of layers in region m Mi = molar mass of component i Nads = absolute adsorbed amount n*ijm = number of moles of component i in all replicas of layer j in region m ni = total number of moles of component i p̂ = number of data points P = bulk phase pressure P0 = reference pressure Pin = pressure calculated from an equation of state r̂ = number of regions r̂conf = number of confined regions rm = number of replicas in region m R = universal gas constant si = volume shift parameter of component i T = system temperature V = total volume of the system Vm = total volume of a region v = molar volume vjm = molar volume of layer j in region m v̂0 = maximum porous volume xB,i = mole fraction of component i in the bulk phase xijm = mole fraction of component i in layer j in region m zjm = distance from the center of the pore to the confining wall

■

* = total amounts of all replicas f = external field contribution fw = far wall of the slit pore ig = ideal gas property in = internal field contribution

REFERENCES

(1) Espósito, R. O.; Castier, M.; Tavares, F. W. Calculations of Thermodynamic Equilibrium in Systems Subject to Gravitational Fields. Chem. Eng. Sci. 2000, 55, 3495−3504. (2) Lira-Galeana, C.; Firoozabadi, A.; Prausnitz, J. M. Computation of Compositional Grading in Hydrocarbon Reservoirs. Application of Continuous Thermodynamics. Fluid Phase Equilib. 1994, 102, 143− 158. (3) Galliero, G.; Montel, F. Nonisothermal gravitational segregation by molecular dynamics simulations. Phys. Rev. E 2008, 78, 041203. (4) Høier, L.; Whitson, C. H. Compositional Grading - Theory and Practice. Soc. Pet. Eng. J. 2001, 4, 525−535. (5) Hamoodi, A. N.; Abed, A. F.; Firoozabadi, A. Compositional modelling of two-phase hydrocarbon reservoirs. J. Can. Pet. Technol. 2001, 40, No. PETSOC-01-04-03, DOI: 10.2118/01-04-03. (6) Montel, F.; Gouel, P. Prediction of Compositional Grading in a Reservoir Fluid Column. SPE Annu. Technol. Conf. Exhib. 1985, DOI: 10.2118/14410-MS. (7) Castier, M.; Espósito, R. O.; Tavares, F. W.; Peçanha, R. P. Calculation of Sedimentation Equilibrium using a Modified Flash Algorithm. Chem. Eng. Sci. 2001, 56, 3771−3779. (8) Cohen, K. The Theory of Isotope Separation; McGraw-Hill: New York, USA, 1951. (9) Cracknell, R.; Golombok, M. Monte-Carlo Simulations of Centrifugal Gas Separation. Mol. Simul. 2004, 30, 501−506. (10) Golombok, M.; Morley, C. Thermodynamic Factors Governing Centrifugal Separation of Natural Gas. Chem. Eng. Res. Des. 2004, 82, 513−516. (11) Golombok, M.; Chewter, L. Centrifugal Separation for Cleaning Well Gas Streams. Ind. Eng. Chem. Res. 2004, 43, 1734−1739. (12) Castier, M.; Tavares, F. W. Centrifugation Equilibrium of Natural Gas. Chem. Eng. Sci. 2005, 60, 2927−2935. (13) Martins, L. S. F.; Tavares, F. W.; Pecanha, R. P.; Castier, M. Centrifugation Equilibrium for Spheres and Spherocylinders. J. Colloid Interface Sci. 2005, 281, 360−367. (14) Tsori, Y.; Leibler, L. Phase-Separation in Ion-Containing Mixtures in Electric Fields. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 7348−7350. (15) Canas-Marin, J. D.; Ortiz-Arango, W. A.; Guerrero-Aconcha, U. E.; Lira-Galeana, C. Thermodynamics of Wax Precipitation under the Influence of Magnetic Fields. AIChE J. 2006, 52, 2887−2897. (16) Canas-Marin, W. A.; Lira-Galeana, C. Calculation of Critical Points in Gas-Condensate Mixtures under the Influence of Magnetic Fields. Ind. Eng. Chem. Res. 2010, 49, 7610−7619. (17) Mosher, K.; He, J.; Liu, Y.; Rupp, E.; Wilcox, J. Molecular simulation of methane adsorption in micro- and mesoporous carbons with applications to coal and gas shale systems. Int. J. Coal Geol. 2013, 109−110, 36−44.

Greek Letters

β = solid heterogeneity parameter, takes on a value of 2 for activated carbon Γ = surface excess δ = Kronecker delta function Δm = interlayer spacing of the adsorbent εs,im = solid−fluid energy interaction parameter ε0i = characteristic energy of component i θijm * = distribution factor L



Article

(18) Alvarez, M.; Levesque, D.; Weis, J. J. Monte Carlo Approach to the Gas-Liquid Transition in Porous Materials. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1999, 60, 5495−5504. (19) Frenkel, D.; Smit, B. Understanding Molecular Simulations: from Algorithms to Applications; Academic Press: San Diego, USA, 2002. (20) Wu, J. Z. In Molecular Thermodynamics of Complex Systems; Lu, X., Hu, Y., Eds.; Springer: Berlin, 2009; pp 1−74. (21) Peng, B.; Yu, Y. A. Density Functional Theory for LennardJones Fluids in Cylindrical Pores and its Applications to Adsorption of Nitrogen on MCM-41 Materials. Langmuir 2008, 24, 12431−12439. (22) Travalloni, L.; Castier, M.; Tavares, F. W.; Sandler, S. I. Thermodynamic Modeling of Confined Fluids using an Extension of the Generalized van der Waals Theory. Chem. Eng. Sci. 2010, 65, 3088−3099. (23) Travalloni, L.; Castier, M.; Tavares, F. W. Phase Equilibrium of Fluids Confined in Porous Media from an extended Peng-Robinson Equation of State. Fluid Phase Equilib. 2014, 362, 335−341. (24) Tan, S. P.; Piri, M. Equation-of-State Modeling of ConfinedFluid Phase Equilibria in Nanopores. Fluid Phase Equilib. 2015, 393, 48−63. (25) Shapiro, A. A.; Stenby, E. H. Potential Theory of Multicomponent Adsorption. J. Colloid Interface Sci. 1998, 201, 146−157. (26) Luo, S.; Lutkenhaus, J. L.; Nasrabadi, H. Effect of Confinement on the Bubble Points of Hydrocarbons in Nanoporous Media. AIChE J. 2016, 62, 1772−1780. (27) Peng, D.; Robinson, D. B. S. New 2-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59−64. (28) Jhaveri, B.; Youngren, G. Three-Parameter Modification of the Peng-Robinson Equation of State To Improve Volumetric Predictions. SPE Reservoir Eng. 1988, 3, 1033−1040. (29) Steele, W. A. The Physical Interaction of Gases with Crystalline Solids I. Gas-Solid Energies and Properties of Isolated Adsorbed Atoms. Surf. Sci. 1973, 36, 317−352. (30) Monsalvo, M. A.; Shapiro, A. A. Study of High-Pressure Adsorption from Supercritical Fluids by the Potential Theory. Fluid Phase Equilib. 2009, 283, 56−64. (31) Li, Z.; Jin, Z.; Firoozabadi, A. Phase Behavior and Adsorption of Pure Substances and Mixtures and Characterization in Nanopore Structures by Density Functional Theory. SPE J. 2014, 19, 1096− 1109. (32) Dubinin, M. M. Generalization of the Theory of Volume Filling of Micropores to Nonhomogeneous Microporous Structure. Carbon 1985, 23, 373−380. (33) Dubinin, M. M.; Astakhov, V. A. Adsorbents, Development of the Concept of Volume Filling of Micropores in the Adsorption of Gases and Vapors by Microporous Adsorbents. Bull. Acad. Sci. USSR, Div. Chem. Sci. 1971, 20, 8−11. (34) Murray, W. Numerical Methods for Unconstrained Optimization; Academic Press: London, 1972; pp 57−71. (35) Michelsen, M. The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1−19. (36) Qiao, S.; Wang, K.; Hu, X. Using Local IAST with Micropore Size Distribution To Predict Multicomponent Adsorption Equilibrium of Gases in Activated Carbon. Langmuir 2000, 16, 1292−1298. (37) Li, Z.; Firoozabadi, A. Interfacial Tension of Nonassociating Pure Substances and Binary Mixtures by Density Functional Theory Combined with Peng-Robinson Equation of State. J. Chem. Phys. 2009, 130, 154108. (38) Dreisbach, F.; Staudt, R.; Keller, J. U. High Pressure Adsorption Data of Methane, Nitrogen, Carbon Dioxide and their Binary and Ternary Mixtures on Activated Carbon. Adsorption 1999, 5, 215−227. (39) Monsalvo, M. A.; Shapiro, A. A. Modeling Adsorption of Binary and Ternary Mixtures on Microporous Media. Fluid Phase Equilib. 2007, 254, 91−100.

M


Phase Equilibrium with External Fields: Application to Confined Fluids

Recommend Documents