Diffusion of hard sphere fluids in disordered porous media: Enskog theory description

,


Introduction
The properties of fluids confined in disordered porous media substantially differ from those of bulk fluids [1]. For the last three decades, starting from the pioneering work of Madden and Glandt [2], great theoretical efforts have been devoted to the study of fluid adsorption in disordered porous media. So far most of these theoretical and simulation-based studies have been focused on the static structural and thermodynamic properties and only a few of them are devoted to investigations of dynamic properties [3][4][5]. According to the model proposed by Madden and Glandt [2], a disordered porous medium is presented as a matrix of quenched configurations of randomly distributed particles. The specificity of such an approach is connected with the double quenched-annealed averages: the annealed average is taken over all fluid configurations while the additional quenched average should be taken over all realizations of the matrix. A standard approach to solve this problem is based on the replica method. Using the replica Ornstein-Zernike (ROZ) integral equation theory [6], the statistical mechanics approach of liquid state was extended to a description of different models of fluids confined in random porous media [7,8] including the chemical reacting fluids adsorbed in porous media [9,10]. However, unlike bulk fluids, no analytical results have been obtained in this approach even for the simplest model such as a hard sphere fluid in a hard-sphere matrix.
In order to solve this problem, Holovko and Dong [11] proposed to extend the classical scaled particle theory (SPT) [12,13] for the description of thermodynamic properties of hard sphere fluids in disordered porous media. During the last decade the SPT approach for hard sphere fluids in disordered porous media was essentially improved and developed [14][15][16][17][18][19]. The approach proposed in [11] and named SPT1, contains a subtle inconsistency appearing when the size of matrix particles is considerably larger than the size of fluid particles. This inconsistency was eliminated in a new approach named SPT2 [16]. Consequently, the first rather accurate analytical expressions were obtained for the chemical potential and pressure of a hard sphere fluid confined in a hard sphere (HS) or of an overlapping hard sphere (OHS) matrix. The obtained expressions include three parameters defining the porosity of the matrix. The first one is related to the bare geometry of the matrix. It is the so-called geometric porosity φ 0 characterizing the free volume, which is not occupied by matrix particles. The second parameter φ is defined by the chemical potential of a fluid in the limit of infinite dilution. It is the so-called probe particle porosity characterizing the adsorption of a fluid in an empty matrix. Usually, φ φ 0 . The third parameter φ * is defined by the maximum value of the fluid packing fraction of a hard sphere fluid in a porous medium. It characterizes the maximum adsorption capacity of a matrix for a given fluid. It is responsible for crowding effects of a fluid in porous media. In [18], a general expression for φ * was proposed which is exact for one-dimensional case and can be considered a good approximation for higher dimensions.
The developed SPT2 approach was generalized for fluids of anisotropic particles [20,21], for a hard sphere mixture [22] and for a mixture of hard sphere and anisotropic particles [23,24] in disordered porous media. The obtained results for the models considered were also used as the reference system taking into account attractive [24,25], associative [26] and Coulombic [27][28][29] interactions.
The SPT2 approach is also useful for the description of static structure of hard sphere fluids in disordered porous media. In particular, in [26] a simple and rather accurate expression was obtained for the contact value of the fluid-fluid pair distribution function which was used for the description of phase behaviour and percolation properties of the patchy colloidal fluids in disordered porous media. The structural information provided by the SPT2 approach similar to the bulk case can be used as a source of input in the Enskog theory for the description of transport properties of hard sphere fluids in porous media (for the bulk case see the reference [30] and references therein for details).
The first attempt to apply the Enskog theory for a hard sphere fluid in a hard-sphere matrix was considered by Yethiraj with coworkers [3] for the self-diffusion coefficient. They considered a fluid in a disordered matrix as a mixture of two components one of which is quenched in space and was treated as the particles with the infinite mass. In such an approach, the contact values of fluid-fluid and fluidmatrix pair-distribution functions are introduced as the input of theory. In [3], they were taken from computer simulations. The comparison between theoretical and simulation results shows a systematic overestimation of theoretical predictions which increase with the increasing fraction of matrix particles. It is clear that in order to improve the theoretical predictions we should modify the contact values of pair distribution functions. In [26], the contact value of the fluid-fluid pair distribution function was obtained in a simple analytical form and was found to be in good agreement with computer simulations. A similar expression can also be found for the contact value of the fluid-matrix pair distribution function. However, in contrast to thermodynamic properties, the obtained expressions for the contact values of both pair distribution functions are described only by the geometric porosity φ 0 and do not include the dependence either on the probe particle porosity φ or on the porosity φ * . We show in this paper that the application of such contact values neglects the effects of fluid particles being trapped by the matrix, and at least the probe particle porosity φ should be included in the Enskog theory for a correct description of the matrix effect.
In this report, we propose such an improvement in a manner similar to the description of thermodynamic properties. We show that such semi-empirical improvement of the Enskog equation predicts correct trends for the effect of porous media on the diffusion coefficient of hard sphere fluids in disordered porous media. The theory developed is applied to the investigation of the diffusion of hard sphere fluids  in different porous media, namely to a hard sphere matrix and to an overlapping hard sphere matrix. The effects of fluid density, matrix porosity, matrix morphology, and fluid to matrix sphere size ratio on the self-diffusion coefficient will be discussed. We also present some comparison with computer simulation data obtained by the group of Yethiraj [3].

