pubs.acs.org/Langmuir © 2010 American Chemical Society
Monte Carlo Optimization Scheme to Determine the Physical Properties of Porous and Nonporous Solids L. F. Herrera, Chunyan Fan, D. D. Do,* and D. Nicholson School of Chemical Engineering, University of Queensland, St. Lucia, Queensland 4072, Australia Received May 19, 2010. Revised Manuscript Received July 26, 2010 A new method, based on a Monte Carlo scheme, is developed to determine physical properties of nonporous and porous solids. In the case of nonporous solids, we calculate the surface area. This surface area is found as the sum of areas of patches of different surface energy on the solid, which is assumed to take a patchwise topology (i.e., adsorption sites of the same energy are grouped together in one patch). As a result of this assumption, we derive not only the surface area, but also the accessible volume and the surface energy distribution. In the case of porous solids, the optimization method is used to derive the surface area and the pore size distribution simultaneously. The derivation of these physical properties is based on adsorption data from a volumetric apparatus. We test this novel idea with the inversion problem of deriving surface areas of patches of different energies for a number of nonporous solids. The method is also tested with the derivation of the pore size distribution of some porous solid models. The results are very encouraging and demonstrate the great potential of this method as an alternative to the usual deterministic optimization algorithms which are known to be sensitive to the choice of the initial guess of the parameters. Since the geometrical parameters are physical quantities (i.e., only positive values are accepted), we also propose a scheme to enforce the positivity constraint of the solution.
1. Introduction The development of computation simulation methods such as density functional theory (DFT) and grand canonical Monte Carlo (GCMC) simulation gave the foundations for new and advanced methods for the characterization of porous solids. These methods are based on the optimization of the matching between an experimental isotherm at a given temperature and the theoretical isotherm calculated from a set of local isotherms at the same temperature. The local isotherms are constructed from a knowledge of the characteristics of the adsorbent solid such as the topology of the pores (slit and cylindrical pores) and the surface structure of the solid (lamellar or slab structure). This a priori information is used to choose an appropriate potential model of interaction between a fluid particle and the solid. Once this is done, a set of local isotherms can be constructed by using tools such as DFT or GCMC. This set could be one with different surface energies or different pore sizes, or in general, it could be one with different unit cells. In this case, each unit cell would represent a characteristic of a solid, for example, surface energy and pore size as just mentioned. This characteristic can be general, and it could represent a particular connectivity of the solid, for example. Once the kernel has been defined, calculated and stored, we can use it in an appropriate mass balance equation to obtain the amount of adsorbate that is dosed into the adsorption cell. This theoretical solution is then matched against the experimental data to derive the parameters of the solids. This matching is usually carried out with a deterministic optimization method, for example, simplex methods or gradient methods. The common denominator in all of these techniques is an initial guess of the solution, and it is well-known that the consequence of the optimization depends critically on this choice of the initial guess, as an incorrect choice would yield incorrect solution or the method always diverges and no solution can be found. *Author to whom all correspondence should be addressed. Fax: þ61-73365-4154; E-mail:
[email protected].
15278 DOI: 10.1021/la102017t
The objective of this paper is to develop a simple means, based on a Monte Carlo scheme, to offer a solution to the characterization problem. The method is simple and robust. We will test its potential with the inversion problem for the characterization of surface area of a nonporous solid, which is composed of patches of various energies. This method is also tested to derive the pore size distribution of some ideal porous solids. It can be used with isotherms at any temperatures, and in this paper, we apply this method to adsorption data at 87.3 and 180 K. The method that we develop here is not only restricted to characterization problems based on adsorption data, but also applicable to any problems that involve deriving parameters from a set of experimental measurements, provided that the required parameters appear in a linear form in the fitting equation (usually mass balance equation).
2. Theory We test the new method with the characterization of structural parameters from adsorption measurement using a volumetric system, which is a common method to obtain reliable equilibrium experimental data. In this method, a known amount of gas is dosed into the adsorption cell, sufficient time is allowed for the system to reach equilibrium and the equilibrium pressure is then recorded. This process is repeated with different amounts dosed into the adsorption cell, and a relationship is obtained between the equilibrium pressure and the dosing amount. Then, the structural parameters of the solid are derived by matching the total amount dosed in the adsorption cell to appropriate theoretical equations. This method is preferable to the use of the excess amount, which requires knowledge of the void volume of the system.1,2 (1) Birkett, G.; Do, D. D. New method to determine PSD using supercritical adsorption: applied to methane adsorption in activated carbon. Langmuir 2006, 22, 7622–7630. (2) Do, D. D., Do, H.D.; Birkett, G. A New Methodology In The Use Of Super-critical Adsorption Data To Determine The Micropore Size Distribution, in Adsorption, Progress In Fundamental and Application Research, Zhou, L., Ed.; World Scientific: Tianjin, 2006; pp 154-173.
Published on Web 09/02/2010
Langmuir 2010, 26(19), 15278–15288
Herrera et al.
Article
Figure 1. Schematic diagram of the adsorption cell with a surface having regions of different surface energy (we show three regions in this figure). The accessible volume is bounded by the dashed line. The particle A residing on the boundary has zero solid-fluid potential energy. The particle B is inside the accessible volume and therefore has a negative solid-fluid potential energy.
2.1. New Method. The amount dosed into the adsorption cell, N (mol), is equal to the sum of the amount in the gas phase and an excess due to adsorption. We assume that the surface is composed of M patches of different characteristic surface energy j*sf (in contrast to the BET theory,3 which assumes that the entire surface has the same energy). We shall adopt a patchwise topology, in which there is no interaction between the adsorbed molecules across the patches.4 With this assumption, we can write the mass balance of the adsorption cell as N ¼
M X
Sj Γj ðPÞ þ Vacc Ff
ð1Þ
j¼1
The linear sum in the RHS of eq 1 results from the assumption of patchwise topology. In this equation, P is the equilibrium pressure corresponding to the total amount dosed (N), Vacc is the accessible volume of the adsorption cell, Sj is the surface area of the patch j, and Γj is the surface excess density (mol/m2) of that patch. Figure 1 shows the schematic diagram of the adsorption cell and the accessible volume (shown as the region bounded by the dashed line). The accessible volume is defined as one in which a fluid particle has a nonpositive solid-fluid potential energy and the accessible surface area is the area of a surface on which the solid-fluid energy is zero. More details about the accessible volume and surface area can be found in earlier publications.5,6 For a given amount dosed (N), we have a corresponding equilibrium pressure P. At this pressure, we can calculate the gas density (F) from some appropriate equation of state and calculate the (3) Brunauer, S.; Emmett, P. H.; Edward, T. Adsorption of gases in multimolecular layers. J. Am. Chem. Soc. 1938, 60, 309–319. (4) Ross, S.; Olivier, J. P. On physical adsorption; Interscience Publisher, Inc.: New York, 1964 (5) Do, D. D.; Do, H. D. Appropriate volumes for adsorption isotherm studies: The absolute void volume, accessible pore volume and enclosing particle volume. J. Colloid Interface Sci. 2007, 316(2), 317–330. (6) Do, D. D., et al. A Monte Carlo integration method to determine accessible volume accessible surface area and its fractal dimension. J. Colloid Interface Sci. 2010, 348(2), 529-536.
Langmuir 2010, 26(19), 15278–15288
surface excess density from the set of local isotherms (the kernel). Knowing these variables, eq 1 is a linear equation in terms of surface areas of all patches and the accessible volume of the adsorption cell. Before we discuss the new optimization technique based on a Monte Carlo scheme, we discuss briefly the construction of a kernel. Each local isotherm is an isotherm for a patch with a particular surface energy (Ej) between the lowest (Emin) and the highest (Emax) surface energies used in the kernel. To obtain this isotherm, we have to make some further assumptions about the structure of the surface. For example, in the case of a solid with covalent bonding among the solid atoms such as alumina and silica gel, we model it as a solid slab of constant volumetric density. If the solid is graphite or has a lamellar structure, which is composed of parallel layers, we model it as a lamellar structure. For a surface having a lamellar structure, the interaction between a fluid particle and a layer of constant surface density is given by the Steele 10-4-3 potential equation7 " # 4 σ 4sf 2 σ sf 10 σsf jSf ¼ jSf 5 z z 3Δð0:61Δ þ zÞ3
jSf ¼ 2πðFS σ2sf Þεsf
ð2Þ
where jsf * is the characteristic energy, Fs is the surface density (number/m2). Δ is the spacing between layers (0.3354 nm) and z is the vertical distance between the fluid particle and the surface. The cross molecular parameters are usually calculated from the Lorentz-Berthelot rule: σsf = (σss þ σff)/2 and εsf = (εssεff)1/2, where the subscripts ff and ss are the parameters for the fluid and solid particles, respectively. After a model of the solid has been selected, grand canonical Monte Carlo (GCMC) simulation is used to calculate the isotherms for surface having specific surface energies. The set of these local isotherms at these specific surface energies is called the kernel, hereafter. Given the experimental data in the form of N versus P, we apply an optimization method to eq 1 to derive the structural parameters: surface areas of all patches of different energies and the accessible void volume of the system. Let K be the number of data points (N, P), that must be greater than the number of parameters (M þ 1) in eq 1 (M surface areas and one accessible void volume); otherwise, the problem is underdefined. The kth data point (Nk, Pk), must satisfy the mass balance equation eq 1, that is Nk ¼ m
M X j¼1
Sj Γj ðPk Þ þ Vacc Ff , k
ð3Þ
for k = 1, 2, ..., K. Equation 3 is linear with respect to the required parameters, M specific surface areas (Sj, j = 1, 2, ..., M), and the accessible volume. The above set of K equations can be cast into a vector-matrix form as follows: N ¼ Ax
ð4Þ
where N is (K 1) vector, A is a K (M þ 1) matrix containing the surface excess density and the density, and x is a (M þ 1) 1 vector of the required geometrical parameters (specific surface areas (Sj, j = 1, 2, ..., M) and accessible volume Vacc). Our task is to find a solution x from this set of K equations, and this is done by minimizing the residual of the sum of squares of relative errors (7) Steele, W. A. Physical interaction of gases with crystalline solids. Surf. Sci. 1973, 36(1), 317–352.
DOI: 10.1021/la102017t
15279
Article
Herrera et al.
between the theoretical value and the corresponding experimental value, i.e. " #2 K NðPk Þ - Nkexp 1X ð5Þ Residual ¼ φ ¼ Kk ¼ 1 Nkexp where Nexp k is the experimental number of moles introduced into the adsorption cell at the equilibrium pressure Pk. An ideal solution to eq 5 is zero, but this is not possible because the surface model that we use to obtain the surface excess density would never describe physical reality perfectly. Therefore, the best one could hope for is that (i) the surface modeled in the kernel is close to the physical surface and (ii) the residual is therefore minimized. Thus, our task here is to minimize the residual of eq 5 for a given choice of the kernel. We now present the details of a new Monte Carlo simulation scheme to achieve this goal. The key advantage of our scheme is that the method is robust and, most importantly, we do not need a guess for the solution to initiate the method. The following algorithm lists the steps taken to minimize the residual (eq 5) for a given kernel of local isotherms. Algorithm. The Monte Carlo optimization scheme is as follows: If J is the number of patches, then we are searching for J values of surface area for these patches and one accessible volume. The total number of parameters to be found is, therefore, J þ 1. 1. First, we choose a value of surface energy at random between the minimum and maximum surface energies used in the kernel, Emin and Emax, respectively. Then, the local isotherm corresponding to this energy is obtained from the kernel. Let this chosen surface energy be E(1). Because we start with only one patch, J=1, this means that the surface is initially assumed to contain only one patch of constant surface energy. Since E(1) is selected at random, it is a continuous value in the domain [Emin, Emax]. Because our task is to find distinct surface energy patches, we give each patch a range of ΔE. This is preferable to dividing the domain [Emin, Emax] into intervals of ΔE a priori if we are searching for the minimum number of patches required to characterize the surface. 2. We select J þ 1 equations at random out of the set of equations in eq 4. We denote this set as Λ, and the remaining equations of eq 4 is denoted by Ω. We write the subset Λ as follows: N 0 ¼ A0 x
ð6Þ
where N0 has a dimension of (J þ 1, 1) and the matrix A0 has a dimension of (J þ 1) (J þ 1). The vector x contains surface areas of J patches and one accessible volume. This set of linear equation of eq 6 is well-defined, i.e., J þ 1 linear equations in terms of J þ 1 unknowns, and therefore, we can readily find the solution by the standard inversion 0 -1
x ¼ ½A
N
0
ð7Þ
We now rewrite the residual of eq 5 by splitting the K experimental data points into (J þ 1) data points, that we just selected randomly to solve for the J þ 1 parameters, and the remaining data points as follows: 8 " #2 " #2 9 þ1 K X NðPk Þ - Nkexp NðPk Þ - Nkexp = 1
