Surface heterogeneity analysis by gas adsorption: Improved calculation of the adsorption energy distribution function using a new algorithm named CAESAR

Surface Heterogeneity Analysis by Gas Adsorption: Improved Calculation of the Adsorption Energy Distribution Function Using a New Algorithm Named CAESAR CORNELIS H. W. VOS AND L U U K K. K O O P A L 1 Laboratory for Physical and Colloid Chemistry, Agricultural University, De Dreijen 6, 6703 BC Wageningen, The Netherlands Received July 30, 1984; accepted October 22, 1984

A new numerical method of solution of the integral equation describing the adsorption of gases on a heterogeneous surface with respect to the adsorption energydistribution function is developed. The derived algorithm, CAESAR (Computed Adsorption Energies, SVD Analysis Result), is based on Singular Value Decomposition(SVD) to control the ill-posedness of the integral equation and negative values of the distribution function are excluded. The experimental error and the estimated error in the local isotherm function are used as criterium to obtain the optimal solution. The method is tested with simulated isotherm data and the effects of random and systematic errors are investigated. Application of the method to N2 adsorption isotherms on four carbon blacks, differing in surface heterogeneity, shows that the method is capable of detecting small differences in adsorption energies. A comparison with results obtained with CAEDMON (Computed Adsorptive Energy Distribution in the MONolayer) and HILDA (Heterogeneity Investigated by a LoughboroughDistribution Analysis) (or ALINDA: Adamson LINg Distribution Analysis) indicates that CAESAR yields more reliable results.

© 1985 Academic Press, Inc.

INTRODUCTION Ordinary solids usually have a surface that is energetically heterogeneous and this heterogeneity will show up in most interfacial properties of such a surface. Of the numerous reasons for heterogeneity we mention local variations in surface structure, variations in crystal structure, nonuniform distribution of surface groups, crystal defects, surface impurities, porosity, etc. In the last decade the role of surface heterogeneity with adsorption has received considerable attention and adsorption models for heterogeneous surfaces have been developed (1, 2). Such models yield the adsorption isotherm when the adsorption energy distribution is known or, vice versa, the adsorption isotherms can be used to calculate the distribution function. Several numerical methods (3-10) have been developed to calculate the adsorption energy distribution. 1TO whom correspondence should be addressed.

A major problem with the existing numerical methods is that the obtained distribution function depends on an arbitrarily chosen parameter which controls the detailedness of the distribution function (11). With H I L D A (Heterogeneity Investigated by a Loughborough Distribution Analysis) (5, 6) presmoothing may be necessary and the distribution function is determined by the number of iterations allowed, with C A E D M O N (Computed Adsorptive Energy Distribution in the MONolayer) (3, 4) the number of "patches" controls the smoothness, and with the regularization methods (7-10) the Lagrange multiplicator controls the smoothness. Another disadvantage of the regularization method developed by House (7, 8) is that the distribution function can obtain negative values. In the regularization method of Brown et aL (9, 10) negative values of the distribution function are excluded. These authors also discussed a method to obtain a value for the Lagrange multiplicator, but the results are not encouraging.

183 0021-9797/85 $3.00 Journal of Colloid and Interface Science, Vol. 105, No. 1, May 1985

Copyright © 1985 by Academic Press, Inc. All fights of reproduction in any form reserved.

184

VOS AND KOOPAL

In this study a new numerical computation method will be presented. The main advantage of this new procedure over the existing ones is that the parameter which controls the fine structure in the calculated distribution is related to the experimental error and the adequacy of the local isotherm function with respect to the system under study. Moreover, negative values in the energy distribution are excluded. The new method will be tested with simulated isotherm data calculated from an analytical expression for the adsorption energy distribution function and the local isotherm equation. To illustrate its value in the case of experimental isotherms we analyzed adsorption data of N2 on four carbon samples (12). The samples differ in degree of graphitization and consequently also in the heterogeneity of the surface. The analysis of these carbons is especially interesting, because it is known that progressive graphitization renders the surface more homogeneous. A preliminary description of the new method, including results for four chemically modified Aerosils has been given elsewhere (13). DESCRIPTION OF THE ANALYSIS

The Integral Adsorption Equation A heterogeneous surface can be envisioned as consisting of patches or regions which individually are energetically homogeneous. We take f(U)dU as the fraction of sites having adsorption energies between U and U + dU and O(U, P) as the local isotherm function. The latter represents the fraction of sites covered, 0, at equilibrium pressure P and energy of adsorption U. Then the coverage, 0T, of the entire surface is described by the equation (2)

O~(P) =

f

Umax

O(U, P)f(U)dU,

[1]

t/Urnin

where U m i n and Umax are the lower and higher extreme of the energy of adsorption of the system under consideration. The integral should be multiplied by Vm (cm 3, Journal of Colloid and Interface Science, Vol. 105, No. 1, May 1985

STP), the monolayer volume, when the overall isotherm is expressed in adsorption units (cm 3, STP). It should be noted that the adsorption energy is a system property and not just a surface property. In other words, the adsorption energy depends also on the adsorbate. In principle the adsorption equilibrium is governed by a free energy of adsorption instead of by the energy of adsorption. In the case of adsorption from the gas phase the thermal entropy contribution to the adsorption free energy is small and to a good approximation the free energy may be replaced by the energy of adsorption as is done above. The local isotherm function to be inserted in the integral adsorption Eq. [1] can be any of the existing adsorption isotherm equations derived for homogeneous surfaces (14). In most of these equations lateral interactions are taken into account. When this is the case, the type of surface heterogeneity determines the way in which the local isotherm function is incorporated in Eq. [1]. Usually, two extremes are considered: patchwise and random heterogeneity. On a patchwise heterogeneous surface the lateral interactions are counted per patch only, interactions between molecules on different patches are neglected. With a random heterogeneous surface the lateral interactions are smeared out over the entire surface. In the analysis of surface heterogeneity from adsorption data a realistic choice for the local isotherm is important. Unfortunately, for most adsorbent-adsorbate combinations, data on the adsorption behavior on the corresponding homogeneous surfaces are not available. Therefore, an arbitrary choice has to be made. For ordinary surfaces such an assumption will usually be less serious than the often-made presumption that the surface is homogeneous, whereafter a (local) isotherm is fitted to the data. In the present work we assess the consequences of a wrong choice of the local isotherm. Moreover, it will be shown that the imperfections of the local isotherm

