6482
Langmuir 2004, 20, 6482-6489
Modeling Mercury Porosimetry Using Statistical Mechanics F. Porcheron and P. A. Monson* Department of Chemical Engineering, University of Massachusetts, Amherst, Massachusetts 01003-9303
M. Thommes Quantachrome Instruments, 1900 Corporate Drive, Boynton Beach, Florida 33426 Received January 7, 2004. In Final Form: April 17, 2004 We consider mercury porosimetry from the perspective of the statistical thermodynamics of penetration of a nonwetting liquid into a porous material under an external pressure. We apply density functional theory to a lattice gas model of the system and use this to compute intrusion/extrusion curves. We focus on the specific example of a Vycor glass and show that essential features of mercury porosimetry experiments can be modeled in this way. The lattice model exhibits a symmetry that provides a direct relationship between intrusion/extrusion curves for a nonwetting fluid and adsorption/desorption isotherms for a wetting fluid. This relationship clarifies the status of methods that are used for transforming mercury intrusion/ extrusion curves into gas adsorption/desorption isotherms. We also use Monte Carlo simulations to investigate the nature of the intrusion and extrusion processes.
I. Introduction The very wide range of uses for porous materials makes determination of their microscopic structure from experimental methods an important problem. Two important techniques for characterization of porous materials, nitrogen adsorption and mercury porosimetry experiments,1-4 rely on the principle of a fluid penetrating the solid void space as a function of a variable external pressure. Given the nature of the fluid and the way it interacts with the solid, the amount of fluid contained in the material as a function of pressure should reflect the solid microstructure. In combination gas adsorption and mercury porosimetry probe the porous material structure over a wide range of pore diameters. At lower pressures, the mercury penetrates the large pores of the material, whereas the smaller ones are filled first during gas adsorption experiments. The pore sizes probed by mercury porosimetry are in the mesopore to macropore range (4 nm to 400 µm) while in gas adsorption this range goes from micropores to mesopores (0.3-300 nm). During gas adsorption experiments, adsorption/desorption isotherms are obtained by monitoring the quantity of adsorbed gas as a function of the relative gas pressure P/P0, where P0 is the bulk saturation pressure.1 Since nitrogen or xenon are wetting fluids for most solid materials, the filling of the porous material occurs before the saturation pressure is reached. The relationship between the amount of adsorbed gas and the pore radius distribution is then obtained via a variety of methods such as application of the Kelvin equation.1 A great deal of molecular modeling research has been focused on understanding these experiments,5-14 and significant progress has been made. * To whom correspondence may be addressed. E-mail:
[email protected]. (1) Everett, D. H. Adsorption Hysteresis. In The Solid-Gas Interface; Flood, E. A., Ed.; Dekker: New York, 1967; Vol. 2. (2) Lowell, S.; Shields, J. E. Powder Surface Area and Porosity, 3rd ed.; Chapman and Hall: New York, 1991. (3) Leo´n y Leo´n, C. Adv. Colloid Interface Sci. 1998, 76-77, 341. (4) Van Brakel, J.; Modry, S.; Svata, M. Powder Technol. 1981, 29, 1. (5) Gelb, L.; Gubbins, K. E. Langmuir 1998, 14, 2097. (6) Gelb, L.; Gubbins, K. E. Langmuir 1999, 15, 305.
Mercury porosimetry is built on an observation made by Washburn that the structure of porous solids could be characterized by forcing a nonwetting liquid (mercury) to penetrate their pores.3 The volume of mercury penetrating the pores can be measured as a function of the applied hydraulic pressure, Ph, from which intrusion/extrusion curves are obtained. The connection to the pore radius distribution is usually made by the Washburn equation
Phr ) -2γ cos θ
(1.1)
where r is the pore radius, γ is the surface tension of mercury, and θ is the contact angle between the mercury and the solid pore. The Washburn equation has the same theoretical status as the Kelvin equation. It is similarly based on the pressure difference developed across a curved vapor-liquid interface and the assumption that the fluids in the pores have essentially bulklike properties. In mercury porosimetry experiments several specific features usually appear.2-4 Due to the nonwetting character of the mercury toward most solid material, the applied hydraulic pressure must be relatively high for the mercury to penetrate the pores. Typically porosimeters are designed to reach pressures of 60000 psi (=4000 bar). Hysteresis is found between the intrusion and extrusion process; i.e., for the same applied hydraulic pressure, the volume of mercury in the porous sample is greater during the extrusion of the mercury meniscus than during intrusion. Also in many cases the extrusion curve does not rejoin the intrusion one as the pressure is lowered. This feature is usually attributed to the entrapment of some mercury within the pore network of the solids.2 The amount of entrapped mercury can sometimes reach =90% of the total volume of the intruded mercury.3 Finally, the (7) Pellenq, R. J.-M.; Rousseau, B.; Levitz, P. E. Phys. Chem. Chem. Phys. 2001, 3, 1207. (8) Woo, H. J.; Sarkisov, L.; Monson, P. A. Langmuir 2001, 17, 7472. (9) Sarkisov, L.; Monson, P. A. Phys. Rev. E 2001, 65, 011202. (10) Woo, H. J.; Monson, P. A. Phys. Rev. E 2003, 67, 041207. (11) Sarkisov, L.; Monson, P. A. Langmuir 2001, 17, 7600. (12) Sarkisov, L.; Monson, P. A. Langmuir 2000, 16, 9857. (13) Kierlik, E.; Monson, P. A.; Rosinberg, M. L.; Sarkisov, L.; Tarjus, G. Phys. Rev. Lett. 2001, 87, 055701. (14) Nicholson, D. J. Chem. Soc., Faraday Trans. 1 1975, 71, 238.
10.1021/la049939e CCC: $27.50 © 2004 American Chemical Society Published on Web 06/16/2004
Modeling Mercury Porosimetry
Figure 1. Two models of intrusion of the mercury within a mesoporous pore network: selective (left) and continuous intrusion (right) (see text).
high hydraulic pressure applied on the material during the intrusion can sometimes lead to the partial collapsing of the pore network or alteration of the sample during the experiments.15 On the basis of a theoretical model of independent cylinders, parameters describing the pore structure of the solid sample can be calculated from the data obtained experimentally.3 For instance, the pore size distribution is directly obtained by inverting the intrusion/ extrusion curves via the Washburn equation. The total pore volume is the total volume of intruded mercury at high pressure providing that the applied pressure is high enough and the total surface area is calculated by performing an integration over the intrusion curve. Several other parameters can be derived from a set of intrusion/extrusion curves.3 Mercury porosimetry experiments are in principle a powerful technique for the characterization of porous solid structure and have been indeed widely used on a variety of materials.3 Various models and theories have been employed to explain intrusion/extrusion hysteresis although not thus far from a molecular perspective. The network and percolation theories16-26 consider hysteresis as to be due to network effects, which reflect the shape of the pores, the distribution of pore sizes, and their connectivity in the porous medium. In contrast the contact angle hysteresis27,28 model attributes hysteresis to different contact angles during intrusion and extrusion. In the so-called energy barrier model29 hysteresis is attributed to energy barriers associated with the formation of a liquid/vapor interface during mercury extrusion. Several questions remain about the nature of the microscopic process of mercury intrusion/extrusion in a mesoporous material. For instance does the liquid meniscus propagate progressively into the network from the bulk liquid or can we observe a selective intrusion of mercury in larger pores at low pressures following the prediction of the Washburn equation (see Figure 1)? In this paper, we model mercury porosimetry experiments with an approach based on statistical thermody(15) Pirard, R.; Blacher, S.; Brouers, F.; Pirard, J. P. J. Mater. Res. 1995, 10, 2114. (16) Orr, J. C. Powder Technol. 1970, 3, 117. (17) Rigby, S. P. J. Colloid Interface Sci. 2000, 224, 382. (18) Rigby, S. P.; Gladen, L. F. Chem. Eng. Sci. 2000, 55, 5599. (19) Rigby, S. P.; Fletcher, R. S.; Riley, S. N. J. Colloid Interface Sci. 2001, 240, 190. (20) Rigby, S. P.; Fletcher, R. S.; Riley, S. N. Ind. Eng. Chem. Res. 2002, 41, 1205. (21) Rigby, S. P.; Fletcher, R. S.; Riley, S. N. Appl. Catal., A 2003, 8526, 1. (22) Rigby, S. P.; Edler, K. J. J. Colloid Interface Sci. 2002, 250, 175. (23) Matthews, G. P.; Ridgway, C. J.; Spearing, M. C. J. Colloid Interface Sci. 1995 171, 8. (24) Salmas, C.; Androutsopoulos, G. J. Colloid Interface Sci. 2001, 239, 178. (25) Zgrablich, G.; Mendioroz, S.; Daza, L.; Pajares, J.; Mayagoitia, V.; Rojas, F.; Conner, W. C. Langmuir 1991, 7, 779. (26) Vidales, A. M.; Miranda, E.; Nazzaro, M.; Mayagotia, V.; Rojas, F.; Zgrablich, G. Europhys. Lett. 1996, 36, 259. (27) Lowell, S. Powder Technol. 1980, 25, 37. (28) Lowell, S.; Shields, J. E. J. Colloid Interface Sci. 1982, 80, 192. (29) Giesche, H. Mater. Res. Soc. Symp. Proc. 1982, 80, 192.
Langmuir, Vol. 20, No. 15, 2004 6483
namics and molecular simulations. From a molecular point of view, the nonwetting nature of the mercury can be treated by a repulsive interaction between the fluid and the solid. As shown by, for example, Page and Monson,30 the adsorption isotherm for such a system will occur at pressures higher than the saturated liquid pressure. Therefore the adsorption/desorption isotherm can be considered as the intrusion/extrusion curve of a nonwetting liquid into the porous material. Thus, in principle all the techniques applied to modeling adsorption/desorption of gases into porous materials can be applied to the intrusion/extrusion of nonwetting liquids and by extension to the modeling of mercury porosimetry experiments. Of course we have to be aware of the differences in behavior that might arise from differences in dynamics in the two situations (the entrapment of mercury on extrusion is an example of this). Dynamic aspects of mercury porosimetry will be addressed in the context of the present modeling approach in a subsequent paper. We describe the system with a lattice model which is simple enough to allow an extensive program of calculations but also sufficiently realistic to describe experiments at a microscopic level. Indeed, for fluids in mesoporous materials characterized by typical domain sizes much larger than molecular length scales, many of the microscopic details of the molecular interactions are expected to be less important to represent the qualitative behavior observed in experiments. Lattice models are very efficient computationally, and the implementation of mean field density functional theory (MF-DFT)8-10 and Grand Canonical Monte Carlo (GCMC) simulations9,10 is quite straightforward. In addition within the framework of this model the transition from a wetting fluid to a nonwetting fluid behavior is well-defined and controlled by the symmetry of the model. We choose to study the behavior of nonwetting liquid in a mesoporous Vycor glass which presents interesting mechanical properties and can be prepared in a relatively well-controlled manner in experiments. Our lattice model is the same as that used in recent studies of gas adsorption in disordered porous materials.8-10,13 The outline of our paper is as follows. In section II we present the main characteristics of our model together with the numerical methods (MF-DFT and GCMC) we have used in this work. Section IV is devoted to the results obtained, and in section V we give a summary of our results and conclusions from the work. II. Model and Methods A. Lattice Model. In our model the porous space is discretized on a face-centered cubic (fcc) lattice. For each lattice site i, the solid occupation variable ti is either 0 if the lattice site is occupied by a solid particle or 1 if the lattice site is available for the occupancy by a fluid particle. In the former case, the lattice site can be either occupied by a fluid particle and the fluid occupation variable ni ) 1 or “free” and then ni ) 0. The Hamiltonian of the lattice model H for a fixed configuration of the solid site can be written as
H ) -ω
∑ ninjtitj - Rω ∑ [niti(1 - tj) + 〈ij〉
〈ij〉
njtj(1 - ti)] - µ
∑i niti
(2.1)
where ω is the fluid-fluid interaction parameter, R is the (30) Page, K. S.; Monson, P. A. Phys. Rev. E 1996, 54, 6557.
6484
Langmuir, Vol. 20, No. 15, 2004
Porcheron et al.
B. Mean-Field Density Functional Theory. Density functional theories (DFTs) have been widely used for the modeling of gas adsorption in confined systems and have played an influential role in the understanding of confined fluids behavior.8-10,33-35 At a mean-field level of description (MF-DFT), this rather simple method is still able to give very good agreement with experimental features such as curves shapes, temperature dependence of the hysteresis, or existence of scanning curves.8 In mean field theory (MFT) the grand potential for our lattice model of the fluid plus solid system can be written as
βΩ[ti] )
∑i [Fi ln(Fi) + (1 - Fi) ln(1 - Fi) - βµFi] βω ∑ [FiFj + RFi(1 - tj) + RFj(1 - ti)] (2.3) 〈ij〉
where Fi is the average fluid density at site i in the system. The minimization of the grand potential then yields the mean-field equation for the local density Figure 2. Mesoporous Vycor glass model. The shaded points represent the solid lattice site and the light ones represent the sites avalaible for a fluid particle. p ) 0.3, a ) 30 Å, facecentered cubic lattice.
ratio of the wall-fluid interaction to the fluid-fluid interaction, and µ is the chemical potential.31 The double summations in eq 2 run over all the distinct nearest neighbors site pairs. We have used configurations of the solid site that represent the pore structure of Vycor glass. The Vycor glass model is built using the Gaussian random field method. The procedure has been fully described in the work of Woo and Monson10 and we only give a short summary here. Gaussian random field methods originate with the Cahn model of spinodal decomposition,32 where an auxiliary random field Φ(r) is used to generate the physical density of the solid via the level cut
Fs(r) ) Θ[Φ(r) - Φc]
(2.2)
where Θ(x) is the Heaviside step function. The lattice variable ti is obtained via the discretization ti ) 1 - Fs(ri). The field cutoff and the spatial correlations of the field can be fixed by the requirement that eq 2.2 correctly yields the statistical properties of the solid density distributions. The spatial fluctuations of the solid density are then designed to match known experimental properties of the material. For the Vycor glass prominent characteristics of the material are the porosity and the experimental structure factor I(q). The Vycor glass model is then characterized by the porosity defined as p ) 1 - Ns/Ntot, where Ntot is the total number of lattice sites and Ns the number of solid sites, and by a lattice constant a setting the scale of resolution of the method. With this method, realistic samples of the Vycor glass can be easily generated within the framework of a lattice model (see Figure 2). We also checked that the pore network displays a bicontinuity throughout the threedimensional space. Periodic boundary conditions are used in all directions of space. A bulk interface can also be added by the explicit modeling of bulk reservoirs along the z direction at the two extremities of the solid material. (31) We find it convenient to include the term involving the chemical potential since the lattice gas Hamiltonian is then isomorphic with that of an Ising type model of a magnet and it permits us to exhibit the symmetry of the model more easily. (32) Cahn, J. W. J. Chem. Phys. 1965, 42, 93.
Fi )
ni 1 + exp[-βµ - βω
[Fj + R(1 - tj)]] ∑ j∈i
(2.4)
which applies for each lattice site i. The equations are easily solved by an iteration scheme, and the convergence is generally rapid except in the regions where the density changes quickly with the chemical potential. At low temperatures these equations have multiple solutions corresponding to local minima of the grand potential. Such behavior is illustrated by the presence of scanning curves observed in MF-DFT calculations of gas adsorption in mesoporous Vycor glasses.8,13 III. Model Symmetry and Transformation of Intrusion/Extrusion Curves into Adsorption/ Desorption Isotherms Early work of Moscou and Lub36 showed that the pore size distribution of a γ-alumina material obtained by nitrogen desorption is similar to the distribution found by inverting the mercury intrusion curve on the same material. In subsequent work Murray et al.37 and Lowell and Shields38 have proposed equations for transforming mercury porosimetry curves into gas adsorption/desorption isotherms. In this section we explore the underlying basis for such transformations in the context of our molecular model. A. Lattice Model Symmetry. We can explore the relationship between gas adsorption and mercury porosimetry more clearly using the symmetry of our lattice model between adsorption/desorption of a wetting fluid and extrusion/intrusion of a nonwetting fluid. We begin by showing the symmetry in our model between a configuration of a wetting fluid in a mesoporous matrix and that of complementary configuration of a nonwetting fluid with respect to µ ) µ0 (the saturation chemical potential) and R ) 1/2. To see this more clearly, it is useful to rewrite the Hamiltonian for our lattice model using the (33) Balbuena, P. B.; Gubbins, K. E. Langmuir 1993, 9, 1801. (34) Woywood, D.; Schoen, M. Phys. Rev. E 2003, 67, 026122. (35) Bock, H.; Schoen, M. Phys. Rev. E 1999, 59, 4122. (36) Moscou, L.; Lub, S. Powder Technol. 1981, 29, 45. (37) Murray, K. L.; Seaton, N. A.; Day, M. A. Langmuir 1999, 15, 8155. (38) Lowell, S.; Shields, J. E. Powder Technol. 1981, 29, 225.
Modeling Mercury Porosimetry
Langmuir, Vol. 20, No. 15, 2004 6485
language of Ising models of magnets.10,13,39 Introducing a variable si which is +1 if site i is occupied by a fluid molecule and -1 if it is empty, then we can write
ω
H′ ) -
4
∑ sisjtitj - ∑i hisiti
(3.1)
〈ij〉
where hi is given by
(µ + ZRω)
hi )
2
ω + (1/2 - R) tj 2 j∈i
∑
(3.2)
and has the form of a magnetic field in an Ising model. The sum runs over the nearest neighbors of the site i and Z is the coordination number of the lattice. H′ differs from H only by terms which are independent of the configuration of fluid molecules. We now further define
R ) 1/2 + δR
(3.3)
µ ) µ0 + δµ
(3.4)
where µ0 is the chemical potential at bulk saturation (µ0 ) -Zω/2 for the lattice gas model with coordination number Z). Positive values of δR are associated with a fluid which wets the solid, and negative values of δR are associated with a nonwetting fluid. Positive δµ corresponds to the liquid state of the bulk fluid and negative δµ to the vapor state of the bulk fluid. It can readily be seen that H′ is unchanged after a transformation that changes the sign of δR, δµ, and all the si values while leaving the ti values fixed. Thus the Hamiltonian of a state of a fluid with positive δR and negative δµ (corresponding to gas adsorption/desorption of a wetting fluid) has the same value as that for a fluid with negative δR and positive δµ of the same magnitudes as in the wetting case and with all the fluid site occupancies reversed (corresponding to intrusion/extrusion of a nonwetting fluid). This establishes the symmetry between gas adsorption and porosimetry in the context of our model. B. Transformation of Intrusion/Extrusion Curves into Adsorption/Desorption Isotherms. We now show that from this microscopic symmetry, mercury intrusion/ extrusion curves can be transformed into equivalent adsorption/desorption curves and we relate the results to those obtained previously. We first relate the hydraulic pressure, Ph, of the mercury to the chemical potential as follows. Integration of the Gibbs-Duhem equation yields
µ l - µ0 )
∫PP VldP h
0
(3.5)
where µl is the chemical potential of the liquid phase, µ0 is the value at saturation, P0 is the saturation pressure, and Vl is the molar volume of the bulk liquid. If the liquid phase is considered to be incompressible and we take note of the fact that the hydraulic pressure in mercury experiments is very large compared with the saturation pressure, eq 3.5 reduces to
µl - µ0 = VlPh
(3.6)
We can rewrite this expression in terms of the activity as (39) Pitard, E.; Rosinberg, M. L.; Stell, G.; Tarjus, G. Phys. Rev. Lett. 1995, 74, 4361.
λl ) exp(PhVl/RT) λ0
(3.7)
where λl is the activity of the liquid phase and λ0 is its value at saturation. Lowell and Shields2,38 have presented an expression that relates the hydraulic pressure to the pressure of the mercury vapor that coexists with the liquid mercury inside the porous material (we can suppose that in a mercury porosimetry experiment the porous material is occupied by either intruded liquid mercury or mercury vapor). We will denote this latter pressure as Pv′. We can derive the Lowell and Shields expression from eq 3.7 as follows. If the mercury vapor inside the porous material has the properties of a bulk ideal gas, then we may calculate its pressure relative to the value at saturation from the activity as follows
Pv′ λv ) P0 λ0
(3.8)
Substituting eq 3.8 into eq 3.7 gives
Pv′ ) exp(PhVl/RT) P0
(3.9)
This is the equation given by Lowell and Shields.2,38 The equation transforms the hydraulic pressure scale for an intrusion experiment into a relative pressure scale for a gas adsorption experiment, but since we have a nonwetting fluid, the pressures are greater than P0. If we now assume that the system follows the same symmetry in activity with respect to saturation as our lattice model, then we can rewrite eq 3.7 as
λ v λ0 ) ) exp(-PhVl/RT) λ0 λl
(3.10)
This transforms the hydraulic pressure scale for an intrusion experiment into a relative activity scale for a gas adsorption experiment. If we further assume that the bulk gas in a gas adsorption experiment is ideal, then we can rewrite eq 3.10 as
Pv ) exp(-PhVl/RT) P0
(3.11)
and we now have an expression relating the hydraulic pressure scale for the porosimetry experiment to a relative pressure scale for a gas adsorption experiment. The expression of Murray et al.37 is similar to this but was derived by combining the Kelvin and Washburn equations and includes a correction to the Kelvin equation for the thickness of an adsorbed film on the pore surface. The remaining part of the transformation for both of these cases is to express the amount of mercury in the pores as an amount of fluid in the pores in a gas adsorption experiment. In the case of our lattice model, that is accomplished using the transformation
Fv ) 1 - Fl
(3.12)
which follows from the hole-particle symmetry of the lattice model. Note that all the expressions in this section were derived without any reference to the Kelvin or Washburn equations. Indeed the only assumptions we have used are the symmetry of the activity or chemical
6486
Langmuir, Vol. 20, No. 15, 2004
Figure 3. Adsorption/desorption curves of a wetting fluid (R ) 1.0, left part) and intrusion/extrusion curves of a nonwetting fluid (R ) 0.0, right part) in Vycor. MF-DFT, T/Tc ) 0.4, βµ0 ) -5. The curves are symmetric with regard to (µ ) µ0, F ) 1/2).
Porcheron et al.
Figure 4. Intrusion/extrusion curves of a nonwetting fluid (R ) 0.0, solid line) as a function of the hydraulic pressure. MFDFT, T/Tc ) 0.4, p ) 0.3.
potential about the saturation value (exact for the lattice model) and the incompressibility of the liquid phase. Thus the accuracy of eqs 3.9 and 3.11 does not require the validity of the Kelvin or Washburn equations. A discussion of the significance of these transformations is in order. The Lowell and Shields transformation gives us literally the adsorption/desorption isotherm for mercury vapor in the porous material. Since the mercury is nonwetting, the relevant values of pressure are P > P0. On the other hand, the approach of Murray et al.37 seeks to transform the intrusion/extrusion isotherm into a gas adsorption isotherm for a wetting fluid such as nitrogen (relevant values of pressure are now P < P0). In the case of our lattice model, that transformation can be accomplished exactly because of the symmetry of the model. In nature that symmetry is only roughly correct. IV. Results We first present results obtained by MF-DFT on a Vycor glass with a model porosity of p ) 0.3 and a lattice parameter a ) 30 Å. The number of fcc unit cells used in all directions of space of the sample is N ) 32 and results are obtained by averaging over 15 realizations of the geometry. Later we present some results from Monte Carlo simulations. A. Isotherms and Model Symmetry. Figure 3 shows the density versus chemical potential isotherms for wetting and nonwetting fluids with R ) 1 and R ) 0. The adsorption/desorption curves of the wetting fluid and the extrusion/intrusion curves of the nonwetting fluid are clearly symmetric with respect to µ ) µ0 and F ) 1/2. In both cases the system exhibits hysteresis between the adsorption/extrusion and desorption/intrusion. The origin of this hysteresis has been studied very recently for the wetting case.8,13 It is associated with the appearance of a very large number of local minima of the grand free energy and may or may not be associated with an underlying phase transition. Even if there is an underlying phase transition, very slow dynamics will make the phase coexistence inaccessible on the time scale of practical experiments.10 The symmetry of the lattice model dictates that the origin of hysteresis in the wetting and nonwetting cases be identical. This leads to the conclusion that in the first instance hysteresis in mercury porosimetry for disordered porous materials has the same origin as
Figure 5. Adsorption/desorption curves of a wetting fluid (R ) 1.0, solid line) and transformed intrusion/extrusion curves of a nonwetting fluid (R ) 0.0, b) in Vycor. MF-DFT, T/Tc ) 0.4, p ) 0.3.
hysteresis in gas adsorption/desorption experiments in the same materials. Figure 4 shows the density of the nonwetting fluid versus bulk (hydraulic) pressure. Values of the bulk pressure in mean field theory are obtained from the bulk equation of state. By using eqs 3.11 and 3.12 we can transform this into an adsorption/desorption isotherm of density versus relative pressure for the wetting fluid. This is shown in Figure 5 together with the results calculated directly, and the agreement is excellent, as it should be given that the only approximation in this case is eq 3.6. As a final illustration of the results from transformations, we show results in Figure 6 from plotting the intrusion/extrusion curve for the nonwetting case versus P′v/P0 using eq 3.9 and using the density transformation of eq 3.12. These results are to be compared with those in Figure 3 of the paper by Lowell and Shields.38 Lowell and Shields interpret these results as pore condensation/evaporation isotherms for the mercury. B. Relationship to Experiments. Mercury porosimetry experiments are performed at ambient temperature T ) 298 K, and the porosimeters are designed to reach high hydraulic pressures up to Ph = 4000 bar.3 The critical temperature of mercury is 1750 K, and the liquid-vapor saturation pressure at T ) 298 K is P0 ) 2.46 × 10-6 bar.
Modeling Mercury Porosimetry
Figure 6. Intrusion/extrusion curves of a nonwetting fluid (R ) 0.0, solid line) as a function of P′v/P0 (see eqs 3.9 and 3.12). MF-DFT, T/Tc ) 0.4, p ) 0.3.
Figure 7. Intrusion/extrusion curves as a function of Ph/P0 obtained on a model of Vycor glass for different wall-fluid interaction parameters, R. From left to right R ) 0.0, -0.5, -1.0. MF-DFT calculations, T/Tc ) 0.17, p ) 0.3, a ) 30 Å.
Thus in a mercury porosimetry experiment the reduced temperature is T/Tc ) 0.17 and relative pressures are typically in the range Ph/P0 ∈ [1-109].40 Let us consider results from our model at this reduced temperature. The critical temperature of the lattice-gas model in MFT is related to the lattice coordination number Z by kTc/ ) Z/4. We can therefore set the temperature in our calculations such that T/Tc ) 0.17 and obtain intrusion/extrusion curves as a function of the chemical potential µ. This is shown for three values of the wall-fluid repulsive interaction R in Figure 7. For small repulsive fields (R ) 0.0) the intrusion of the mercury within the whole Vycor occurs in a sharp step once the chemical potential has reached a given threshold. As the repulsive field is increased (i.e., R decreases), the intrusion of mercury in the Vycor is more continuous as it is shifted to higher µ values and spreads over a wider range of chemical potential. On the other hand, the extrusion of mercury (40) It is worth noting that to build a detailed off-lattice model of mercury porosimetry, an intermolecular potential capable of modeling the small ratio of triple point temperature to critical temperature for mercury would be necessary. In particular, the reduced temperature T/Tc ) 0.17 lies below the triple point temperature for the LennardJones 12-6 potential, so this potential would not be suitable. The lattice model does not present this problem since there is no solid-phase intervening at low temperatures.
Langmuir, Vol. 20, No. 15, 2004 6487
from the Vycor always displays a continuous shape regardless of the value of R. Once again, the extrusion range is larger for smaller values of R. All the curves present a hysteresis region which narrows as we make the wall-fluid interaction more repulsive. The pressure range where the intrusion/extrusion processes take place is roughly located between Ph/P0 ∈ [105-106]. This is somewhat lower than the highest relative pressures seen in mercury intrusion/extrusion observed in porosimetry experiments. This difference could reflect the difference between the volumetric properties of the lattice fluid model and those for mercury. The shapes of the calculated curves are also in good qualitative agreement with experimental results obtained on different kinds of porous materials. Indeed experiments performed on materials such as sandstone,41 dolomite,42 calcium-silicate glasses,43 or sol-gel silica spheres22 show the existence of a hysteresis and a steeper intrusion curvature as observed in our simulations. C. Scanning Curves. Scanning curves are an important feature of gas adsorption experiments performed on mesoporous solids in the hysteresis region.1 These curves are obtained by reversing the direction of the applied external pressure in the hysteresis region before the complete adsorption/desorption of the gas into/out of the solid material. Recent theoretical work has shown that these curves can be generated in MF-DFT calculations performed on lattice models of adsorption/desorption in Vycor8 and other systems.13 Scanning curves are observed during mercury porosimetry experiments as well. For instance in the work of Wardlaw and McKellar42 several cycles of mercury intrusion/extrusion/reintrusion were performed on different types of sedimentary rocks. The experiments showed that when the pressure is suddenly reversed, scanning curves following different paths from the complete intrusion/extrusion curves are obtained. On the intrusion branch these curves display an upward curvature whereas the curvature is downward when the pressure is suddenly increased from an extrusion process. We can model this phenomenon by reversing the direction of the evolution of the chemical potential within the boundaries of the hysteresis during an extrusion or an intrusion path (see Figure 8). We refer to the scanning curves as extrusion scanning when the chemical potential is decreased from a state on the intrusion curve and to intrusion scanning when µ is increased from an extrusion process. The curvatures of extrusion scanning curves are upward, and those of intrusion scanning curves are downward. The curvatures are more pronounced when the scanning curves are obtained from a deeper penetration in the hysteresis region. The shapes of scanning curves obtained from our simulations are in good agreement with the ones obtained during mercury porosimetry experiments. Despite the simplicity of our model, we are able to reproduce at least qualitatively the main features observed during a mercury porosimetry experiment in mesoporous materials. In fact, the main difference lies in the mercury retention obtained during the extrusion process in many experiments, and at very low pressure the density in the extrusion branch does not reach zero.3 This behavior is attributed to the entrapment of mercury within the network of the mesoporous solids. Alternative interpretations have been proposed such as in the work of Pirard et al.15 These authors suggested that the apparent mercury retention at very low pressures reflects mechan(41) Chatzis, I.; Dullien, F. A. L. Powder Technol. 1981, 29, 117. (42) Wardlaw, N. C.; McKellar, M. Powder Technol. 1981, 29, 127. (43) Saravanapavan, P.; Hench, L. L. J. Non-Cryst. Solids 2003, 318, 14.
6488
Langmuir, Vol. 20, No. 15, 2004
Porcheron et al.
Figure 8. Extrusion (left) and intrusion (right) scanning curves. Mean-field DFT calculations averaged over 10 realizations of the Vycor glass, R ) -0.5, T/Tc ) 0.17, p ) 0.3.
ical deformation and collapse of the pore network for the silica gel material they studied. Rigby et al.21 demonstrated, using light microscopy experiments on a different silica gel material, that pore collapse does not occur, and therefore mercury is really entrapped within the material. Moreover, recent neutron scattering experiments by Makri et al.44 confirmed that mercury porosimetry experiments performed on porous glasses lead to the entrapment of mercury within the pore network. Entrapment cannot be predicted by our MF-DFT calculations since the states obtained by solving the DFT represent at least a local minimum of the grand potential for the system and complete extrusion must occur in the limit of low pressure. Entrapment at low pressure must presumably be a dynamic phenomenon. In the next section, we further investigate the intrusion/extrusion processes and the issue of mercury retention within the Vycor glass by studying the equilibration processes in Monte Carlo simulations. D. Monte Carlo Simulations of the Intrusion/ Extrusion Processes. We have performed GCMC calculations on a Vycor sample with a geometry slightly different than the one used so far. The geometry is organized as N × N × nN along the x, y, and z directions respectively where N ) 32 is the number of fcc unit cells along each direction and n ) 5 is the number of cubic units in the z direction. Periodic boundary conditions are imposed on the x and y directions, whereas the Vycor solid in the z direction is in direct contact with two bulk reservoirs of size N at each extremities of the sample. The fcc cubic unit cell size used is a ) 15 Å. Again the curves are calculated in a sequential way, and at a given chemical potential the configuration is equilibrated for 5000 GCMC steps (each step involves an attempt at changing successively the occupancy of each fluid site in the lattice) and then 5000 further GCMC steps are performed for averages. This simulation can be view as a lattice equivalent of the geometry used in recent studies of the dynamics in the hysteresis region.11,12 In Figures 9 and 10 we show snapshots representative of the mercury intrusion and extrusion, respectively. From a completely empty sample (including the bulk reservoirs), we increase the chemical potential µ. The process first starts with a filling of the bulk reservoirs, immediately followed by the intrusion of mercury within the pores in direct contact with the bulk liquid. When the chemical potential further increases, the fluid spreads throughout the Vycor sample. Preferentially larger pores are filled first. However some of these large (44) Makri, P. K.; Stefanopoulos, K. L.; Mitropoulos, A. Ch.; Kanellopoulos, N. K.; Treimer, W. Physica B 2000, 276-278, 479.
Figure 9. Representative snapshots of the mercury intrusion within a model of Vycor glass in explicit equilibrium with bulk reservoirs. GCMC simulations, T/Tc ) 0.17, R ) -0.5, p ) 0.3, a ) 15 Å, n ) 5.
Figure 10. Representative snapshots of the mercury extrusion within a model of Vycor glass in explicit equilibrium with bulk reservoirs. GCMC simulations, T/Tc ) 0.17, R ) -0.5, p ) 0.3, a ) 15 Å, n ) 5.
pores can remain empty until higher chemical potentials are reached and pores with smaller diameters can even be filled first. After these large pores are filled, the fluid continues to intrude in the smallest pores of the material. Evidently, two processes are at work here. The first is the
Modeling Mercury Porosimetry
Langmuir, Vol. 20, No. 15, 2004 6489
graduate penetration of liquid front into the system from the external surfaces. The second is transfer of fluid within the material. In the GCMC simulations these processes are primarily modeled by insertion/deletion of molecules, and the time evolution of the system follows Glauber dynamics.10 In nature, these processes would occur by diffusive mass transfer although there is evidence from simulations of gas adsorption/desorption in the hysteresis regime that these mechanisms can be similar.12 In work in process, we are studying the intrusion/extrusion process using Kawasaki dynamics simulations of the model which gives more realistic treatment of diffusive mass transfer. Once the sample is completely filled, we progressively decrease the chemical potential to model the extrusion. As seen in ref 10 the extrusion process is rather different from intrusion. We see in particular that the liquid becomes progressively fragmented as the pressure is lowered and is removed from the system in a rather nonuniform manner. It is tempting to identify the appearance of pockets of liquid in regions of the system at low pressures on extrusion with the entrapment processes seen on extrusion in mercury porosimetry. However, in the GCMC simulations the fluid is completely extruded at low pressure. Our Kawasaki dynamics simulations should reveal more information about the extent to which we can model the entrapment process. As has been shown in recent modeling of gas adsorption/desorption in Vycor,10 Glauber dynamics does give a qualitatively accurate picture of the dynamics in the hysteresis region. As an illustration of this in the present context, we have calculated the relaxation time for density fluctuations in our GCMC simulations. The relaxation time τg is defined as10
τg )
∫0t
z
C(t) dt
(4.1)
where tz is the first zero of C(t) and
C(t) )
〈
〉
δF(0)δF(t) δF(0)δF(0)
(4.2)
The density fluctuation is defined as
δF(t) ) F(t) - 〈F〉
(4.3)
where F(t) is the instant total density at a time t (i.e., a given step of the GCMC simulations) and 〈F〉 is the average density. Practically, we proceed by performing long GCMC simulations for different chemical potentials. Then we calculate the corresponding autocorrelation function C(t) and calculate the relaxation time by performing an integration of the correlation function over different time origins. In Figure 11, we plot the relaxation time τg during the extrusion process. Clearly, the relaxation time sharply increases when µ reaches the hysteresis boundaries in the Vycor. This behavior can be seen as the “mirror” process of the relaxation time evolution observed in the work of Woo and Monson10 for the case of a wetting fluid adsorption in a Vycor glass. The increase in relaxation time on extrusion in the hysteresis regime is at least suggestive of the origins of mercury entrapment, especially when coupled with the visualizations in Figure 10. The dynamic nature of the entrapment is also suggested by an observation of Lowell and Shields2 during mercury porosimetry
Figure 11. Relaxation time τg (b) obtained during the extrusion process together with the intrusion/extrusion curves (dashed line). GCMC simulations, T/Tc ) 0.17, R ) -0.5, p ) 0.3.
experiments stating that at the completion of an intrusion-extrusion cycle, mercury will slowly continue to extrude sometimes for hours. Conclusions We have performed a set of mean-field density functional theory and grand canonical Monte Carlo simulations to model mercury porosimetry experiments on a molecular model of Vycor. The Vycor samples are generated with the Gaussian random field method which provides a realistic description of the pore network.8 We decomposed the system into a lattice model where the interaction parameter R between a fluid and a solid sites is set to repulsive values (R < 0.5) in order to model the nonwetting nature of the mercury with respect to many solids. The results obtained describe qualitatively the main features observed in mercury porosimetry experiments. For instance, the bulk pressures of intrusion/extrusion process are very high, in agreement with the experiment, and we are also able to qualitatively describe the origin of the scanning curves observed in some experiments. The intruding fluid tends to fill the larger pores at low pressures even though the Washburn equation is not completely verified since some small pore diameter can be filled before larger ones. Our MF-DFT calculations and Monte Carlo simulations do not explicitly describe the mercury entrapment observed in experiments (when the extrusion curve does not close the intrusion one at low pressure). However visualizations from GCMC simulations of the model system suggest that liquid fragments can remain within the Vycor system before the complete extrusion. The calculation of the relaxation time τg then shows that the time needed for these fragments to extrude from the Vycor may increase dramatically, and this could be responsible for the entrapment of mercury observed in experimental situations. As a work in progress, we are examining the dynamic aspects of these systems in more detail using Kawasaki dynamics simulations. A closer comparison with experiment will be possible once we have evaluated the importance of the relaxation dynamics in both extrusion and intrusion. Acknowledgment. This work was supported by the National Science Foundation (CTS-0220835). LA049939E