The Enskog theory
The Enskog theory is based on the assumption that each collision between hard spheres is completely independent and instantaneous [30]. In theory, only binary collision is considered and all higher multiple collisions are neglected. Moreover, it is assumed that the frequency of the binary collision increases by an amount proportional to the probability of one particle finding its neighbours. The motion of particles of type µ in the fluid can be described in terms of the velocity autocorrelation function where ϑ µ (t) is the velocity of the particle type µ at time t and . . . is the equilibrium ensemble average. The time evolution of ψ µ (t) is well described by the generalized Langevin equation [31], where the role of memory kernel is played by the friction coefficient ξ µ (t). The frequency dependent ξ µ (t) is given by the Green-Kubo formula [31] ξ µ (z) = 1 3kT where F µ (t) is the force between a fixed particle of type µ and the surrounding particles at a time t, k is the Boltzmann constant, T is the absolute temperature. We also note that F µ (t)F µ (0) is not the usual force-force time correlation function since the time evolution of F(t) is given by a Liouville projection operator where iL N is the Liouville operator for an N-particle system, Q = 1 − P and P is the Mori-Zwanzig projection operator defined as [31] The self-diffusion coefficient D µ for particles of type µ is related to the corresponding friction coefficient ξ µ via the Einstein relation as where ξ µ = ξ µ (z = 0) is the friction coefficient in the stationary limit. For a mixture of hard spheres, the Enskog theory leads to the following expression [3,32] where m µ , m ν are the masses of particles of types µ and ν, respectively, η ν = 1/6πρ ν σ 3 νν is the packing fraction of particles of type ν, ρ ν is the number density, is the hard sphere interaction diameter between species µ and ν, and g µν (σ µν ) is the pair distribution function at contact between species µ and ν. Similar to [3], we mimic a hard sphere fluid in a hard sphere matrix by a binary mixture where one component is infinitely massive. As a result, we have