SURFACE HETEROGENEITY

A N A L Y S I S B Y GAS A D S O R P T I O N

equation can to some extent be accounted for in the determination off(U). A second problem is that Eq. [ 1], a Fredholm equation of the first kind, is ill-conditioned with respect to solution o f f ( U ) . The first numerical methods to calculate f ( U) (35) have not addressed special attention to this problem. In the regularization methods (7-10) due attention is paid. Also in the present method the ill-posedness of Eq. [1] is explicitly considered. Singular Value Decomposition For the numerical evaluation Eq. [1] has to be approximated by a summation:

i=1

1,...,m)

volves decomposition of matrix A according to A = UQV t [4] in which U i s a m × n matrix, V t is the transpose of the n × n matrix V, and Q = diag(q~, q2 . . . . .

qn),

[5]

where qi's are assumed to be in decreasing order, qi's are called the singular values of A. If A is exact and of rank r (that is to say, r linearly independent rows or columns) the values of qi are qi > 0

for

i = 1, 2 . . . . .

q~=0

for

i=r+l,r+2

r .....

n.

In practice A is never exact or, qi will never be exactly zero and only the sequence of the singular values can be obtained

OT(Pj) = U ~ O(U~, Pj)f(U~)AU (j=

185

[2]

in which n (~
ql >~ q2 >1 q3 >~" • ">~ qn. Consequently, the numerical rank r of A has to be defined in practice as the n u m b e r of rows or columns of matrix A which m a y be considered as linearly independent in the numerical analysis. The value of r will depend on the actual situation, that is, on the machine precision, the experimental error in the measured isotherm, and on the quality of the local isotherm function. Once the rank of A has been established, A can be replaced by a nonsingular matrix Ap using the SVD procedure, mpequals A v = UQpV t

[6]

in which

as

0T = Af,

[3]

where 0T is a vector with elements 0T(Pj), f is a vector with elements f(Ui), and A is a matrix with elements aj.i = O(Ui, Pj)AU. Owing to the ill-posed character of Eq. [1] the matrix A is near-singular. Inversion of A using a standard method is therefore not leading to a unique solution of A -~. Special precautions are required to overcome this problem. We applied the so-called Singular Value Decomposition method (15). This in-

Qp = diag(ql, q2 . . . . .

qr, oo . . . . .

oo).

In practice the largest singular value q~ is used instead of oo. The difference between A and Apis that too small singular values have been replaced by a very large value. As a result the matrix Ap is nonsingular. The equation OT = Apf [7] can now in principle be solved with respect to f using a standard method. Use is made of a nonnegative least-squares (NNLS) proJournal of Colloid and Interface Science, Vol. 105, No. 1, May 1985

186

VOS AND KOOPAL

cedure (16) which has the advantage that all elements of f will be positive. Morrison and Sacher (4) have given a brief discussion of this NNLS procedure and used it in the refined version of CAEDMON. For sake of clarity it should be remarked that the SVD procedure also can provide an inverse of Ap: Ap 1=

VQ~lU t

[8]

Ap ~ can be used to obtain f: f

= AplOT

[9]

.

In general this solution may contain negative values o f f ( U i) which are physically unrealistic. To avoid this difficulty the SVD method is only used to replace A by Ap, whereafter f is solved using the NNLS procedure.

Rank of Matrix A The remaining problem is to determine the numerical rank of matrix A. This is done following a method proposed by Hanson (17). As stated above, the actual rank of A will depend on the numerical inaccuracy, the experimental error in 0T and on the approximation that is inherent to the choice of a specific local isotherm function. For the moment, it will be assumed that the machine precision is very high and that the local isotherm equation is exact. In that case the rank of A should be chosen in accordance to the experimental error in 0T. To this order f is first determined for a given minimal rank using Eq. [7]. The obtained value of f, f~, is now multiplied by the original matrix A to obtain the calculated overall isotherm 0T,c =Afc.

[101

To check this solution the deviations between 0T(Pj) and 0Tx(Pj) have to be considered. The solution f is required such that

~ ]OT(PJ)~ OT'c(PJ) = m, j=l

[11]

m would be asking for more accuracy than is inherently present in 0T(Pj). The rank of A is found by an iterative calculation. In each step r is increased by one until the value of r is found which corresponds with Eq. [11] or when a specified maximum is reached (usually r = n). To account for the machine precision (~ ~ 10 -19) the restriction has been made that qr should not be less than 10 - 6 . The latter value is established by trial and error using simulated isotherm data obtained for several distribution functions. In Eq. [ ! 1] the absolute difference between 0T(Pj) and 0T,c(Pj) divided by aj is used, whereas Hanson suggested to use the square of this difference divided by the variance. Small and large errors are weighted more equally by Eq. [11] and in practice Eq. [11] turned out to give better results than Hanson's approach. In practical situations the local isotherm function is only an approximation of reality. Therefore, it is appropriate to include in oi also the error occurring in Oxx owing to the approximate nature of 0(U, P). The errors in OT and Ox,~ are independent and operate simultaneously. If it is assumed that both deviations are random, then the overall deviation, aj, can be written as o-j = [{0"j(0T)}2 -4- {O'j(0T,c)}2]°5

