- Email: [email protected]
Calculation of the adsorption energy distribution from the adsorption isotherm by singular value decomposition
Colloids and Surfaces, 14 (1985) Elsevier Science Publishers B.V.,
87-95 Amsterdam
-Printed
87
in The Netherlands
CALCULATION OF THE ADSORPTION ENERGY DISTRIBUTION FROM THE ADSORPTION ISOTHERM BY SINGULAR VALUE DECOMPOSITION
LUUK
K. KOOPAL
and KEES
VOS
Laboratory for Physical and Colloid Chemistry, 6703 BC Wageningen (The Netherlands) (Received
6 July
1984;
accepted
in final form
Agricultural 25 September
University,
De Dreijen
6,
1984)
ABSTRACT
An algorithm has been developed to obtain the energy distribution function from an adsorption isotherm. The ill-posedness of the integral adsorption equation is taken into account and negative values in the calculated distribution function are excluded. Computations are performed with computer program CAESAR (computed adsorption energies, SVD analysis result). The method is tested with simulated isotherm data and the results are very promising. The energy distribution functions of four carbon-black samples confirm that graphitization at high temperatures leads to nearly homogeneous surfaces. An analysis of different Aerosils shows that chemical modifications of the surface are reflected in the energy distribution. Comparison of the results obtained with CAESAR with those of other methods indicates that the same trends are observed.
INTRODUCTION
Most powders and precipitates important in industrial processing have non-uniform surfaces. For a quick characterization of these solids, gas-adsorption techniques are frequently used, whereafter surface properties are calculated using equations based on homogeneous (uniform) surfaces. Clearly, the use of such models is contradictory to the heterogeneous character of most surfaces. In the last decade, progress has been made with models which explicitly include surface heterogeneity [l], and numerical computation methods have become available for the analysis of the adsorption isotherm, which give both the surface area and the adsorption-energy distribution function over the sample surface [2] . A problem with these numerical analyses is that good quality data are required and that data manipulation is needed, or that the obtained distribution function depends on an arbitrarily chosen parameter which controls the stability of the calculated distribution function [ 31 .
0166-6622/85/$03.30
0 1985
Elsevier
Science
Publishers
B.V.
88
In this paper, a brief description is given of a new numerical method in which the amount of detail seen in the energy distribution function is linked to the quality of the experimental data and to the precision of the local isotherm function. To illustrate the analysis, use will be made of isotherm data based on known adsorption-energy distribution functions and local isotherm equations. Furthermore, experimental isotherms of N2 adsorption on carbon samples differing in degree of graphitization [4] , and on chemically modified Aerosils [5] will be analyzed. DESCRIPTION
OF THE ANALYSIS
Adsorption on a heterogeneous surface can be envisioned as occurring on patches which individually are energetically homogeneous. Let a fraction of sites having adsorption energies between U and U + dU be f(U) dU and let O(U,P) be the local isotherm function which relates the fraction of sites covered, 0, at equilibrium pressure P and energy of adsorption U. Then the coverage, BT(P), of the total surface is described by the following equation
[21 “max eT(P)
=
s u min
e
(up)f(U)dU
(1)
This is a Fredholm equation of the first kind which is ill-conditioned with respect to the solution of f(U) when BT(P) and 0 (U,P) are known. The first numerical methods used to calculate f(U) have not addressed special attention to this problem, in the regularization method of House due attention is paid [2,3]. However, the latter method allows f(U) to be negative. In the present method both the ill-posedness of (1) and the non-negativity of f(U) are considered explicitly. For numerical evaluation (1) has to be approximated by a summation. The set of linear equations derived from this summation is written in matrix notation eT = Af
(2) where eT = @(Pj); f = f(Ui); A = A(i,j), with A(i,j) the value of AU 6 (Ui,Pj), and NJ is the quadrature coefficient. As a result of the ill-posed character of (1) the matrix A_ is near-singular and cannot be inverted using a standard procedure. To overcome this problem, we applied the singular value decomposition (SVD) method [6] . This method can be used to replace A by a non-singular matrix A’ provided the rank, r, of &is known. Equation (2) with A’ instead of 11_can be used to compute f using a least squares method which does not allow f(Ui) to be negative. We used the NNLS procedure described by Lawson and Hanson [7]. The problem is now reduced to establishing the rank of A. The numerical rank of A can be defined as the number of rows or columns of matrix A
89
which can be considered as linearly independent in the numerical analysis. The rank is calculated as follows. Starting with a lower limit for r the distribution vector -f is calculated and multiplied by the original matrix _A to give the calculated isotherm, 8~,~ eT,c
= 4
(3)
fc
The solution f, is checked solution is obtained when m
BT (pj ) -
eT,c(Pj
)
c
by comparing
@!?Twith i!$,c, and the optimal
=m
(4)
j=l
O.i where Uj is a measure of the error in eT(pj). The value of r is adjusted until 8T,c(pj) is in accordance with (4). This type of calculation method has been proposed by Hanson [8] . In order to perform the calculations a computer program named CAESAR (computed adsorption energies, SVD analysis result) has been written in simula. RESULTS
The method was tested by using isotherm data based on a given distribution function and local isotherm. Results obtained for a bi-Gaussian distribution function f(U) and the Hill-De Boer isotherm as 0 (U,P) are shown in Fig. 1, where the experimental error in eT(p) is varied. It follows that, in the case of an exact local isotherm and errors in f+#) of the order of lo%, valuable information is still obtained on the distribution function with the present method. An accurate determination of f(U) is, however, only possible with accurate data (at low pressures!) as follows from Fig. la. Besides making experimental errors, a wrong choice of the local isotherm function can also be made. For instance, overestimating the lateral inter-
Fig. 1. Effect of a random 0.1%; (b) 10%.
error
in the experimental
isotherm.
Standard
deviation
(a)
90
actions leads to an energy distribution function which is shifted to the lowenergy side and somewhat flattened. In the case of an underestimation of the lateral interactions, the obtained energy distribution is shifted to somewhat larger energy values and the peaks become narrower and higher. In practical situations, the uncertainty in f(U) is not only determined by the error in eT(p), but also by the inaccuracy or the inadequacy of 0 (U,P). It should be realized that all models for adsorption on homogeneous surfaces are only an abstraction of the reality. Therefore, with practical systems the inadequacy of the local isotherm function is estimated by minimizing the difference between !?T and $T,c, irrespective of the shape of the distribution function. In other words, Eqn (2) is solved with respect to f using the NNLS procedure and setting r = n, or giving r a value derived from the machine precision. The obtained value of f is used to calculate aT,c and U(8T,c), the average deviation due to the error in 0 (U,P) (5)
This error final t_hat
procedure only works well with accurate adsorption data, where the made with the local isotherm outweighs the experimental error. In the calculation, f(U) is obtained on the basis of Eqn (4) with the difference oj is now determined by the deviation due to the experimental error in
BT, oj(eT)
and (J@T,d.
Using this procedure, adsorption isotherms [4] of four carbon samples differing in graphitization, viz. Spheron, Spheron-800, Spheron-2700 (Graphon) and Sterling-2700, have been analyzed. From the literature, it is well known that the Spheron and Sterling surfaces become more homogeneous if heated above lOOO”C, due to graphitization. Heating to 2700°C produces a surface which is frequently used as a standard homogeneous surface. Figures 2a, b and c show the result of our analysis for Spheron, Spheron800 and Graphon. The Hill-~-De Boer equation with the lateral interaction parameter equal to 5.5 has been used as the local isotherm [9] . The results indicate that the surface heterogeneity of Spheron has changed upon heating to 800°C but the differences are not spectacular. The differences between Spheron-800 and Graphon are very pronounced, Graphon indeed is much more homogeneous. Figure 2d shows results for Sterling-2700. Comparing Graphon with Sterling-2700 shows that Graphon is slightly less homogeneous. In both cases there are two large and one small peak present. The small peak at the lower energy end may be due to multilayer adsorption. The two large peaks indicate two types of adsorption sites with adsorption energies for N2 differing by about 1 RT. The same experimental data have also been analyzed with CAEDMON [lo] . House and Jaycock [ll] have analyzed Spheron and Spheron-800 using HILDA. The observed trends are the same with all methods. With CAEDMON, a very precise comparison is not possible because in that case
00
U
u,
Fig. 2. Adsorption-energy distribution, f(U), as a function of the adsorption energy, calculated with CAESAR for Spheron (a); Spheron-800 (b); Graphon (c); and Sterling2700 (d) (note the difference in scale of f(U)). The adsorption energies are given in units of RT, with In A’ as reference point. For N, adsorption using the HB-equation, A’ = 7.1 10e9 m2 N-‘.
the virial equation has been used as the local isotherm. HILDA shows more detail, but most of this detail is due to experimental error. The adsorption isotherms of four chemically modified Aerosils [5] have also been used to calculate f(U) using CAESAR. The Aerosil types and their BET and HILDA monolayer capacities for N, [5] are indicated in Table 1. For comparison, the distribution functions have also been calculated with an algorithm which is equivalent to HILDA 12, 111, it is called ALINDA, Adamson-Ling distribution analysis. House has been so kind as to send us a TABLE
1
Characteristics Type
and monolayer
of Aerosil
capacities
(mg g-‘)
of the different
Aerosil
Code
BET
HILDAa
CAESAR
ALINDA
Si-G-Si-GH Si-NH-Si-GH Si+-S i-( CH 3 ) 3 Si--NH-Si-( CH 3 ),
A B C D
44.4 34.6 42.7 38.3
45.5 40.2 48.5 38.6
51.7 39.9 57.8 45.1
49.9 38.8 62.3 47.4
aMultilayer
applied
correction
[ 51.
samples
92
tape of HILDA, but we have f ound it more convenient to write a new program based on the same algorithm. H I L D A and ALINDA differ only with respect t o numerical details and can be seen as 'twin' programs. ALINDA uses a quasi Hermite spline t o represent the integral distribution function instead o f a cubic spline. In t he iterations, the square root of the squared deviations between ~T a n d ~T,c is minimized instead of the rms deviation. The spline subroutines used in ALINDA were taken from IMSL [12]. For the present analysis, the Hill--De Boer isotherm with the lateral interaction parameter equal to 5.5 was used in both programs. Results o f the analyses with CAESAR and ALINDA for the different types o f Aerosil are shown in Figs 3a -- d and Figs 4a - - d, respectively. The m o n o l a y e r capacities obtained using the present methods, d e n o t e d as Vm (CAESAR) and Vm (ALINDA), are summarized in Table 1. It follows t hat the obtained m o n o l a y e r capacities are larger than Vm (BET) and Vm (HILDA). The difference with the latter is probably due to the fact that House et al. [5] applied a multilayer correction based on the BET equation. If, for the m o m e n t , it is assumed that the Vm (BET) values are correct (although the BET model will be in error), t he n the high values obtained with CAESAR
CE ®
CE ®
f(U f(u)
"fro..... 3.oo
10.00
15.00
2o.oo
5. O 0
2 0 •O 0 U
U a
CE ®
CE ©
f(U)
4
f(U
,2 5.00
•
. U
20.O 0
5.
O0
~
20
O0
U
Fig. 3. Adsorption-energy distributions functions calculated with CAESAR (CE) for chemically modified Aerosils. Figures a, b, c and d refer to Aerosil A, B, C and D, respectively.
87-95 Amsterdam
-Printed
87
in The Netherlands
CALCULATION OF THE ADSORPTION ENERGY DISTRIBUTION FROM THE ADSORPTION ISOTHERM BY SINGULAR VALUE DECOMPOSITION
LUUK
K. KOOPAL
and KEES
VOS
Laboratory for Physical and Colloid Chemistry, 6703 BC Wageningen (The Netherlands) (Received
6 July
1984;
accepted
in final form
Agricultural 25 September
University,
De Dreijen
6,
1984)
ABSTRACT
An algorithm has been developed to obtain the energy distribution function from an adsorption isotherm. The ill-posedness of the integral adsorption equation is taken into account and negative values in the calculated distribution function are excluded. Computations are performed with computer program CAESAR (computed adsorption energies, SVD analysis result). The method is tested with simulated isotherm data and the results are very promising. The energy distribution functions of four carbon-black samples confirm that graphitization at high temperatures leads to nearly homogeneous surfaces. An analysis of different Aerosils shows that chemical modifications of the surface are reflected in the energy distribution. Comparison of the results obtained with CAESAR with those of other methods indicates that the same trends are observed.
INTRODUCTION
Most powders and precipitates important in industrial processing have non-uniform surfaces. For a quick characterization of these solids, gas-adsorption techniques are frequently used, whereafter surface properties are calculated using equations based on homogeneous (uniform) surfaces. Clearly, the use of such models is contradictory to the heterogeneous character of most surfaces. In the last decade, progress has been made with models which explicitly include surface heterogeneity [l], and numerical computation methods have become available for the analysis of the adsorption isotherm, which give both the surface area and the adsorption-energy distribution function over the sample surface [2] . A problem with these numerical analyses is that good quality data are required and that data manipulation is needed, or that the obtained distribution function depends on an arbitrarily chosen parameter which controls the stability of the calculated distribution function [ 31 .
0166-6622/85/$03.30
0 1985
Elsevier
Science
Publishers
B.V.
88
In this paper, a brief description is given of a new numerical method in which the amount of detail seen in the energy distribution function is linked to the quality of the experimental data and to the precision of the local isotherm function. To illustrate the analysis, use will be made of isotherm data based on known adsorption-energy distribution functions and local isotherm equations. Furthermore, experimental isotherms of N2 adsorption on carbon samples differing in degree of graphitization [4] , and on chemically modified Aerosils [5] will be analyzed. DESCRIPTION
OF THE ANALYSIS
Adsorption on a heterogeneous surface can be envisioned as occurring on patches which individually are energetically homogeneous. Let a fraction of sites having adsorption energies between U and U + dU be f(U) dU and let O(U,P) be the local isotherm function which relates the fraction of sites covered, 0, at equilibrium pressure P and energy of adsorption U. Then the coverage, BT(P), of the total surface is described by the following equation
[21 “max eT(P)
=
s u min
e
(up)f(U)dU
(1)
This is a Fredholm equation of the first kind which is ill-conditioned with respect to the solution of f(U) when BT(P) and 0 (U,P) are known. The first numerical methods used to calculate f(U) have not addressed special attention to this problem, in the regularization method of House due attention is paid [2,3]. However, the latter method allows f(U) to be negative. In the present method both the ill-posedness of (1) and the non-negativity of f(U) are considered explicitly. For numerical evaluation (1) has to be approximated by a summation. The set of linear equations derived from this summation is written in matrix notation eT = Af
(2) where eT = @(Pj); f = f(Ui); A = A(i,j), with A(i,j) the value of AU 6 (Ui,Pj), and NJ is the quadrature coefficient. As a result of the ill-posed character of (1) the matrix A_ is near-singular and cannot be inverted using a standard procedure. To overcome this problem, we applied the singular value decomposition (SVD) method [6] . This method can be used to replace A by a non-singular matrix A’ provided the rank, r, of &is known. Equation (2) with A’ instead of 11_can be used to compute f using a least squares method which does not allow f(Ui) to be negative. We used the NNLS procedure described by Lawson and Hanson [7]. The problem is now reduced to establishing the rank of A. The numerical rank of A can be defined as the number of rows or columns of matrix A
89
which can be considered as linearly independent in the numerical analysis. The rank is calculated as follows. Starting with a lower limit for r the distribution vector -f is calculated and multiplied by the original matrix _A to give the calculated isotherm, 8~,~ eT,c
= 4
(3)
fc
The solution f, is checked solution is obtained when m
BT (pj ) -
eT,c(Pj
)
c
by comparing
@!?Twith i!$,c, and the optimal
=m
(4)
j=l
O.i where Uj is a measure of the error in eT(pj). The value of r is adjusted until 8T,c(pj) is in accordance with (4). This type of calculation method has been proposed by Hanson [8] . In order to perform the calculations a computer program named CAESAR (computed adsorption energies, SVD analysis result) has been written in simula. RESULTS
The method was tested by using isotherm data based on a given distribution function and local isotherm. Results obtained for a bi-Gaussian distribution function f(U) and the Hill-De Boer isotherm as 0 (U,P) are shown in Fig. 1, where the experimental error in eT(p) is varied. It follows that, in the case of an exact local isotherm and errors in f+#) of the order of lo%, valuable information is still obtained on the distribution function with the present method. An accurate determination of f(U) is, however, only possible with accurate data (at low pressures!) as follows from Fig. la. Besides making experimental errors, a wrong choice of the local isotherm function can also be made. For instance, overestimating the lateral inter-
Fig. 1. Effect of a random 0.1%; (b) 10%.
error
in the experimental
isotherm.
Standard
deviation
(a)
90
actions leads to an energy distribution function which is shifted to the lowenergy side and somewhat flattened. In the case of an underestimation of the lateral interactions, the obtained energy distribution is shifted to somewhat larger energy values and the peaks become narrower and higher. In practical situations, the uncertainty in f(U) is not only determined by the error in eT(p), but also by the inaccuracy or the inadequacy of 0 (U,P). It should be realized that all models for adsorption on homogeneous surfaces are only an abstraction of the reality. Therefore, with practical systems the inadequacy of the local isotherm function is estimated by minimizing the difference between !?T and $T,c, irrespective of the shape of the distribution function. In other words, Eqn (2) is solved with respect to f using the NNLS procedure and setting r = n, or giving r a value derived from the machine precision. The obtained value of f is used to calculate aT,c and U(8T,c), the average deviation due to the error in 0 (U,P) (5)
This error final t_hat
procedure only works well with accurate adsorption data, where the made with the local isotherm outweighs the experimental error. In the calculation, f(U) is obtained on the basis of Eqn (4) with the difference oj is now determined by the deviation due to the experimental error in
BT, oj(eT)
and (J@T,d.
Using this procedure, adsorption isotherms [4] of four carbon samples differing in graphitization, viz. Spheron, Spheron-800, Spheron-2700 (Graphon) and Sterling-2700, have been analyzed. From the literature, it is well known that the Spheron and Sterling surfaces become more homogeneous if heated above lOOO”C, due to graphitization. Heating to 2700°C produces a surface which is frequently used as a standard homogeneous surface. Figures 2a, b and c show the result of our analysis for Spheron, Spheron800 and Graphon. The Hill-~-De Boer equation with the lateral interaction parameter equal to 5.5 has been used as the local isotherm [9] . The results indicate that the surface heterogeneity of Spheron has changed upon heating to 800°C but the differences are not spectacular. The differences between Spheron-800 and Graphon are very pronounced, Graphon indeed is much more homogeneous. Figure 2d shows results for Sterling-2700. Comparing Graphon with Sterling-2700 shows that Graphon is slightly less homogeneous. In both cases there are two large and one small peak present. The small peak at the lower energy end may be due to multilayer adsorption. The two large peaks indicate two types of adsorption sites with adsorption energies for N2 differing by about 1 RT. The same experimental data have also been analyzed with CAEDMON [lo] . House and Jaycock [ll] have analyzed Spheron and Spheron-800 using HILDA. The observed trends are the same with all methods. With CAEDMON, a very precise comparison is not possible because in that case
00
U
u,
Fig. 2. Adsorption-energy distribution, f(U), as a function of the adsorption energy, calculated with CAESAR for Spheron (a); Spheron-800 (b); Graphon (c); and Sterling2700 (d) (note the difference in scale of f(U)). The adsorption energies are given in units of RT, with In A’ as reference point. For N, adsorption using the HB-equation, A’ = 7.1 10e9 m2 N-‘.
the virial equation has been used as the local isotherm. HILDA shows more detail, but most of this detail is due to experimental error. The adsorption isotherms of four chemically modified Aerosils [5] have also been used to calculate f(U) using CAESAR. The Aerosil types and their BET and HILDA monolayer capacities for N, [5] are indicated in Table 1. For comparison, the distribution functions have also been calculated with an algorithm which is equivalent to HILDA 12, 111, it is called ALINDA, Adamson-Ling distribution analysis. House has been so kind as to send us a TABLE
1
Characteristics Type
and monolayer
of Aerosil
capacities
(mg g-‘)
of the different
Aerosil
Code
BET
HILDAa
CAESAR
ALINDA
Si-G-Si-GH Si-NH-Si-GH Si+-S i-( CH 3 ) 3 Si--NH-Si-( CH 3 ),
A B C D
44.4 34.6 42.7 38.3
45.5 40.2 48.5 38.6
51.7 39.9 57.8 45.1
49.9 38.8 62.3 47.4
aMultilayer
applied
correction
[ 51.
samples
92
tape of HILDA, but we have f ound it more convenient to write a new program based on the same algorithm. H I L D A and ALINDA differ only with respect t o numerical details and can be seen as 'twin' programs. ALINDA uses a quasi Hermite spline t o represent the integral distribution function instead o f a cubic spline. In t he iterations, the square root of the squared deviations between ~T a n d ~T,c is minimized instead of the rms deviation. The spline subroutines used in ALINDA were taken from IMSL [12]. For the present analysis, the Hill--De Boer isotherm with the lateral interaction parameter equal to 5.5 was used in both programs. Results o f the analyses with CAESAR and ALINDA for the different types o f Aerosil are shown in Figs 3a -- d and Figs 4a - - d, respectively. The m o n o l a y e r capacities obtained using the present methods, d e n o t e d as Vm (CAESAR) and Vm (ALINDA), are summarized in Table 1. It follows t hat the obtained m o n o l a y e r capacities are larger than Vm (BET) and Vm (HILDA). The difference with the latter is probably due to the fact that House et al. [5] applied a multilayer correction based on the BET equation. If, for the m o m e n t , it is assumed that the Vm (BET) values are correct (although the BET model will be in error), t he n the high values obtained with CAESAR
CE ®
CE ®
f(U f(u)
"fro..... 3.oo
10.00
15.00
2o.oo
5. O 0
2 0 •O 0 U
U a
CE ®
CE ©
f(U)
4
f(U
,2 5.00
•
. U
20.O 0
5.
O0
~
20
O0
U
Fig. 3. Adsorption-energy distributions functions calculated with CAESAR (CE) for chemically modified Aerosils. Figures a, b, c and d refer to Aerosil A, B, C and D, respectively.