23605-3
Hereafter, we use the conventional notations [6][7][8], where the index "1" is used to denote the fluid component and the index "0" denotes the matrix particle, τ = σ 11 /σ 00 . In accordance with (2.5), we have the following expression for self-diffusion of a hard sphere fluid in a disordered porous medium where D 0 1 = (kT σ 2 11 /m 1 ) 1/2 . For further calculations, the expressions for the contact values of the pair distribution functions g 11 (σ 11 ) and g 10 (σ 10 ) are needed. To this end, we develop and improve the SPT2 approach presented in [26]. We start from g 11 (σ 11 ). In accordance with the SPT2 approach, the contact value of a small scaled particle and a fluid particle can be presented in the form where σ 1s = 1/2(σ 11 + σ ss ), λ s = σ ss /σ 11 , p 0s (λ s ) is the probability of finding a cavity created by the scaled particle in the matrix in the absence of fluid [17]. Here, the index "s" is used for the scaled particle. For the hard sphere matrix In order to obtain the expression for g 11 (σ 11 ), we expand g 1s (σ 1s ) as where G (0) 1s , G (1) 1s and G (2) 1s are found from the continuity of g 1s (σ 1s ) and the first and second derivatives with respect to λ s at λ s = 0. Subsequently, we can put in (2.12) λ s = 1. As a result, we have where φ 0 = 1 − η 0 is the geometric porosity. We note that in the calculation G (2) 1s the derivative was taken only from the dominator of the first derivative from (2.10) and the numerical coefficient 9 4 was changed to 1 2 in order to describe correctly Carnahan-Starling correction [21]. The expression for the contact value g 10 (σ 10 ) can be found in a similar manner and can be presented in the form (2.14) Figure 1 depicts the prediction from the Enskog theory versus computer simulation data taken from [3] for the dependence of the self-diffusion coefficient D 1 on the fluid packing fraction η 1 for the case τ = 1 at different values of the packing fraction of matrix particles η 0 . As we can see, with an increasing η 1 or η 0 , the self-diffusion coefficient D 1 monotonously decreases. We have a very good agreement between the theory and computer simulations for the bulk case and for a rather low value η 0 = 0.05. However, similar to [3] figure 1 demonstrates that with an increasing η 0 , the theory greately overestimates the value of D 1 by an order of magnitude for η 0 = 0.20. Yethiraj with coworkers [3] noted that such a discrepancy between the Enskog theory prediction and computer simulation results demonstrates that with an increasing η 0 , configurations of the fluid and the disordered matrix are very different from those of the equilibrium mixtures, which is the basic assumption in the Enskog theory. We remark that in the equilibrium case, the presence of matrix particles hinders the dynamics of fluids only by the geometric porosity φ 0 = 1 − η 0 , which according to (2.13) and (2.14) defines the contact values g 11 (σ 11 ) and g 10 (σ 10 ). However, in the case of immobile particles of the matrix, there are additional effects from the matrix due to the geometric constrains induced by the obstacles, which is not present in the equilibrium case. These additional effects become significant with an increasing η 0 , slowing down the translation diffusion of fluids. In this paper we show that this effect is strongly connected with the probe particle porosity φ introduced by us for the description of thermodynamic properties of hard sphere fluids [16,17], which, however, is not present in the description of contact values of the pair distribution functions g 11 (σ 11 ) and g 10 (σ 10 ).

Revision and extension of the Enskog theory to hard sphere fluids in disordered porous media
We remember that the original Enskog equation for a hard sphere fluid was formulated nearly a hundred years ago as the generalization of the Boltzmann kinetic equation to high densities of hard sphere fluids [31,33]. It includes, as a multiplier the pressure term (p 1 /kT ρ 1 − 1), which due to the virial theorem [34] p 1 kT ρ 1 − 1 = 4η 1 g 11 (σ 11 ) (3.1) can be presented in a well-known form [30,32] via the contact value of the fluid-fluid distribution function g 11 (σ 11 ), while the expression for self-diffusion of a hard sphere fluid can be presented in the well-known form [η 1 g 11 (σ 11 )] −1 .

(3.2)
However, for a hard sphere fluid in a porous medium, no simple virial expressions like (3.1) for the pressure p 1 exist and it is somewhat problematic to write an expression like (2.9) for self-diffusion of a hard sphere fluid in a disordered porous medium. Of course, we can write the expression (2.9) for the D 1 /D 0 1 but in this expression we cannot consider g 11 (σ 11 ) and g 10 (σ 10 ) as the contact values of the fluid-fluid and fluid-matrix distribution functions. In general, they are some thermodynamic properties defined by equation (2.9). In this paper, we do not modify the expression (2.13) for g 11 (σ 11 ). We consider and modify only the expression (2.14) for g 10 (σ 10 ). We start from the infinite dilution. When η 1 → 0 In the equilibrium case, according to (2.14) (3.4) In this paper, we change the expression (3.4) to

