Langmuir 1999, 15, 8235-8245
8235
Calculation of Pore Size Distributions of Activated Carbons from Adsorption Isotherms G. M. Davies,† N. A. Seaton,*,‡ and V. S. Vassiliadis† Department of Chemical Engineering, University of Cambridge, Pembroke Street, Cambridge, CB2 3RA, United Kingdom, and School of Chemical Engineering, University of Edinburgh, King’s Buildings, Mayfield Road, Edinburgh, EH9 3JL, United Kingdom Received March 5, 1999. In Final Form: July 26, 1999 The ability to determine whether a pore size distribution can be fitted to a given set of experimental data is relevant to the development of pore network models of the internal structure of activated carbons and other amorphous adsorbents. A closely related problem that arises when pore size distributions can be fitted to a given set of data is to determine suitable methods to calculate reliable, stable, and representative pore size distributions. This paper addresses both of these issues by presenting a procedure to determine whether a pore size distribution based on a specified model pore geometry can be fitted to a set of experimental data and then using this procedure, in conjunction with regularization techniques that stabilize the calculations, to determine statistically significant pore size distributions. This approach has been demonstrated using the adsorption of methane at 308 K onto BPL-activated carbon as a case study.
Introduction A pore size distribution (PSD) is frequently used to describe the internal structure of activated carbons and other amorphous microporous adsorbents. As such it assumes that the complex internal structure of adsorbents can be described in terms of an equivalent collection of regularly shaped model pores. The equivalence between the real and model structures is based on their adsorption behavior and the PSD is therefore calculated from the adsorption integral equation
N(Pi) )
∫0∞F(w, Pi) f(w) dw
i ) 1...n
(1)
where N(Pi) is the experimentally determined adsorption at pressure Pi, F(w, Pi) is the adsorption density in a model pore of characteristic dimension w, f(w) is the pore size distribution, and n is the total number of adsorption measurements used in the analysis. A number of aspects need to be taken into account when calculating a PSD from the adsorption integral equation. For example, a suitable shape for the model pores needs to be specified as well as an appropriate method to calculate the adsorption within these pores, suitable adsorption data needs to be identified and the overall validity of a PSD (as a suitable method to characterize the internal structure of microporous materials) needs to be addressed.1-9 This paper addresses two issues that arise once the shape of the model pores, the method to calculate the adsorption within these pores and the experimental data to be included in the analysis have been specified. The first issue is how to determine whether at least one pore † ‡
University of Cambridge. University of Edinburgh.
(1) Davies, G. M.; Seaton, N. A. Langmuir 1999, 15, 6263. (2) Bojan, M. J.; Steele, W. A. Carbon 1998, 36 (10), 1417. (3) Yin, Y. F.; McEnaney, B.; Mays, T. J. Carbon 1998, 36 (10), 1425. (4) Turner, A. R.; Quirke, N. Carbon 1998, 36 (10), 1439. (5) Olivier, J. P. Carbon 1998, 36 (10), 1469. (6) Davies, G. M.; Seaton, N. A. Carbon 1998, 36 (10), 1473. (7) McEnaney, B.; Mays, T. J.; Chen, X. Fuel 1998, 77 (6), 557. (8) Chen, X. S.; McEnaney, B.; Mays, T. J.; Alcaniz-Monge, J.; CazorlaAmoros, D.; Linares-Solano, A. Carbon 1997, 35 (9), 1251. (9) Quirke, N.; Tennison, S. R. Carbon 1996, 34, 1281.
size distribution based on the chosen model pores exists that can be fitted to the data. This is important because it can be used to identify pore network models of the internal structure that are not sufficiently realistic to account for the adsorption data being considered. The first section of this paper therefore presents a method that can be used to determine whether any pore size distributions exist that can be fitted to a given set of data. As we demonstrate, when at least one PSD can be fitted to a set of experimental data, an infinity of PSDs can be fitted. The second issue addressed in this paper is therefore how to choose a suitably representative PSD when an infinity of PSDs can be fitted to the data. As a case study we have calculated a PSD for BPLactivated carbon based on slit-shaped model pores and the adsorption of methane at 308 K (data from Gusev et al.10). In this study, the adsorption in individual pores with a range of widths was calculated using grand canonical Monte Carlo simulation; details of the simulation have been published elsewhere.1 Determining Whether a PSD Can Be Fitted to Experimental Adsorption Data Two aspects significantly complicate determining whether any PSD can be fitted to a given set of experimental data. The first is that we have no a priori knowledge of the forms of the PSD, if any, that are likely to fit the data. The second is that we have only a finite amount of data at our disposal. Frequently, analytical functions are proposed for the PSD. These functional forms are usually motivated by physical arguments and the fact that the parameters of these functions can be determined in a statistically significant manner (i.e., there are more data points than parameters). An example of a physical argument is that a real material is likely to exhibit a relatively smooth variation in pore sizes and that most of the pores will be centered about a few dominant pore sizes. Arguments of this type usually lead to n-modal gamma or log-normal functional forms of the PSD being adopted (see for example, (10) Gusev, V. Y., O’Brien, J. A.; Seaton, N. A. Langmuir 1997, 13 (10), 2815.
10.1021/la9902643 CCC: $18.00 © 1999 American Chemical Society Published on Web 09/18/1999
8236
Langmuir, Vol. 15, No. 23, 1999
refs 9 and 11-15). There is however a danger in adopting an analytical form for the PSD. The danger is that we are no longer determining whether any PSD can be fitted to the data but rather whether a PSD conforming to the proposed functional form exists. Therefore, if the functional form cannot be fitted to the data to a desired accuracy, we cannot exclude the possibility that other forms might have been successful. Increasing the flexibility of the adopted functional form (i.e., increasing the complexity or the number of adjustable parameters) would limit the possibility of not being able to calculate a PSD that fits the data due to it being constrained to an inappropriate functional form. Unfortunately, this approach introduces new difficulties. The adsorption integral equation cannot, in general, be inverted analytically, and one therefore needs to resort to numerical inversion. Numerical inversion relies on a discrete representation of the adsorption integral equation and solves the equation exactly only at those discrete points. The behavior of the function away from these points is not accounted for in the inversion procedure, which leads to the danger of highly oscillatory functions being calculated. (A similar problem is encountered when trying to fit a piecewise cubic spline to an arbitrary set of data points.) A recent study that increases the flexibility of the functional form adopted by representing the PSD as a piecewise cubic spline is that of Jagiello.16 In contrast to the representation of a PSD as an analytical function, a discrete representation is more flexible and does not attempt to quantify the PSD away from the quadrature points. Recent PSD analyses that have been based on discrete distributions include those of McEnaney et al.,7 Gusev et al.,10 Gusev and O’Brien,17 Samios et al.,18 Sosin et al.,19 and Sosin and Quinn.20 In this paper we extend the development of such an approach and show how it can be used to establish whether any PSD exists that can be fitted to the data and how it can be used to calculate a PSD representative of the real solid. (Of course, the true PSD in an amorphous (or only locally crystalline) solid such as an activated carbon is a continuous function. To be realistic, a discrete PSD would have to have a sufficient density of points to give a good approximation to a continuous function.) An initial problem that is encountered when discretizing the adsorption integral equation (which would also arise if a functional form of the PSD were adopted) is that the integration extends, in principle, to an infinitely large pore size. We suggest that the integration to infinity can be successfully incorporated into the analysis by splitting the pore size domain into two regions. The first region is the main “region of interest” and contains one or more quadrature intervals. The second region extends from the upper limit of the region of interest to an arbitrarily large pore size and contains only one quadrature interval. This single quadrature interval in the second region is merely used to account for all the adsorption that occurs in the (11) Wang, K.; Do, D. D. Langmuir 1997, 13 (23) 6226. (12) Do, D. D.; Do, H. D. Stud. Surf. Sci. Catal. 1994, 87, 641. (13) Aukett, P. N.; Quirke, N.; Riddiford, S.; Tennison, S. R. Carbon 1992, 30, 913. (14) Seaton, N. A.; Walton, J. P. R. B.; Quirke, N. Carbon 1989, 27, 853. (15) McEnaney, B.; Mays, T. J.; Causton, P. D. Langmuir 1987, 3, 695. (16) Jagiello, J. Langmuir 1994, 10, 2778. (17) Gusev, V. Y.; O’Brien, J. A. Langmuir 1997, 13, 2822. (18) Samios, S., Stubos, A. K.; Kanellopoulos, N. K.; Cracknell, R. F.; Papadopoulos, G. K.; Nicholson, D. Langmuir 1997, 13, 2795. (19) Sosin, K. A.; Quinn, D. F.; MacDonald, J. A. F. Carbon 1996, 34 (11), 1335. (20) Sosin, K. A.; Quinn, D. F. J. Porous Mater. 1995, 1, 111.
Davies et al.
pores larger than those in the region of interest. Adopting evenly spaced quadrature points in the region of interest we obtain
δwj )
wu - wl m-1
j ) 1...m - 1
(2)
where δwi is the quadrature interval, wu and wl are the upper and lower limits of the region of interest, and m is the total number of quadrature intervals in the PSD. The single quadrature interval of the second region is
δwm ) wmax - wu
(3)
where wmax is the maximum limit used in the analysis. The pore sizes at which the PSD is calculated, w/j , are defined to be the average pore size in each quadrature interval:
w/j ) wj +
δwj 2
(4)
The motivation for splitting the domain of the pore size distribution into two intervals is based in part on the “window of reliability” introduced by Gusev et al.10 The window of reliability identifies a pore size above which it is difficult to reliably distinguish one pore from another in the PSD analysis. This happens when the pores are so large that the adsorbed layers on opposing pore walls are so far apart that they do not interact with each other; as the upper limit of the window is approached, the isotherms begin to look increasingly similar to each other and to the isotherm for adsorption on a single surface. The upper limit of the window of reliability is a function of the adsorbate used and the temperature and the maximum pressure at which adsorption was measured. For the adsorption of methane at 308 K and up to a maximum pressure of approximately 30 bar, this pore size occurs between 15 and 20 Å. The window of reliability impacts on the region of interest that is chosen. The region of interest must extend to at least the upper limit of the window of reliability. In the case study presented in this paper this means that the region of interest needs to extend to at least 15 Å. If this were not the case, the single quadrature point extending to an arbitrarily large pore size would encompass a range of adsorption behavior (and isotherm shapes), the neglecting of which would distort the PSD in the region of interest. Although the region of interest in the analysis can extend beyond this limit, the exact locations of any peaks beyond the window of reliability are less well defined than those within it. We return to this point later. Having defined a suitable quadrature, the adsorption integral equation, eq 1, can be approximated as
N(Pi) ≈
F(w/j ,Pi)f(w/j )δwj ∑ j)1
(5)
Equation 5 can be written compactly in matrix notation as
N ) AWf where
(6)
Calculation of Pore Size Distributions
Langmuir, Vol. 15, No. 23, 1999 8237
N ) (N(Pi))i)1...n
(7)
A ) (F(w/j ,Pi))i)1...n,j)1...m
(8)
W ) diag(δwj)j)1...m
(9)
f ) (f(w/j ))j)1...m
(10)
Equation 6 in turn can be written as
N ) Af*
(11)
f* ) Wf
(12)
f ) W-1f*
(13)
where
and
It is useful to classify the calculation of PSDs according to the number of data points and quadrature intervals used in the analysis: (1) n > m, there are more data points than quadrature intervals; (2) n ) m, there are as many data points as there are quadrature intervals; and (3) n < m, there are more quadrature intervals than data points. The first case is unlikely to occur in practice because of the expense involved in acquiring experimental data. We therefore only investigate the second and third cases. The approach that we have adopted is to investigate an extremely straightforward example that involves initially fitting a PSD to only two data points. We use this example to highlight the important principles in determining PSDs and then generalize the results to include n data points. The two initial data points that we have used involve the adsorption of methane at 308 K onto BPL carbon at pressures of 1.14 and 20.33 bar.10 The experimental adsorption vector corresponding to these conditions is
N)
[ ] 0.80 3.79
mmol/g
(14)
For case 2, we can fit a discrete PSD that contains two pore sizes. Since we have already proposed splitting the pore size domain into two regions, defining the quadrature of the PSD merely involves identifying the region of interest and a suitable upper limit to use in the analysis. Because we are primarily interested in characterizing the microporous region, we have chosen the lower and upper limits of the region of interest to be the smallest pore in which adsorption takes place, which for methane is approximately 6.1 Å, and the pore width that defines the limit of the microporous region, which by IUPAC convention is 20 Å. (Note that 20 Å extends beyond the window of reliability.) A suitable upper limit to use in the analysis is a model pore in which the excess adsorption (in terms of the density of adsorbate in the pore) becomes negligible, relative to the adsorption in smaller pores. In our analysis, the adsorption of methane is already negligibly small for pores larger than 100 Å, and we therefore used this pore size as the upper limit in the analysis. The excess adsorption in slit-shaped model pores at 1.14 and 20.33 bar, calculated using grand canonical Monte Carlo molecular simulations,1 is shown in Figure 1. Using these results, the average excess adsorption between 6.1 and 20.0 Å for the low and high pressures form the A11 and A12 entries of the adsorption matrix, respectively. Similarly, the average excess adsorption between 20.0 and 100.0 Å form the A12 and A22 entries. (We have
Figure 1. Excess adsorption of methane at 308 K in slit-shaped model pores calculated using grand canonical Monte Carlo molecular simulations.
calculated the average excess adsorption using a mean value theorem approach. i.e., A11δwi, yields the exact area under the profile for the low-pressure adsorption between 6.1 and 20 Å. This approach minimizes the error introduced by discretizing the PSD.) Therefore, in this instance:
A)
[
1.65 0.12 7.37 1.53
]
mmol/cm3
(15)
From linear algebra we know that as long as the columns (or equivalently the rows) of A are linearly independent then a unique solution will exist to eq 11 which in turn will yield a unique f. However, this solution will be physically meaningful only if the two components of f* are nonnegative (since this will result in f having positive components). A “feasible” PSD can therefore be defined as one for which
f: ) f(w/j ) g 0
j ) 1...m
(16)
Solving eqs 11 and 13 with A defined by eq 15 yields the feasible solution
f)
[
0.034 0.003
]
cm3/gÅ
(17)
Our understanding of when a feasible PSD will exist can be improved by examining the column space of A. In Figure 2 we have plotted the two column vectors of A as well as N in R2. Plotted in this way the horizontal axis is the adsorption at the low pressure (1.14 bar) and the vertical axis is the adsorption at the high pressure (20.33 bar). Geometrically it is easy to see when we will be able to calculate feasible PSDs: only the points that lie between the positive projections of the two column vectors can be expressed as positive multiples of these vectors. We now address the case of when the experimental data point does not lie between the projections of the column vectors and it is therefore not possible to calculate a feasible PSD based on the given adsorption matrix. As an example, we will consider what would have happened if the adsorption at the higher pressure of 20.33 bar had been 3.50 mmol/g instead of 3.79 mmol/g. This leads to the following adsorption vector:
N)
[ ] 0.80 3.50
mmol/g
(18)
8238
Langmuir, Vol. 15, No. 23, 1999
Davies et al.
Figure 2. The two columns of the adsorption matrix (eq 15) and the experimental data points plotted in the adsorption space of 1.14 and 20.33 bar.
Figure 3. The columns of the adsorption matrix (eq 15 and 19) and the experimental data points plotted in the adsorption space of 1.14 and 20.33 bar.
In this case we would not have been able to calculate a feasible 2-pore PSD. In cases such as these, we need to investigate whether there exists a more flexible PSD, which contains more pore sizes, that could be fitted to the data. A slightly more flexible PSD is one that contains three pore sizes, in which case there are two quadrature intervals within the region of interest. Now the A matrix is not square since there are two data points and three pore sizes. This situation corresponds to case 3 identified earlier. As such, we immediately recognize that there will be an infinity of purely algebraic solutions (since there are more free parameters than equations). However, by including more column vectors we are able to extend the subset within R2 that can be expressed as positive scalar multiples of these vectors. If this subset can be extended to include the experimental data point, then at least one feasible PSD exists. The components of the new A matrix corresponding to the 3-pore PSD are as follows:
A)
[
2.82 0.48 0.12 9.01 5.74 1.53
]
mmol/cm3
(19)
Figure 3 compares the original two-column vectors (and the original data points of Gusev et al.10) and the three new column vectors (and the slightly perturbed data point at 20.33 bar). Note that the last column vector in the new A matrix corresponds to the last column vector in the
original A matrix. Also note that the second and third column vectors in the new A matrix are almost linearly dependent (i.e., they are very nearly collinear). The most important point in Figure 3 however is that the slightly perturbed data point, which does not lie between the two original column vectors (i.e., the case 2 scenario), does lie between the three new column vectors (i.e., the case 3 scenario). This means that at least one feasible 3-pore PSD exists. Unfortunately there is not only one feasible 3-pore PSD; there is an infinite number of such distributions. This fact can be used to prove that either infinity or no PSDs can be fitted to a given set of experimental points. This is because the only time a “unique” PSD can be calculated is when there are as many quadrature intervals as there are data points, that is, for a case 2 scenario. However, all case 2 scenarios can be converted to case 3 scenarios by merely including an additional quadrature interval into the analysis. Later we consider the implications of these infinite sets of solutions on the characterization. Having established that including more pore sizes into the analysis extends the subset within R2, we now investigate whether there is a number of pore sizes beyond which the subset in R2 is not significantly increased. Figure 4 shows the relationship between the subset in R2 and the number of pore sizes used in the analysis. In this figure we have included a trace running along the tops of the column vectors so that it is easier to identify the order in which they occur. This figure clearly shows that the subset within R2 is not increased significantly beyond 50 pore sizes. Such a subset can be termed the “limiting envelope” of reachable feasible solutions. It is also interesting to note in this figure that the adsorption in very small pores is equivalent to the adsorption in some larger pores; this is the reason the column vectors cross over each other. The preceding analysis can be used effectively to determine whether at least one PSD can be fitted to two experimental data points. We now consider how to determine whether at least one PSD can be fitted to n data points. In this case, the column space of the adsorption matrix lies in Rn. As before, feasible PSDs will only exist for sets of experimental data points that belong to the subset within Rn that is bounded by the positive projections of the column vectors of the adsorption matrix. In n dimensions we require a nongraphical method to determine whether the experimental data point lies within this subset. Although it may seem logical first to check whether a feasible PSD consisting of n pore sizes can be fitted to n data points (i.e., a case 2 scenario) using linear algebra, the nature of the problem often precludes this approach. Very often it is impossible to numerically calculate a reliable PSD because the adsorption matrix has become too ill-conditioned. This means that the column vectors of the adsorption matrix, A, are nearly linearly dependent, too much so to calculate a reliable solution. (In two dimensions, the near linear dependency of the column vectors for the large pore sizes can be seen in Figure 4. Notice however that the linear dependency is not quite as bad in Rn as it may seem in these figures. This is because there are more dimensions in Rn that allow the vectors to be “further away” from each other.) As an example of when this approach fails, we were unable to determine a PSD based on 19 quadrature intervals and the complete methane isotherm of 19 data points10 using a NAG library21 (21) NAG Fortran Library Manual - Mark 16, NAG Ltd, Wilkinson House, Jordan Hill Road, Oxford, OX2 8DR, United Kingdom, 1993.
Calculation of Pore Size Distributions
Langmuir, Vol. 15, No. 23, 1999 8239
Figure 4. Plots of the column vectors of the adsorption matrix in the adsorption space of 1.14 and 20.33 bar.
linear solver, routine F04ATF, due to the ill-conditioning of the adsorption matrix. (Strictly, the solver failed due to the ill-conditioning of the AW matrix.) It is possible however, for each of the three cases, to determine whether any feasible PSD exists using optimization routines. This is achieved by determining whether the following residual, subject to eq 16, can be reduced to below a specified tolerance: n
R)
m
[Ni(Pi) - ∑F(w/j ,Pi)δwj f(w/j )]2 ∑ i)1 j)1
(20)
The residual defined by eq 20 can be written more compactly as
R ) (N - AWf)T(N - AWf)
(21)
If, for a given set of data, a minimization routine is unable to reduce the residual to below a specified tolerance for a small quadrature interval, we can then conclude that no PSD exists that can be fitted to the data. It is worth noting that it is not possible for these routines to become trapped in a local minimum (and hence terminate before finding the global minimum) due to the fact that the problem is convex. More specifically, R is a positive semidefinite quadratic function in f. This implies that only one global minimum or several degenerate minima (all of equal value) exists. To demonstrate this approach we used a nonlinear programming routine22 interfaced through the GAMS (22) MINOS 5.0, Systems Optimization Laboratory, Department of Operations Research, Stanford University, CA, 1988.
mathematical programming package23 to minimize the residual for the complete methane isotherm10 for two quadrature schemes: one containing 19 quadrature intervals (a case 2 scenario) while the other contains 50 quadrature intervals (a case 3 scenario). In both cases the residual could be reduced to below a solver tolerance of 10-6, and we can therefore conclude there an infinity of PSDs that are compatible with the data exist. Two PSDs that we calculated are presented in Figure 5. In each case the fit to the data is also shown. Figure 5 exhibits an important property that is frequently encountered when calculating discrete pore size distributions. This property is the spikiness of the PSD. (In Figure 5 the first peak in the PSD based on 50 quadrature intervals extends almost to a value of 2 cm3/g Å.) In general, the spikiness of a PSD increases as the number of quadrature intervals included in the analysis increases. This effect can be accounted for by reconsidering the adsorption space. As an example, we reconsider the adsorption space as a function of the number of the pore sizes included in the analysis for the original problem that involved fitting a PSD to two data points. This adsorption space was shown in Figure 4. Note that in each case (i.e., for different number of quadrature intervals), many solutions can be found that involve only two pore sizes. Also, note that all of the column vectors are of a similar order of magnitude. This means that the components of f* in each case will be of the same order of magnitude. However, as the number of quadrature intervals increases, the actual size of the quadrature interval decreases. Therefore, the components of f* are (23) GAMS Development Corporation, 1217 Potomac Street, NW, Washington, DC 20007, 1996.
8240
Langmuir, Vol. 15, No. 23, 1999
Davies et al.
Figure 5. Two PSDs based on 19 and 50 quadrature intervals and their corresponding fits to the experimental data.
significantly magnified as the number of quadrature intervals increase as can be seen from eq 13 and
W-1 ) diag(1/δwj)j)1...m
(22)
Equations 13 and 22 highlight that the spikiness of a PSD is predominantly a mathematical artifact. In other words, the spikiness of a PSD is predominantly related to the number of quadrature intervals used in the analysis and not to any physical attribute per se of the internal structure of an activated carbon. This observation is important in the following section whe we consider how to calculate suitable representative PSDs of the internal structure of adsorbents. The PSDs calculated in this section should not be viewed as representative of the internal structure of an adsorbent without further consideration. Two examples serve to highlight why further analysis is required. First, although we were able to fit a PSD to a full isotherm consisting of 19 data points using a minimization routine, it is significant that a linear solver could not calculate a result at all (based on 19 quadrature intervals) because the problem was too ill-conditioned. A consequence of an illconditioned system is that any small perturbations in the experimental data could lead to substantially different PSDs. Therefore, before any confidence can be placed in the result, we need to incorporate criteria into the solution strategy to ensure that the resulting PSD is stable with respect to perturbations in the experimental data. The second example demonstrates the danger of there being an infinity of possible PSDs. The minimization routine terminates as soon as it finds a minimum for the residual. The actual solution that it reports is therefore strongly affected by the exact algorithm of the routine and starting
point used. Figure 6 shows four substantially different PSDs fitted to the original two data points that we considered. These distributions were obtained by merely changing the initial estimate of the PSD. This example clearly demonstrates the need to determine suitable techniques to identify a representative PSD. Determining a Representative PSD We now consider the problem of obtaining a PSD, based on a given set of data, that is representative of the structure of the real solid. Two requirements of a representative PSD are (i) that it be unique and (ii) that it be relatively insensitive to small perturbations in the data. The PSD must be unique in the sense that there must not be two or more dissimilar PSDs that are compatible with the data. Since unique PSDs can only be calculated in case 1 or case 2, we can immediately conclude that representative PSDs can only be calculated when there are at least as many data points as there are quadrature intervals in the analysis. This implies that although it is acceptable to include more quadrature intervals than data points into an “existence analysis”, this is inappropriate when attempting to calculate a representative PSD. Therefore, representative PSDs cannot be calculated in instances where more quadrature intervals need to be introduced into the existence analysis before at least one feasible PSD can be found. It follows that more data are required in order to calculate PSDs of higher resolution (i.e., PSDs with smaller quadrature intervals). The second requirement arises because PSDs that vary significantly due to small perturbations in the data are not reliable enough to be of practical use. The fact that the amount of data limits the resolution of the PSD is an important reason for adopting a uniform
Calculation of Pore Size Distributions
Langmuir, Vol. 15, No. 23, 1999 8241
Figure 6. Four different PSDs fitted to the original two data points summarized in eq 14 obtained by varying the initial estimate of the PSD.
term that accounts for the smoothness of the distribution
quadrature interval within the region of interest. By defining a suitable “rule” to determine the quadrature we were able to prevent the quadrature spacing being included with the unknown parameters in the analysis. Although a variable quadrature scheme could be adopted, this approach would require almost twice as much data to obtain PSDs of a similar resolution to the ones determined using the fixed quadrature scheme. Having established that only as many quadrature intervals as there are data points should be included when attempting to calculate representative PSDs, we now address the problem of stabilizing the numerical computation. The most commonly used method to stabilize the result is to incorporate additional constraints that are based on the smoothness of the PSD. This method, termed regularization, has been described in detail by Wilson,24 Szombathely et al.,25 Merz,26 and Whaba.27 Regularization stabilizes the result by not allowing the fitting routine to fit the data more closely than is warranted on the basis of the conditioning of the system.26 It also enables us to calculate more realistic PSDs by reducing the spikiness of the result, which has already been shown to be a mathematical artifact. The assertion that smoother PSDs are more realistic corresponds to an argument that we noted in the Introduction: real materials are more likely to exhibit a continuous distribution of pore sizes centered about a few dominant pore sizes. Regularization is incorporated into the PSD analysis by modifying the residual defined in eq 20 to include a
RReg )
(24) Wilson, J. D. J. Mater. Sci. 1992, 27, 3911. (25) Szombathely, M.; Brauer, P.; Jaroniec, M. J. Comput. Chem. 1992, 13, 17. (26) Merz, P. H. J. Comput. Phys. 1980, 38, 64. (27) Whaba, G. SAIM J. Numer. Anal. 1977, 14, 651.
This approximation to the second derivative requires / ) and δwm+1. Following a values for f(w/0), δw0, f(wn+1 prescription similar to that of Wilson,24 we define these quantities as follows:
RReg ) (N - AWf)T(N - AWf) + RS
(23)
where R is a strictly nonnegative smoothing parameter and S is a suitable discrete representation of a function that measures the smoothness of the PSD. Here we follow standard practice and define S as the integral of the square of the second derivative of the PSD. Rewriting eq 23 as a summation using this measure of smoothing we obtain n
∑ i)1
m
[Ni(Pi) -
F(w/j ,Pi)δwj f(w/j )]2 + ∑ j)1 m
f′′(w/j )2δwj ∑ j)1
R
(24)
The second derivative of the PSD can be approximated using a finite difference formula / ) - f(w/j ) f(wj+1 1
f′′(w/j ) )
/2δwj+1 + 1/2δwj 1
-
/ f(w/j ) - f(wj-1 ) 1
/2δwj + 1/2δwj-1
/2δwj+1 + 1/2δwj
(25)
8242
Langmuir, Vol. 15, No. 23, 1999
Davies et al.
f(w/0) ) 0
(26)
δw0 ) δw1
(27)
/ ))0 f(wm+1
(28)
δwm+1 ) δwm
(29)
Reintroducing vector notation, eq 25 can be written as
f′′(w/j ) ) Df
(30)
where the components of D are obtained directly from eq 25. Equation 24 can now be written more compactly as
RReg ) (N - AWf)T(N - AWf) + RfTDTWDf
(31)
Wilson24 has demonstrated that minimizing the residual of eq 31 corresponds to minimizing
RReg ) -2NTAWf + fT(WTATAW + RDTWD)f (32) Equation 32 can be solved for a given smoothing parameter using standard nonlinear minimization routines. However, to calculate a representative PSD, a method is required to determine an optimal smoothing parameter to be used in the analysis. Choosing a suitable method to determine the optimal degree of smoothing is not at all straightforward. Both Wilson24 and Szombathely et al.25 highlight that this is a subjective choice and that different methods may result in quite different PSDs. Therefore, to avoid conclusions that are specific to a particular smoothing criterion, it is worth considering more than one criterion. The three criteria used most frequently are the so-called “L curves”, “aim functions”, and the generalized cross-validation method. L curves (see for example ref 16) are a plot of some measure of the error of the fit to the data against the smoothing parameter. In general, the error of the fit to the data increases as the value of the smoothing parameter increases. However, below a threshold value of the smoothing parameter the increase in the error is often negligible, while above the threshold the error increases rapidly. Plots of the error against the smoothing parameter therefore resemble an “L” lying on its side. These curves are used to identify this threshold value which is taken to be the optimal extent of smoothing. A suitable measure of the error of the fit, E, based on the least squares residual is n
∑ i)1
E(R) ) 1/n
m
(N(Pi) -
F(w/j )δwj f(w/j ))2 ∑ j)1
(33)
The “aim functions”, as described by Szombathely et al.,25 are closely related to the L curve approach. This approach also defines a suitable measure of the error of fit to the data but, instead of identifying the threshold value of the smoothing parameter as the most suitable value, they adopt the value which results in a fitted error that is comparable to an estimate of the experimental error. Szombathely et al.25 have shown that the estimate of the error defined in eq 33 increases monotonically as the smoothing parameter is increased; this criterion thus leads to a unique value of the smoothing parameter. We have found this and similar aim functions are conservative estimates of the optimal smoothing parameter in that they
lead to very smooth PSDs. We have also found that this approach is prone to allowing systematic errors to develop in the fit to the data, especially in regions where low extents of adsorption take place. The observation that a good choice of the smoothing parameter is one that most accurately enables us to predict any one of the n experimental data points from a PSD that is determined using the remaining n - 1 data points leads to the generalized cross-validation method. A significant consideration in applying this technique is that it can usually only be implemented in a straightforward manner when calculating unconstrained PSDs (i.e., PSDs in which the distribution value is not restricted to being strictly nonnegative). Although a suitable form of the generalized cross-validation score function has been derived for the case when the PSD is constrained to being nonnegative,27 the evaluation of this function demands considerable computation.24 For this reason Wilson24 proposed an approximation to the generalized crossvalidation score function that is appropriate for constrained PSDs. Assuming that all the data points are of equal importance, Wilson’s approximation to the crossvalidation function is
GCV(R) ) 1
/n(N - AWf)T(N - AWf)
[1 - 1/nTr(AW(WTATAW + RDTWD)-1WTAT)]2
(34)
where f is determined from eq 32. The optimal smoothing parameter is determined as that one which minimizes the cross-validation function. These three criteria have been used to determine a representative PSD for BPL on the the basis of the complete methane isotherm at 308 K consisting of 19 data points.10 The region of interest was specified, as before, to be the microporous region, and 19 quadrature intervals were included in the analysis. Figure 7 presents the L curve and cross-validation score function for this system. Two appropriate smoothing factors are R ) 10-1 and R ) 10-5, which are based on the L curve criterion and the generalized cross-validation score function, respectively. We note that, even for an extremely high value of R ) 10, the average estimate of the error of the fit is still well below an estimate of the experimental error of about 1.5%.28 The two PSDs based on R ) 10-1 and R ) 10-5 are presented in Figure 8. We also show the PSD based on a smoothing factor of 10 for comparison. Note that the first two PSDs exhibit a pronounced peak at the upper boundary of the region of interest. This indicates that it may be worth extending the region of interest to determine this peak more completely. We have therefore reanalyzed the data by extending the upper limit of the region of interest from 20 to 30 Å. The corresponding L curve and the generalized cross-validation score functions for the PSDs on the basis of the extended region of interest are presented in Figure 9. From Figure 9 we obtained R ) 10-1 on the basis of the L curve and R ) 5 × 10-3 on the basis of the generalized cross-validation score function. The corresponding PSDs are presented in Figure 10 with the PSD-based R ) 10 for comparison. The first two PSDs in Figure 10 are in remarkable agreement. The difference between the third PSD and the two previous distributions highlights the somewhat subjective choice of an appropriate smoothing factor and the ill-posed nature of the adsorption integral (28) Gusev, V. Y.; O’Brien, J. A. Personal communication, 1997.
Calculation of Pore Size Distributions
Langmuir, Vol. 15, No. 23, 1999 8243
Figure 7. L curve and the generalized cross-validation score function for the fits to all the data using 19 quadrature intervals and an upper limit to the region of interest of 20 Å.
Figure 8. PSDs based on all 19 experimental data points and smoothing factors of 10-1, 10-5, and 10. (Maximum limit of the region of interest is 20 Å.)
equation. This is particularly evident in Figure 11 which shows that all three PSDs fit the experimental data extremely well. Finally, we reconsider the confidence that can be placed on the exact locations of the peaks that have just been presented. As we have already noted, Gusev et al.10 introduced the concept of a window of reliability, which identifies the pore size region that can be reliably characterized using specific experimental data. The exact locations of peaks at large pore sizes are difficult to determine due to the fact that the isotherms of the larger model pores (which are calculated using grand canonical Monte Carlo simulation) become essentially indistinguishable. This problem was implicit in our analysis and
can be identified as the nearly linearly dependent columns of the adsorption matrix that corresponded to the adsorption in the larger pores. The almost collinear column vectors can be identified easily in Figure 4. Although we have incorporated regularization into the analysis to stabilize the calculated result it is important to recognize that there is still a certain degree of uncertainty in the exact location any peaks that occur above pore sizes of about 15-20 Å. Confidence in the peaks at these larger pore sizes can be improved by using a more strongly adsorbing adsorbate, measuring the isotherm at lower temperatures, or extending the measurements to include the adsorption at higher pressures.
8244
Langmuir, Vol. 15, No. 23, 1999
Davies et al.
Figure 9. L curve and the generalized cross-validation score function for the fits to all the data using 19 quadrature intervals and an extended region of interest. (Upper limit to the region of interest of 30 Å.)
Figure 10. PSDs based on all 19 experimental data points and smoothing factors of 10-1, 5 × 10-3, and 10. (Maximum limit of the region of interest is 30 Å.)
Conclusions A rigorous procedure to determine whether there are any feasible PSDs that can be fitted to a given set of experimental data has been presented. This procedure is particularly powerful in that it is not limited to determining whether PSDs corresponding to specific functional forms exist; instead it searches the “limiting envelope” of all feasible solutions. Determining a representative PSD based on a given set of data is significantly more difficult than determining whether there are any feasible PSDs that are compatible with the data. Two factors complicate the analysis. The first factor is concerned with being able to attach a statistical significance to the PSD that is calculated. This
is related to the amount of data that are available. The amount of data determines the “resolution” of the PSD that can be calculated in that only as many quadrature intervals as there are data points can be included in the analysis. The second factor is the ill-posed character of the adsorption integral equation. We have demonstrated that standard stabilization techniques that are used in conjunction with the controlled quadrature lead to remarkably consistent characterizations of the internal structure. We have demonstrated that a PSD can be fitted to the adsorption of methane at 308 K onto BPL-activated carbon, which is in agreement with previous work.1,10 The adsorption of methane onto BPL-activated carbon was,
Calculation of Pore Size Distributions
Langmuir, Vol. 15, No. 23, 1999 8245
amounts of data, since a suitable model of the internal structure should be capable of being consistent with a wide range of data. Such tests could involve establishing whether a PSD can be calculated on the basis of several methane isotherms at different temperatures, or isotherms of different pure components or even pure and multicomponent adsorption data. As more data is included in the analysis it becomes more difficult to calculate a PSD. In these instances it is important to be able to establish whether a PSD exists that is compatible with the data (even although it might be difficult to calculate a representative PSD). This is where the techniques presented in this paper become particularly useful. By establishing whether a PSD is compatible with an appropriate range of data we are able to assess the potential utility of PSDs in engineering calculations. The results of such an investigation are reported by Davies and Seaton.1 Figure 11. Fit to the experimental data of the PSDs shown in Figure 10.
however, used only as a case study to demonstrate the techniques presented in this paper. The same techniques can be used to test the realism of a PSD as a model of the internal structure of any activated carbon. This can be achieved by determining whether a PSD for a given activated carbon can be calculated from increasing
Acknowledgment. G. M. Davies gratefully acknowledges financial support from the Bradlow Foundation. The donors of the Petroleum Research Fund, administered by the ACS, are acknowledged for partial support of this work. The authors also thank Matthias Heuchel for helpful comments during the preparation of this paper. LA9902643