[121

with aj(0x) accounting for the error in 0T(Pj) and aj(0T,c) for the error resulting from the approximate nature of the local isotherm. In principle, aj(0T) can be obtained when the experimental error in 0T(P) is known. The deviation aj(0T,c) is difficult to establish and it may not be a random but a systematic deviation. Below a way to obtain an estimate of an average value of a(0T,c), which is the same for all j, will be discussed. Due to the averaging a(OT,¢) should be treated as a random deviation.

O'j

CAESAR where, in this case, ai is the mean deviation of 0T(Pj) due to experimental errors. Demanding the summation to be smaller than Journal of Colloid and Interface Science, Vol. 105, No. 1, May 1985

In order to perform the calculations a computer program named CAESAR (Com-

187

SURFACE HETEROGENEITY ANALYSIS BY GAS ADSORPTION

puted Adsorption Energies, SVD Analysis Result) has been written in Simula. CAESAR offers five options for the local isotherm. For mobile adsorption on patchwise heterogeneous surfaces the Hill-de Boer (HB) and virial equation are available. The HB equation can be written as (14) 0 1-0

(0 exp 1 - 0

CO)

= krnP e x p ( - g / R T ) ,

[131

where C is the lateral interaction parameter and km is a constant depending on the vibrational frequencies of the adsorbed atoms. The lateral interaction parameter C can be derived from three-dimensional gas properties (18, 19) if it is assumed that these interactions are an adsorbate property only. For C = 0 the Volmer equation is obtained. The virial equation developed by Ross and Morrison (3, 18) and written in closed form is 2O ln0+---21n(1 1-0

-0)+2(B*-2)0

+ ~3( C * - 2 ) 0 2 + g 4 ( D * -- 2)0 3 + ¼(E* - 2)04 + 6(F* - 2)05 = ln(kmP),

[141

where B* to F* are reduced virial coefficients. In the case of N2 at 77.5 K these coefficients are

where C' is the lateral interaction parameter and kL a constant characteristic for localized adsorption. The F F G equation can be applied for patchwise and random heterogeneous surfaces by taking the lateral interactions per patch only (0" = 0) or by averaging the interactions over the entire surface (0" = 0T). C' can be approximated on the basis of threedimensional gas properties of the adsorbate and the lattice structure of the surface (19). With C' = 0 the F F G isotherm reduces to the Langmuir equation. In this case the type of heterogeneity is irrelevant. RESULTS AND DISCUSSION

Simulated Isotherm Data In order to test the reliability of the method, isotherm data based on an exactly known local isotherm and distribution function are required. As there are no experimental results available for which both requirements hold, it is necessary to test the methods with generated isotherm data. To this order we used the HB equation with C = 5.5 as local isotherm and a bi-Gaussian distribution function: 1 (Odl

f(U) = ~

O/2~ -1

-~l + -~2] [a,exp/31(U1

- Uo,1)2 + c~2exp f12(U2 - U0,2)2], [16] B* = - 1.9346 C* = 3.43929 D* = 2.0646 E* = 2.0322 F* = 1.9050.

Again, these values are ideal 2 - D coefficients and take no account of the gas-solid perturbation terms. For localized adsorption the Langmuir and F r u m k i n - F o w l e r - G u g g e n h e i m (FFG) equations (14) are available. The F F G equation can be written as

0

- 1-0

exp(C'0*) =

kLP exp(-U/RT),

[15]

where Uo,i is the position of the ith maximum, ai is the relative weight of peak i, and fli is a measure of its width. The overall isotherm for a patchwise heterogeneous surface was calculated using U0,1 = 0.80, U0,2 = 3.80, al = a2 = 0.5, fll = 6.0, and r2 = 3.0. The obtained 21 data points, OT(Pj), are equally distributed on a logarithmic pressure scale. F r o m this isotherm and the known local isotherm equation the distribution function has been recalculated. In Fig. 1 the original and the calculated distribution function are compared for the case that both 0(U, P) and 0T(P) are without error. In this case the numerical inaccuracy determines the value of qr. As long as qr >/ 10 - 6 the rank r equals its m a x i m u m value, Journal of Colloid and Interface Science,

Vol. 105, No. 1, May 1985

188

VOS AND KOOPAL g

f(U) was evaluated. Results are shown in

f(u)

.00

1.00

-2."00- -

3.00

U

4.00

5.00

FIG. 1. Comparisonof the computed adsorption energy distribution function (- - O - -) with the true distribution function (--). The overall isotherm, O-r(P),and the local isotherm, O(U,P), were free of errors. In all figures the adsorption energies are in units of RT.

Figs. 2a to d. The standard deviation in the isotherm data varied from 0.1 to 10%. In the case of a deviation of 0.1% the distribution is still very well reproduced. Deviations of 1 to 5% lead to a flattening of the high-energy peak. In the case of a deviation of 10% the distribution curve is smoothed considerably, however the position of the peaks is well preserved and valuable information on f(U) is still obtained. Table I shows the values of the mean relative deviation, 6, between 0 T and 0T,c, where 6 is defined as

Em iOT(Pj ) _ 0 T , c ( 5 ) [17] 6=5 I • r = n. This is almost always the case if a mj=t local isotherm is used in which lateral interactions are taken into account. For the Lang- It follows that 6 is increasing with the expermuir or the Volmer isotherm qr values smaller imental error, this is an inherent property of than 10-6 frequently occur. This is due to our calculation method. When large experithe fact that in these equations P and O(U, mental errors are present it is vicious to P) are almost linearly dependent for small pursue great accuracy in the calculated disvalues of O(U, P), consequently matrix A tribution function. The distribution function becomes very ill conditioned. Figure 1 shows shown in Fig. 2d is the best attainable result that the agreement between the calculated in the case of a standard deviation of 10% in and the true distribution function is excellent. the experimental data. A better fit would This is partly due to the fact that negative lead to details in f(U) which would be due values of f(U) have been excluded, f(U) to experimental errors in 0V, rather than to values calculated with an ordinary least- physical characteristics of the system. In Table squares method, which allows f(U) to be I 6 is notated as 6(0T) in order to emphasize negative, did not lead to such good results: that 6 is the result of an error in 0T. Another effect of the error in OTis that the the region between two peaks became negacalculated specific surface area, Asp, differs tive and the peak height decreased. If the from its true value. Asp can be calculated Langmuir equation governs the local isousing N (see Eq. [2]) and assuming a crosstherm, the agreement between the calculated sectional area for the adsorbate molecules. and true distribution function is also not as perfect as in Fig. 1. It m a y be concluded that In Table I the error in Asp is tabulated. It is in the presence of lateral interactions between not surprising that a large error in OT leads the adsorbate molecules the true distribution to a less accurate value of Asp. The error in function can be found with CAESAR, pro- Asp is mainly determined by the experimental vided the experimental isotherm is free of error in the adsorption data close to the monolayer coverage. errors. Thus far we considered randomly distribIn practice experimental errors will always uted errors in 0T(P). In the case of systematic be present. To imitate this situation normally errors, the obtained distribution function will Gaussian distributed errors have been superimposed on the values of OT(Pa) (Pj values also deviate systematically from the true were assumed to be accurate), whereafter function. To illustrate this the distribution Journal of Colloid and InterfaceScience,Vol. 105, No. 1, May 1985

SURFACE HETEROGENEITY ANALYSIS BY GAS ADSORPTION

0.1%

189

1%

f(U)

f(U)

m

--1.50

.00

1.50

3.00

U

4.50

6.00

-1,50

.00

1.50

3.00

4.50

6.00

U' 4.50

6.00

~U

10%

5%

flu)

f(U)

--1,50

.00

13 0

3.00

U

4.50

6.00

--1.50

.00

1.50

3.00

FIG. 2. Effect of the presence of random errors in the overall isotherm data on the calculated adsorption energy distribution function. The magnitude of errors is indicated.

function was recalculated using various local isotherm equations. In this way we simultaneously show how a wrong choice o f the local isotherm affects the results. In the first place we varied the lateral interactions parameter C in Eq. [13]. That is to say, the distribution function has been calculated with C -- 6.75 and 3.9 instead o f with C -- 5.5. In these calculations we used the exact values of 0T(Pj) (no experimental error). Consequently the rank o f A is determined by the machine precision. With the TABLE I The Mean Relative Deviation 6(0T)and the Error in A~r, as a Function of the Experimental Error in 0T (%)

(%)

(~o)

0.0 0.1 1.0 5.0 10.0

0.002 0.030 0.277 1.39 3.41

o. 1 0.2 0.2 3.0 3.7

HB equation as local isotherm it turns out that r = n. The resulting distribution functions are shown in Figs. 3 and 4. Figure 3a shows that overestimation of the lateral interactions leads to a flattening o f f ( U ) and to a shift o f the peaks to lower energy values. This is due to the fact that part o f the adsorbate-surface interactions is assumed to be of lateral nature. In the case o f an underestimation o f C, the assumed local isotherm increases less with pressure than the true local isotherm and reaches a lower pseudosaturation level. As a consequence, the obtained energy distribution is shifted to higher energy values and the peaks have become more narrow (see Fig. 3b). Instead o f using a wrong value for the lateral interaction, it is also possible that a wrong choice is m a d e for the local isotherm function. Figures 4a and b show the effect o f this type o f error. The distribution function in Fig. 4a is calculated with the virial equation for N2 as local isotherm and in Fig. 4b f ( U ) is evaluated using the F F G equation with Journal of Colloid and Interface Science, Vol. 105, No. l, May 1985

190

VOS AND KOOPAL

b

i'

m

f(U)

iI II

flu)

'I

I I I I I t

8

~d ~

-1.50

I'

3,.=

.00

1.50

i'~ 3 O0

'~' U

4.50

-1.50

6.00

.00

1.50

3.00

U

4.50

6.00

FIG. 3. Effect of systematic errors in the local isotherm on the calculated adsorption energy distribution. As local isotherm the HB equation with C = 6.75 (a) and C = 3.9 (b) is used, instead of with C = 5.5 the correct value for the lateral interaction parameter.

lateral interaction parameter C' = - 4 . 1 . Both this value and C = 5.5 in the H B equation are reasonable values for adsorption o f N2 (19). Figure 4a shows two narrow peaks instead o f the true flU). Apparently the lateral interactions are underestimated. This conclusion is supported by the fact that the virial isotherm on an h o m o g e n e o u s surface corresponds best with the H B equation with C = 3.8 (22). Use o f the F F G equation splits the peak at high energy values in two or three peaks and flattens the peaks somewhat. F r o m these examples it follows that systematic errors affect the shape o f the obtained distribution curves in various ways and that an accurate prediction is not feasible if the local isotherm function used differs too m u c h f r o m the true 0(U, P). Yet, the basic features

e,i

.o

;;

a

o

o f the distribution function are present: in all cases two m a i n peaks were f o u n d at about the correct position on the energy interval. Systematic errors, like experimental ones, can lead to deviations in the calculated value for Asp, Table II illustrates this. Overestimating the lateral interaction leads to too small values o f Asp, whereas underestimation gives rise to too large values o f Asp. The variations in Asp(calc) are m u c h larger in the case o f systematic errors than for r a n d o m experimental errors in 0T, even if the latter a m o u n t to 10% (compare with Table I). In Table II also values o f ~ (see Eq. [17]) are given. In this case 6 is a measure o f the m e a n relative deviation introduced by the error in the local isotherm function. We therefore denote this ~ as 6(0). W h e n the

,~

+1 jI I I jt

e,i f{U)

==

- I .50

t;

II

ii

I L

I I

[ I

.00

1.50

3.00

U

b f(u)

4.50

6.00

-1.50

.00

1.50

3.00 U

4.50

6.00

FIG. 4. Effect of a wrong choice of the local isotherm on the calculated adsorption energy distribution. (a) The virial equation (Eq. [ 14]) with coefficientsof N2 is used instead of the HB (C = 5.5) equation. (b) The FFG equation (Eq. [15]) with C = -4.1 is used instead of the correct local isotherm. Journal of Colloid and Interface Science, Vol. 105, No. 1, May 1985

SURFACE

HETEROGENEITY

ANALYSIS

TABLE II Relative Deviation in Asp and the Mean Relative Deviation 5(0) Owing to a Wrong Choice of the Local Isotherm

HB, C = 5.5 (correct O(U, P)) HB, C = 6.75 HB, C = 3.9 Virial equation (N2) FFG, C' = -4.1

AA,~

6(0)

(%)

(%)

0 -3 + 13 +33 -22

0.002 0.163 2.079 1.394 0.869

correct isotherm (HB, C = 5.5) is used 6 is very small, whereas a 100- to 1000-fold increase in 6 is observed if the wrong local isotherm function is chosen. Results obtained with the virial equation compare well with those using HB, C = 3.9, as expected. The HB equation with strong lateral interactions (C = 6.75) leads to a relatively small value of 6. In the extreme case of the condensation approximation, where the local isotherm is replaced by a step function (20), the fit between 0T and 0T,o can be made excellent. In this physically unrealistic situation the patches are filled subsequently and f ( U ) can be adjusted such that 6 will depend on the numerical inaccuracy only. In general, in the absence of experimental errors, 6 will be gradually decreasing if the lateral interaction, C, as accounted for in the local isotherm increases. When the correct function for 0 is used, the 6(C) function must have a very narrow minimum at the optimum value of C (see Table II). In practice it is very difficult to establish this narrow minim u m because of the presence of the experimental errors and the approximate nature of the isotherm equation used. Setting aside these practical difficulties, optimization of 6 offers in principle a means to obtain the correct local isotherm function. When the optimization procedure leads to a lateral interaction which is stronger than that to be expected on the basis of three-dimensional behavior of the adsorbate gas, the result should be considered with reservations. Provided that it is well known, the specific

BY GAS ADSORPTION

191

surface area may constitute a further selection criterium for the local isotherm function. Too small values of Asp indicate an overestimation of the lateral interaction, too high values indicate an underestimation. Comparing 6(0), due to errors in 0, with 6(0~), and due to errors in 0T (see Tables II and I, respectively), shows that 6(0) easily outweighs 6(0T) when the error in 0T is small In other words, if the experimental error in 0T is small, then large values of 6 indicate that the local isotherm function used is rather approximate. In that case 6 is a measure of the average relative deviation introduced by the error in the local isotherm equation. In the next section, where experimental data will be analyzed, the adequacy of the local isotherm equation will be assessed in this way. The result is taken into account in the determination of r. Hence, it influences the maximum detail in f ( U ) that may be considered realistic. Apart from the reported distribution functions several other bi-Gaussian and singlepeaked functions have been studied briefly. Similar results as reported here were obtained. Distribution functions with more than two peaks have not been investigated. The illposedness of Eq. [1] is more likely to be due to too much smoothness of O(U, P) than to the nature off(U). We therefore believe that CAESAR controls ill-posedness of Eq. [1] also i f f ( U ) contains several peaks.

Experimental Isotherms of N2 on Carbon Black Surfaces Excellent surfaces to test the heterogeneity analysis are graphitizable carbon blacks, such as Spheron and Sterling (CABOT Corp.). By heating for several hours at 2700°C in an inert atmosphere these carbons are graphitized to a high extent, resulting in a homogeneous or nearly homogeneous surface (13, 21). Sterling MT-2700 attains the highest state of graphitization, due to its large particle size. Both Sterling MT-2700 and Spheron2700 (usually called Graphon) are frequently used as standard uniform surfaces. Journal of Colloid and Interface Science, Vol. 105, No. 1, May 1985

192

vos AND KOOPAL

Another advantage of these carbon blacks is that detailed low-pressure N2 gas adsorption isotherms at 77.5 K are available (13) on Spheron, Spheron-800 (heated for 12 h at 800°C), Graphon, and on Sterling MT-2700. In order to analyze these isotherms a decision has to be made regarding the local isotherm. There is evidence that the adsorption of N2 on Graphon and Sterling MT2700 at 77.5 K is mobile (13, 22). Both the HB and the virial equation (18) apply to this situation. The virial equation is physically more sophisticated than the HB equation. The second and third virial coefficients are based on a Lennard-Jones interaction potential. The HB equation is derived from the two-dimensional van der Waals equation of state, it incorporates pair interactions on the basis of a Sutherland potential. The constants in both equations can be obtained from three-dimensional N2 gas properties (18, 19). It should be noted that, by using 3-D gas properties to obtain these constants, effects of gas-solid perturbations have been assumed negligible. Because of the mobile character of the adsorbate, we regarded the surfaces as patchwise heterogeneous. The adsorption of N2 on carbon blacks is not restricted to monolayer coverage. Multilayer adsorption may take place on highenergy sites before completion of the monolayer. To correct for multilayer adsorption a procedure has been developed (23) based on multilayer adsorption theories (24, 25). For the present results, the corrections were small if only adsorption values below the BETmonolayer capacity were used. Therefore, multilayer corrections have not been applied. The upper and lower limit of the energy interval to be considered is initially estimated by means of the condensation approximation (20, 3). In the case of Spheron and Spheron800 this interval has been enlarged somewhat in the definite calculations, whereas in the case of Graphon and Sterling MT-2700 ultimately a somewhat smaller interval than the condensation interval has been used. It should be realized that the choice of the actual energy interval is always somewhat Journal of Colloid and Interface Science, Vol. 105, No. 1, May 1985

arbitrary. The theoretical limits to the energy range, viz. U m i n = 0 and U r e a x : 0(3, usually are impractical. The interval obtained by the condensation approximation should be seen as a first approximation of a practical energy range. In the case of a wide distribution it is likely that the condensation approximation interval is somewhat too small, leading to distortions o f f ( U ) near the energy extremes considered. With Spheron and Spheron-800, both having relatively wide distribution functions, the distortions at the ends of f(U) could be considerably reduced by taking the energy interval somewhat larger than the condensation approximation interval. For Graphon and Sterling MT-2700, both having a rather narrow energy distribution, a wide energy range makes the solution o f f ( U ) less stable, because at most values of U in the energy range f(U) is close to zero. In such a case limiting the energy range to values below the condensation approximation range gives smaller distortions at the ends of the curves. In general it will be a good practice to check upon the distortions at the ends of the curves by allowing some variation in the condensation approximation interval. The rank of the matrix A has been established as outlined in the previous sections. The experimental error in the measured isotherms is small, on average about 0.3%. Hence, the error due to the approximate nature of the local isotherm can be estimated under the assumption that the experimental isotherm is correct. To calculate cr(0x,c) firstly the rank of A is determined. Because we consider for the moment that 0x(Pj) is without errors, r is determined by the machine precision and can be calculated. Having r and the given local isotherm equation ff(OT,c) c a n be calculated as 1

~(0T,o) = -

m

~

[0T(ej) -

0~,c(~)l.

mj=l

In the present case the adsorption is given as volume V (cm s, STP) and the monolayer capacity, Vm (cm 3, STP), is obtained in the

SURFACE HETEROGENEITY ANALYSIS BY GAS ADSORPTION

calculation as the normalization constant. Therefore, the error o(V~) is obtained instead of a(0Tx). 1

m

o-(Vc) = - - Z VmIOT(Pj)- 0T,~(Pj)[, m j=l

where Vc is defined as Vm0w,c(Pj). The thus obtained value of ~r(V~)is an average quantity independent of Pj. Relative errors ~r(Vc)/Vm as obtained for the four carbons using either the HB (C = 5.5) or the virial equation (Eq. [14]; N2) as local isotherm are shown in Table III. The virial equation gives slightly larger values for cr(Vc) than the HB (C = 5.5) equation. Due to the greater lateral interaction accounted for in the HB (C = 5.5) equation, the values of a(V~) for the HB should be considerably smaller than those obtained using the virial equation. However, this is not true. Therefore, and because the virial equation is based on a slightly more sophisticated physical model, the latter equation is mainly used in the calculations presented below. The experimental errors essentially follow from the experimental accuracy. They depend on the pressure and adsorption value attained and vary from 0.1 to about 0.7%. On average the experimental value is about 0.3%. It should be realized that a part of the experimental error is incorporated in ~r(V~), owing to the manner in which a(V~) has been calculated. This part is small if the experimental error is much smaller than the error due to the choice of O(U, P) and large in the

TABLE III Estimated Mean Relative Errors, a ( V c ) / V m , that Originate from the Approximate Nature of the Local Isotherm ~r(Vo)lVm (%)

sample

HB (C = 5.5)

Virial

Spheron Spheron-800 Graphon SteflingMT-2700

2.85 4.51 2.52 2.64

4.42 4.67 3.77 3.95

193

case of an accurate local isotherm. Hence, the values of o-j(V) to be introduced in Eq. [12] should be smaller than the actual value of aj(V). In the present case where a(Vc)/Vm is roughly 5 to 10 times as large as the relative experimental error, aj(V)/Vj, we have assumed that about 15% of the experimental error is accounted for in o(Vc). Therefore the remaining average relative experimental error is about 0.25%. For sake of simplicity this average relative error has been used in the further calculations. The overall error is obtained according to Eq. [12] as

o~ = [{~j(v)}2 + {~(v~)}2] °-5. The thus obtained values of ~rj are used to find the final rank of A, whereafter the calculation o f f ( U ) is completed. The choices made for the calculation of at introduce a certain element of arbitrariness. Nevertheless, the advantage of the suggested procedure still is that the rank of matrix A is connected to the overall error. When different data sets are analyzed the method ensures that this is done in a consistent way. Distribution functions calculated with CAESAR using the virial equation as the local isotherm are presented in Figs. 5 and 6. Table IV gives the various parameters used and calculated. The final rank r is in all cases higher than the rank which was obtained using the HB equation. This suggests that the virial equation indeed is a better local isotherm. The values obtained for Vm correspond very well with those obtained by Wesson (12). From electron microscopy it is known that graphitization at 2700°C decreases the particle size. The present analysis confirms this. Figure 5 shows the energy distribution functions in the case of Spheron and Spheron800. Spheron-800 is slightly less heterogeneous than Spheron. Graphitization of the entire particles usually starts above 1000°C, however it is not surprising that the surface structure is already affected somewhat below 1000°C. In Fig. 6 the distributions of Graphon and Sterling MT-2700 are shown. Clearly these surfaces are much more hoJournal of Colloid and Interface Science, Vol. 105, No. 1, May 1985

194

VOS AND KOOPAL

f(U) Spheron

I

8

10

12

14

16

18

20

U

22

Spheron

8O0

f(O)

10.00 12.00

14.00 16.00

U

18.00 20.00

FIG. 5. Adsorption energy distributions of N2 adsorbed (77.5 K) on Spheron and Spheron-800 calculated with CAESAR using Eq. [14] as local isotherm. The adsorption energies are given in units of R T , at 77.5 K R T equals 0.64 kJ/mol.

mogeneous than Spheron-800. The energy distribution in the case of Graphon still contains two large peaks, about 1 R T apart on the energy axis. This indicates that the surface is not completely homogeneous. The distribution of Sterling MT-2700 presents one small and one very large peak. Apparently the Sterling MT-2700 surface is more homogeneous than that of Graphon. The positions of the peaks are very similar for both surfaces, as expected. The adsorption energy in the case of Sterling is just in between the energies of the two peaks of Graphon. In another paper (13) the same carbons have been analyzed using the HB (C = 5.5) equation as local isotherm. The energy distribution of Spheron as obtained with both local isotherm equations is very much the same, but the peaks in the "HB-distribufion" are somewhat lower. For the other surfaces the HB-distributions are somewhat wider and the peaks are lower than for the "virialJournal of Colloid and InterfaceScience, Vol. 105,No. 1, May 1985

distributions." This is due to the fact that the lateral interaction as accounted for in the HB (C = 5.5) equation is larger than that in the N2-virial equation (compare with Fig. 3). As an example, the HB-distributions for Graphon and Sterling MT-2700 are shown in Fig. 7. In the case of Sterling MT-2700 the virial-distribution shows one large peak (Fig. 6) whereas in the HB-distribution this peak is replaced by two smaller peaks (Fig. 7). Accepting that Sterling MT-2700 is a homogeneous surface supports to our conclusion that, in the present case, the virial equation gives a better description of the local isotherm than the HB (C = 5.5) equation. The same experimental data have also been analyzed by Wesson (12) using the first version of CAEDMON. In the present case the distributions show considerably more detail and reveal the differences between the samples much better. Sacher and Morrison (4) using the revised version of CAEDMON analyzed the data of Spheron, Spheron-800,

g f(U)

i Graphon

12.00

14.00

3]'1

16.00 U 18.00

20.00

Sterling MT

f(U},~.~ !

g 12.00

14.00

16.00

U

18.00

20.00

FIG. 6. Adsorption energy distribution functions of N2 adsorbed (77.5 K) on Graphon and Sterling MT-2700 calculated with CAESAR using Eq. [14] as local isotherm.

SURFACE HETEROGENEITY ANALYSIS BY GAS ADSORPTION

195

TABLE IV Parameters Used in the Calculations off(U) and the Resulting Values of ~ and Vm

Spheron Spheron-800 Graphon Sterling MT-2700

m

n

Um~.

U~

r

(%)

(cm3, STP/g)

39 60 56 49

39 60 56 49

9.0 9.5 12.0 10.9

21.0 19.5 18.6 19.3

25 38 20 27

0.33 0.54 1.18 0.95

26.7 25.5 30.2 3.29

and Graphon. The results compare reasonably well for Spheron and Spheron-800, but not so well for Graphon. This is partly due to the fact that a small n u m b e r of "patches" has been chosen for the calculation with C A E D M O N . Automatically only a small n u m b e r of points of the distribution function can be calculated, but the solution is stabilized against the ill-posedness of Eq. [ 1] (11). House and Jaycock (5) have used the experimental data on Spheron and Spheron800 to calculate the distribution functions

f(u)

12.00

14.00

16.00 U 18.00

20.00

f(U) (

12.00

14,00

16.00

U

18.00

20.00

FIG. 7. Adsorption energy distribution functions of N2 adsorbed on Graphon (a) and Sterling MT-2700 (b); calculated with CAESAR but with the HB (C = 5.5) equation as local isotherm (compare with Fig. 6).

with HILDA. The distribution function for Spheron-800 as calculated with H I L D A (using the virial equation) is very similar to that given in Fig. 5. For Spheron more detail is shown. This is due to the fact that with H I L D A the difference between 0T and OT,cis minimized, regardless of the errors due to the experimental inaccuracy and the local isotherm. Consequently, some unrealistic fluctuations in flU) m a y appear. In (5) distribution functions were also calculated using the HB (C = 5.5) equation as local isotherm. In that case both distributions tend to show more fine structure than those obtained using CAESAR with the HB-equation. The HBdistributions obtained with H I L D A also show more fine structure than the virial-distributions found with HILDA. Probably this is due to the fact that a larger lateral interaction leads to a closer fit between OT and OT,c.As a consequence, experimental errors in 0T(Pj) tend to show up more strongly in the HBdistribution function. We have written our own version of HILDA, called A L I N D A (Adamson LINg Distribution Analysis). ALINDA and H I L D A differ only with respect to numerical details (13). Results for Graphon and Sterling M T 2700 as obtained with A L I N D A are shown in Fig. 8. The HB equation is used as local isotherm. Comparing Figs. 7 and 8 shows that the results are very similar, but that A L I N D A tends to transform experimental errors into fine structure o f f ( U ) . It m a y be concluded that CAESAR offers a straightforward computation method for the adsorption energy distribution. In general there is a good agreement between CAESAR, C A E D M O N , and H I L D A (and ALINDA). Journal of Colloid and Interface Science, Vol. 105, No. 1, May 1985

196

VOS AND KOOPAL

a flu)

°o. ~ _ ~ 12,00

14.00--" 16.00 U 18,00

12.00

14.00

%00

b f(U)

16.00

U

18.00

20.00

FIG. 8. Adsorption energy distribution functions of N2 adsorbed on Graphon (a) and Sterling MT-2700 (b), calculated with ALINDA (comparable to HILDA) and using the HB (C = 5.5) equation as local isotherm (compare with Fig. 7).

T h e m a i n a d v a n t a g e o f C A E S A R is t h a t t h e c h o i c e o f r is c o n s i s t e n t w i t h t h e o c c u r r i n g errors, so t h a t f ( U ) is p r e s e n t i n g j u s t t h e right a m o u n t fine structure. I n a s u b s e q u e n t p a p e r it will b e s h o w n t h a t t h e strategy d e v e l o p e d for C A E S A R , c a n also b e a p p l i e d to C A E D M O N a n d to t h e regularization method. The thus evoluated m e t h o d s will b e c o m p a r e d w i t h C A E S A R and ALINDA. ACKNOWLEDGMENTS Discussions with Dr. J. Papenhuijzen and Dr. J. M. H. M. Scheutjens are gratefully acknowledged. REFERENCES 1. Jaroniec, M., Adv. Colloid Interface Sci. 18, 149 (1983).

Journal of Colloid and Interface Science, Vol. 105, No. 1, May 1985

2. House, W. A., "Colloid Science," Specialist Periodical Reports, Vol. 4, Chap. 1. Chem. Soc. London, 1983. 3. Ross, S., and Morrison, I. D., Surf. Sci. 52, 103 (1975). 4. Sacher, R. S., and Morrison, I. D., J. Colloid Interface Sci. 70, 153 (1979). 5. House, W. A., and Jaycock, M. J., J. Chem. Soe. Faraday Trans. 1 73, 942 (1977). 6. House, W. A., and Jaycock, M. J., Colloid Polym. Sci. 256, 52 (1978). 7. House, W. A., J. ColloM Interface Sei. 67, 166 (1978). 8. House, W. A., J. Chem. Soc. Faraday Trans. 1 74, 1045 (1978). 9. Britten, J. A., Travis, B. J., and Brown, L. F., AIChE Symp. Ser. 79, 7 (1984). 10. Brown, L. F., and Tmvis, B. J., in "Fundamentals of Adsorption" (A. L. Myers and G. Belfort, Eds.), Engineering Foundation, New York, 1984. 11. Papenhuijzen, J., and Koopal, L. K., in "Adsorption from Solution" (R. H. Ottewill, C. H. Rochester, and A. L. Smith, Eds.), p. 211. Academic Press, New York]London, 1983. 12. Wesson, S. P., Ph.D. thesis. Rensselaer Polytechnic Institute, Troy, New York, 1975. 13. Koopal, L. K., and Vos, C. H. W., submitted for publication. 14. Young, D. M., and Crowelt, A. D., "Physical Adsorption of Gases." Butterworths, London, 1962. 15. Golub, G. H., and Reinsch, C., in "Handbook for Automatic Computation, Vol. 2: Linear Algebra" (J. H. Wilkinson and C. Reinsch, Eds.), p. 134. Springer Verlag, Berlin 1971; Numer. Math. 14, 403 (1970). 16. Lawson, C. L., and Hanson, R. J., "Solving Least Squares Problems." Prentice-Hall, Englewood Cliffs, N. J., 1974. 17. Hanson, R. J., SIAMJ. Numer. Anal. 8, 616 (1971). 18. Morrison, I. D., and Ross, S., Surf Sci. 39, 21 (1973). 19. Ross, S., and Olivier, J,P., "On Physical Adsorption." Wiley-lnterscience, New York, 1964. 20. Harris, L. B., Surf Sci. 10, 129 (1968). 21. Smith, W. R., "Encyclopedia of Chemical Technology," 2nd ed., Vol. 4, p. 243. Interscience, New York, 1964; Schaeffer, W. D., Smith, W. R., and Polley, M, H., Ind. Eng. Chem. 45, 1721 (1953). 22. Van Dongen, R. H., Ph.D. thesis. Technical University Delft, The Netherlands, 1972. 23. Adamson, A. W., and Ling, I., Adv. Chem. Ser. 33, 51 (1961). 24. Hill, T. L., Z Chem. Phys. 14, 441 (1946). 25. Brunauer, S., Emmett, P. H., and Teller, E., J. Amer. Chem. Soc. 60, 309 (1938).