23605-5
which can be considered as the ratio of the probability of finding a cavity created by a fluid particle in an empty matrix to the probability of finding a cavity created by a point particle in an empty matrix. φ 0 and φ are the geometric porosity and the probe particle porosity, respectively. In accordance with [17], for a hard sphere fluid in a hard sphere matrix In the presence of fluid particles, similar to thermodynamic properties [17][18][19], we put (3.7) The first term in (3.7) leads to a divergence at η 1 = φ and similar to thermodynamic consideration, we change 1/(1 − η 1 /φ) to [17][18][19] (3.8) Consequently, for g 10 (σ 10 ), we have Now, for the diffusion coefficient we again have the expression (2.9) in which, however, g 10 (σ 10 ) is given by (3.9) and for g 11 (σ 11 ) we have the previous form (2.13). The results obtained from this expression for τ = 1 are presented in figure 2. As we can see, the improved version of the Enskog theory is much better than the standard version (see figure 1). It agrees with the computer simulation data and quite accurately reproduces the trend of change in the dependence of matrix particles on the packing fraction.
In figure 3, we compare the prediction for the ratio D 1 /D 1 (η 0 = 0) as a function of η 0 against computer simulations results from [3], where D 1 (η 0 = 0) is the diffusion of the fluid at the same value of η 1 as for D 1 but with η 0 = 0. As we can see, after extension, the Enskog theory reproduces the trends of the computer simulation results. The dependence of D 1 on η 0 is more or less similar for all values η 1 . Only our theory does not reproduce the effects of underestimation of D 1 /D 1 (η 0 = 0) for η 1 = 0.05 observed in  computer simulations [3] and slightly overestimates the values D 1 /D 1 (η 0 = 0) at a higher fluid density η 1 = 0.20. Now, we discuss the effect of the fluid to matrix sphere size ratio τ on the self-diffusion coefficient of a hard sphere fluid in a disordered porous medium. Figure 4 illustrates the diffusion coefficient of a hard sphere fluid in a disordered porous medium normalized by the diffusion coefficient of hard spheres of the same size in the absence of the matrix, D 1 /D 1 (η 0 = 0), as a function of η 1 for different τ. Similar to [3], the total packing fraction η = η 1 + η 0 is fixed at η = 0.20 for all cases. For comparison, in figure 4 the computer simulation data taken from [3] are presented as well. As we can see, theoretical and simulations results are different. In all cases, the theoretical predictions are lower than the computer simulations data. However, both approaches qualitatively show the same trend of dependence on the size ratio which is opposite to the theoretical results obtained from the standard Enskog theory [3]. In the limiting case, when τ → 0, D 1 → D 1 (η 1 /φ 0 ) and the ratio D 1 /D 1 (η 0 = 0) cannot be larger than unity contrary to the simulation data [3]. In the opposite case, when τ → ∞ and η 0 → 0 in such a way that τ 3 η 0 = 1/6πρ 0 σ 3 11 is fixed [19], we have As we can see, the ratio D 1 /D 1 (η 0 = 0) decreases with a decreasing η 1 and with an increasing τ. As noted by Yethiraj with coworkers [3], the size-dependent behaviour of diffusion of a fluid in a porous medium can be qualitatively explained using the concept of free volume [35]. Replacing a big fluid particle by several immobile matrix particles of smaller sizes increases the excluded volume for the particles and this slows down the translation motion of the fluid particles. In the opposite case, replacing small fluid hard spheres by an immobile particle of a bigger size decreases the excluded volume for the diffusing particles. As a result, diffusion of fluid particles becomes faster.
Finally we discuss the effect of morphology of a porous medium on the behaviour of self-diffusion coefficient of hard sphere fluids. To this end, we consider a hard sphere fluid in two different matrices, namely a hard sphere matrix and an overlapping hard sphere matrix. For the self-diffusion coefficient of a hard sphere fluid in a porous medium, we use the expression (2.9), where g 11 (σ 11 ) and g 10 (σ 10 ) are given by the expression (2.13) and (3.9), respectively. The porosities φ and φ 0 for a hard sphere matrix are given in (3.6) and for the overlapping hard sphere matrix in [17] φ = exp −(1 + τ) 3 η 0 and φ 0 = exp(−η 0 ). (3.12) In figure 5, we compare the self-diffusion coefficients D 1 of a hard sphere fluid as functions of η 1 in these two different matrices at different η 0 . As we can see, the self-diffusion coefficient D 1 in an overlapping hard sphere matrix at the same η 0 and η 1 for all cases is higher than in a hard sphere matrix and this difference increases with an increasing value of parameter η 0 .

Conclusions
In this paper, we extended the Enskog theory to hard sphere fluids in disordered porous media and used it for the description of the self-diffusion coefficient of hard spheres. At the beginning, we modelled a hard sphere fluid in a disordered porous medium by an equilibrium two-component mixture in the limit m 0 → ∞ for the species corresponding to the matrix particles. In such an approach, the contact values of the fluid-matrix and fluid-fluid pair distribution functions are introduced as the input of the Enskog theory. For their calculation, we used the scaled particle theory previously extended by us in [11,[14][15][16][17][18][19] for the description of thermodynamic properties of hard sphere fluids in disordered porous media. The expressions obtained for the contact values of the fluid-matrix and fluid-fluid pair distribution functions are described only by the geometric porosity φ 0 and do not include the dependence on the probe particle porosity φ and porosity φ * defined by the maximum value of the fluid packing fraction of a hard sphere fluid in a porous medium. All three types of porosity are important for the description of thermodynamic properties of a hard sphere fluid in a disordered medium. We showed that the application of such contact values neglects the effects of trapping of fluid particles by the matrix and at least the probe particle porosity φ should be included in the Enskog theory for a correct description of the matrix effect. Since the contact values of the fluid-fluid and fluid-matrix pair distribution functions have proven to be in very good agreement with computer simulations data [26], we consider that the Enskog theory or any other theory based only on static correlations between particles cannot correctly describe the effect of a static disordered matrix. Such a conclusion was also reached by Yethiraj with coworkers [3] from the analysis of the computer simulations data.

23605-8
In the present paper, we extend the Enskog theory by modifying the contact values of the fluidmatrix and fluid-fluid pair distribution functions with new properties which include the dependence not only on the geometric porosity φ 0 but also on the probe particle porosity φ. In this procedure, we consider that in the limit φ → φ 0 these fictitious contact values coincide with the corresponding real contact values obtained by us in the framework of the scaled particle theory. We should note that the correction of the fluid-matrix contact value and the correction of the fluid-fluid contact value play different roles. The correction of the fluid-matrix contact value is very important for the description of the self-diffusion coefficient at small fluid densities and the correction of the fluid-fluid contact value can be important for the correction of the density dependence of the self-diffusion coefficient. Due to this, in this paper as the first step of such correction we consider only the modification of the fluid-matrix contact value. In the infinite dilution limit η 1 → 0, we have the expression (3.5) for the fluid-matrix contact value. Subsequently, at a finite value of η 1 , we have the expression (3.7) which corresponds to the approximation SPT2b for the description of thermodynamic properties [16,17]. Then, we use the approximation (3.8) which leads to the expression (3.9) corresponding to the approximation SP2b1 in the description of thermodynamic properties [17][18][19]. In order to describe the fluid-fluid contact value, we left the expression (2.13) without any modification. We will discuss the possibility of modification of g 11 (σ 11 ) in a separate paper.
We showed that such a semi-empirical improvement of the Enskog theory predicts the correct trends for the influence of a porous medium on the diffusion coefficient of a hard sphere fluid. We obtained a good agreement with computer simulations data. We discussed the effects of fluid density, fluid to matrix sphere size ratio, matrix porosity and matrix morphology on the self-diffusion coefficient. We also plan to consider some other transport coefficients such as shear and bulk viscosities. In addition, we intend to generalize the results obtained to binary hard sphere mixtures in porous media.