Langmuir 1997, 13, 1327-1331
1327
Modeling of Porous Media and Surface Structures: Their True Essence as Networks† Vicente Mayagoitia,* Fernando Rojas, Isaac Kornhauser, and Hermilo Pe´rez-Aguilar Departamento de Quı´mica, Universidad Auto´ noma Metropolitana-Iztapalapa, Apartado Postal 55-534, Me´ xico 13, D.F. 09340 Mexico Received September 29, 1995X Networks possess closed paths or cycles connecting their elements. Until now indeed, modeling of porous media and surface structures has been performed disregarding their intricate nature as networks. A qualitative and quantitative adequate simulation of heterogeneous media by Monte Carlo methods is proposed. This allows a deeper understanding of the nature of such media and the prediction of a variety of processes taking part in them.
Introduction If n elements are initially correlated to create a chain, and if now the nth element is connected with the first one to form a cycle, then the previous correlations would lead to a substantial inconsistency, as not only these two new neighboring elements would be scarcely if at all correlated but the whole set of elements would become in a wrong state of interaction. Conceiving cycles having all members properly correlated is the essential problem dealt with in this contribution. Networks surpass other simpler structures (e.g., trees) in that they possess closed paths or cycles connecting their elements. Much evidence (e.g., by optical methods, replicas, etc.) about the morphologies of porous media and surface structures strongly reveals that they rather appear as networks. Although this is not surprising, the amazing point is that, until now, Monte Carlo methods used previously to simulate such structures have been set up disregarding completely the intricacies of this network nature. Fundamentals of the Twofold Description of Porous Media The “twofold description” 1 is particularly suitable to (i) conceive,2,3 represent,4 and classify5 porous network morphologies, (ii) understand and predict the mechanisms of capillary processes5-8 (capillary condensation and evaporation, mercury intrusion, imbibition, immiscible displacement, etc.), (iii) determine the texture of porous † Presented at the Second International Symposium on Effects of Surface Heterogeneity in Adsorption and Catalysis on Solids, held in Poland/Slovakia, September 4-10, 1995. X Abstract published in Advance ACS Abstracts, September 15, 1996.
(1) Mayagoitia, V.; Rojas, F.; Kornhauser, I. Langmuir 1993, 9, 2748. (2) Mayagoitia, V.; Kornhauser, I. In Principles and Applications of Pore Structural Characterization; Haynes, J. M., Rossi-Doria, P., Eds.; Arrowsmith: Bristol, 1985; p 15. (3) Mayagoitia, V.; Cruz, M. J.; Rojas, F. J. Chem. Soc., Faraday Trans. 1 1989, 85, 2071. (4) Cruz, M. J.; Mayagoitia, V.; Rojas, F. J. Chem. Soc., Faraday Trans. 1 1989, 85, 2079. (5) Mayagoitia, V.; Rojas, F.; Kornhauser, I. J. Chem. Soc., Faraday Trans. 1 1988, 84, 785. (6) Mayagoitia, V.; Gilot, V.; Rojas, F.; Kornhauser, I. J. Chem. Soc., Faraday Trans. 1 1988, 84, 801. (7) Mayagoitia, V. In Characterization of Porous Solids II; Rodrı´guezReinoso, F., Rouquerol, J., Sing, K. S. W., Unger, K. K., Eds.; Elsevier: Amsterdam, 1991; p 51. (8) Mayagoitia, V.; Rojas, F.; Kornhauser, I.; Zgrablich, G.; Riccardo, J. L. In Characterization of Porous Solids III (Studies in Surface Science and Catalysis, Vol. 87); Rouquerol, J., Rodrı´guez-Reinoso, F., Sing, K. S. W., Unger, K. K., Eds.; Elsevier: Amsterdam, 1994; p 141.
S0743-7463(95)00812-2 CCC: $14.00
materials,9,10 and (iv) account for the peculiarities of heterogeneous surfaces of adsorbents (adsorption equilibrium and surface diffusion, chemisorption, reaction rate in the adsorbed phase, and surface characterization).11 It also has been applied to the description of other kinds of disordered systems such as arboreous12 and dense13 aggregates resulting from the sol-gel transition, as well as to the assessment of the morphology of polymers.14 This approach could be extended indeed to treat many other different kinds of disordered media. The terminology relevant to porous media is the one employed in this contribution. Nevertheless for surface structures, etc., correspondent parameters can be readily introduced, while mathematical expressions would be practically the same. Every porous network can be visualized as constituted by two kinds of alternated elements: sites (antrae, cavities) and bonds (capillaries, passages). The connectivity, C, is the number of bonds meeting at a site, while every bond is delimited by two sites. For simplicity, the size of each entity is expressed by means of a unique quantity, R, defined in the following way: for sites, considered as hollow spheres, RS is the radius of the sphere, while for bonds, idealized (due to their function of passages) as hollow cylinders open at both ends, RS is the radius of the cylinder. Instead of considering only the size distribution of voids without regard to the kind of element (site or bond) to which they belong, it would seem more appropriate to establish a double distribution of sizes. For this, FS(R) and FB(R) are the normalized size distribution functions, for sites and bonds, on a number of elements basis, so that the probabilities to find a site or a bond having a size R or smaller are, respectively
S(R) )
∫0RFS(R) dR;
B(R) )
∫0RFB(R) dR
(1)
From the very definitions of “site” and “bond” emerges, in a natural way, the following “construction principle”: (9) Mayagoitia, V.; Rojas, F. In Fundamentals of Adsorption III; Mersmann; A. B., Scholl, S. E., Eds.; The Engineering Foundation: New York, 1991; p 563. (10) Mayagoitia, V. Catal. Lett. 1993, 22, 93. (11) Mayagoitia, V.; Rojas, F.; Pereyra, V. D.; Zgrablich; G. Surf. Sci. 1989, 221, 394. Mayagoitia, V.; Rojas, F.; Riccardo, J. L.; Pereyra, V. D.; Zgrablich, G. Phys. Rev. B 1990, 41, 7150. Mayagoitia, V.; Kornhauser, I. In Fundamentals of Adsorption IV; Suzuki, M., Ed.; Kodansha: Tokyo, 1993; p 421. (12) Mayagoitia, V.; Domı´nguez, A.; Rojas, F. J. Non-Cryst. Solids 1992, 147&148, 183. (13) Mayagoitia, V.; Rojas, F.; Kornhauser, I. J. Sol-Gel Sci. Technol. 1994, 2, 259. (14) Kuznetsova, G. B.; Mayagoitia, V.; Kornhauser, I. Int. J. Polym. Mater. 1993, 19, 19.
© 1997 American Chemical Society
1328 Langmuir, Vol. 13, No. 5, 1997
Mayagoitia et al.
the size of any bond is always smaller than the size of the site to which it leads. Two self-consistency laws guarantee the fulfillment of the above principle. The “first law” establishes that bonds must be sufficiently small as to accommodate together with the sites belonging to a given size distribution:
first law
B(R) g S(R)
for every R
(2)
A “second law” is still required since when there is a considerable overlap between the size distributions, there appear topological size correlations between neighboring elements. Thus, the events of finding a size RS for a site and a size RB for a given one of its bonds are not independent. The probability density for this joint event to occur is
F(RS∩RB) ) FS(RS) FB(RB) φ(RS,RB)
(3)
The “second law” can be expressed as:
φ(RS,RB) ) 0
second law
for RS < RB
(4)
If the randomness in the topological assignation of sizes is raised up to a maximum, while complying with the restriction imposed by the construction principle, the “most verisimilar” form of φ for the correct case RS g RB is obtained3
∫
dS exp - S(R ) B B - S φ(RS,RB) ) ) B(RS) - S(RS)
(
S(RS)
patches”, a picture made famous by Ross and Olivier15 in the context of heterogeneous surfaces); a case, more general and realistic perhaps, of an intermediate degree of correlation between element sizes, first considered by Ripa and Zgrablich.16 Naive Methods to Describe Heterogeneous Media. As stated above, in media displaying an intermediate degree of overlap between the site- and bond-size distributions, it is absolutely necessary to observe the “construction principle”. Many ingenuous and very simple methods used by former authors could achieve this purpose. Some of them are: fully random size assignation to sites (that will be termed FRS), followed by size assignation to bonds, as random as possible, while fulfilling the construction principle; a similar procedure in which bonds rather than sites are first imposed (FRB); and size assignation according to an alternated sequence of sites and bonds (ASB). In this case the size of the previous element constrains the election of that of the next one. In all these methods the original size distribution F(R), remains unchanged except that improper sizes are rejected, leading to new effective distributions, F(R)Ψ. In the FRS method, if a bond size is to be chosen, one takes the smallest value, RS, among these of their two delimiting sites, and then Ψ becomes
{
Ψ
)
(
exp -
B(R ) dB ∫B(R ) B - S)
B(RB) - S(RB)
R e RS (6)
for 0
R < RS
while for the FRB method, the corresponding function Ψ for a site having its bond’s greatest value of RB is
S
B
1 B(RS)
(5)
Topological size correlations promote a “size segregation effect”, consisting in that sites and bonds of the bigger size join together to form regions of large elements, while elements of the smaller size reunite to constitute alternated regions of small entities. This effect becomes more important as the overlap increases. Its consequences on the development of capillary processes are of the utmost importance. A pure analytical description of a network over an extended domain of elements would be a mathematical nightmare and this approach is quite impossible for the moment, so one has to resort to Monte Carlo methods, which provide a straightforward solution to this intricate problem. Present Status of Monte Carlo Methods for the Description of Porous Networks Monte Carlo methods have been extensively used to follow the course of capillary processes in porous media as well as surface processes, as adsorption, surface diffusion, and thermal desorption, on heterogeneous surfaces. The main objectives of such studies are to gain an insight about the morphology of the adsorbed phase as well as to obtain quantitative results on static and dynamic behavior. All along these previous works, porous media and surface structures have been conceived in a crude manner. With regard to porous media, three situations have been currently visualized: element sizes are distributed fully at random across the network, that is, size correlations between neighboring elements are absent; correlations are so strong that the network is constituted by an ensemble of homogeneous substructures (“homottatic
{
Ψ
1 1 - S(RB)
R g RB (7)
for 0
R < RB
these procedures being combined in the ASB method. The above observations correspond to a primary recognition of size topological correlations. It is evident from now on that they introduce a distortion in the size distributions of the most constrained elements (bonds in FRS and sites in FRB, in ASB distortion is in some way shared alike by these two kinds of entities). Clearly one observes also that FRS and FRB lead automatically and artificially to networks endowed of a remarkable randomness at the level of a few tens of elements as all members of a certain kind (sites or bonds respectively) display totally random sizes. This extreme randomness is in fact absent in real overlapped structures, which possess, by virtue of the “segregation effect”, alternated domains of considerable extent, constituted by entities of either small or large values of R, grouped together. This morphology has in turn an enormous influence on the structuralization of the diverse phases or fluids during the course of capillary processes.4,17 The “Self-Consistent” (SC) Method.4,17 It is based on the arguments of the twofold description stated in the previous section and allows a better representation of the network, from the point of view of the development of the above mentioned morphology. Now it is useful, for the discussion of examples of networks that will be presented afterward, to treat the (15) Ross, S.; Olivier, J. P. On Physical Adsorption; Interscience: New York, 1964. (16) Ripa, P.; Zgrablich, G. J. Phys. Chem. 1975, 79, 2118. (17) Mayagoitia, V.; Cruz, M. J.; Rojas, F.; Kornhauser, I.; Zgrablich, G.; Pereyra, V. D. Gas Sep. Purif. 1992, 6, 35.
Modeling of Porous Media and Surface Structures
Langmuir, Vol. 13, No. 5, 1997 1329
Figure 1. Precise element-size assignment order schemes for (a) the original self-consistent method (SC) and (b) the improved (ISC) method.
case of the probability density to have a site of size RS having two among its (distinguishable) delimiting bonds of sizes RB1 and RB2, respectively:
F(RS∩RB1∩RB2) ) FS(RS) FB(RB1) φ(RS,RB1) FB(RB2) φ(RS,RB2) (8)
parallel lines are determined to form the surface. Each element gets its size from modified distribution functions having the information about the number and sizes (RB1, RB2 for subsequent sites; RS1 for posterior bonds) of previously, first-order, neighboring settled elements. These modified functions (for the whole of the network) are
Adopting the convention of that RB1 g RB2, and if RS takes any value:
FS(RS) φ(RS,RB1) φ(RS,RB2) FB(RB) φ(RS1,RB)
F(RB1∩RB2) )
∫R
FB(RB1) FB(RB2)
∞ B1
φ(RS,RB1) φ(RS,RB2) FS(RS) dRS (9)
In the event of values RB1 and RB2 being already set, the conditional probability density to find RS is
F(RS/RB1∩RB2) ) F(RS∩RB1∩RB2)/F(RB1∩RB2) ) FS(RS) φ(RS,RB1) φ(RS,RB2)/
∫R∞ φ(RS,RB1) φ(RS,RB2) FS(RS) dRS B1
(10)
Prior to size assignation to elements it is worth making some considerations on the precise way in which these elements are linked together to form the network. For discussion Figure 1a presents a very simple regular 2D network having C ) 4, where only some of the sites of interest are represented. They are termed R, β, γ, δ, ...; the bonds linking them will be named as Rβ, βγ, γδ, .... From Figure 1a it is clear that, e.g., there are two types of γ sites: γ1 possesses only one β-first-order neighboring site, while γ2 has two of them. Thus their branching situation is very different, so that, once a particular site R is distinguished in the network, they are no longer equivalent. The same holds for δ, ‚‚‚, ζ, ‚‚‚. Topological correlations spread their influence hierarchically (from nearer to farther elements). If the size of R is imposed, it influences the four Rβ bonds, then the four β sites, and afterward the βγ, γ, γδ, ‚‚‚ elements. In Figure 1a it is possible to appreciate the influence trajectories that the R site spreads over its neighbors. These sequences are carefully respected in the original SC method when size assignment is undertaken (this is a Markovian process, as the size of an element is constrained only by those of all immediately preceding neighbors): the size of R is imposed randomly. Assignment for elements Rβ, β, βγ1, ... located on a first line is made, and then subsequent
FS(RS) φ(RS,RB1)
for sites
(11)
for bonds sites (first row)
Note that every bond has only one first-order neighboring site which constraints its size. A random number, RND, is generated from a uniform, normalized distribution. This number is taken in correspondence with the pertinent distribution function in eq 11 as follows: when sites on the first row (having only one constraining bond) are being fixed, each site size R is obtained by means of
RND )
∫RR F(RS/RB) dRS ) ∫RR φ(RS,RB1) FS(RS) dRS B1
B1
(12)
In turn for sites of the remaining rows, as every site is delimited by two bonds of values already settled:
∫RR F(RS/RB1∩RB2) dRS ∫R φ(RS,RB1)φ(RS,RB2) FS(RS) dRS ) ∞ ∫R φ(RS,RB1) φ(RS,RB2) FS(RS) dRS
RND )
B1
B1
(13)
B1
RB1 being again bigger than RB2. For bonds, which always are constrained by only one site, their size R is calculated by:
RND )
∫RR
S1
F(RB/RS) dRB )
∫RR
S1
φ(RS1,RB) FB(RB) dRB (14)
The number of elements to be settled must be large enough to reach a satisfactory statistical representation. This depends drastically on the degree of overlap. Shortcoming of the SC Method. Previously, we have drawn attention about systematic deviations, from overlap
1330 Langmuir, Vol. 13, No. 5, 1997
Mayagoitia et al.
Table 1. Distortions of Size Distributions Obtained by the Original SC Method with Respect to Requested Gaussian Values sites
bonds
R h S (Å)
σS (Å)
skewness
requested SC method
600 600.4
50 49.3
0 -0.023
requested SC method
500 501.2
50 48.9
requested SC method
450 466.4
50 44.0
requested SC method
420 506.2
50 39.0
R h B (Å)
σB (Å)
skewness
kurtosis
Overlap ) 0.0455 3.0 2.96
400 400.5
50 49.9
0 0.008
3.0 2.96
0 -0.006
Overlap ) 0.3173 3.0 2.96
400 400.9
50 49.8
0 0.0003
3.0 2.97
0 0.057
Overlap ) 0.6170 3.0 2.98
400 411.7
50 48.4
0 -0.062
3.0 3.01
Overlap ) 0.8414 3.0 3.61
400 480.2
50 43.0
0 -0.300
3.0 3.27
0 -0.44
kurtosis
Table 2. Pore Twofold Size Distributions Obtained by the Improved SC Method against Requested Gaussian Values sites
bonds
R h S (Å)
σS (Å)
skewness
requested ISC method
420 426
50 46.6
0 0.07
requested ISC method
405 403
25 24.6
0 -0.01
kurtosis
R h B (Å)
σB (Å)
Overlap ) 0.8414 3.0 2.94
400 406
50 46.7
Overlap ) 0.9521 3.0 2.28
400 399
25 24.61
skewness
kurtosis
0 0.06
3.0 2.93
0 -0.06
3.0 2.29
Table 3. Associative Energy for Aggregation Distributions Obtained by the Improved SC Method against Requested Gaussian Values sites jS × 10-2/(J mol-1)
bonds
σS
skewness
jB × 10-2/(J mol-1)
kurtosis
σB
skewness
kurtosis
requested ISC method
450 448
50 50.6
0 0.13
Overlap ) 0.6170 3.0 2.98
400 397.6
50 49
0 0.075
3.0 3.06
requested ISC method
420 402.8
50 54.2
0 0.145
Overlap ) 0.8414 3.0 2.86
400 381.5
50 53.9
0 0.149
3.0 2.84
requested ISC method
405 411.7
50 24.2
0 0.308
Overlap ) 0.9521 3.0 1.91
400 406.7
50 24.1
0 0.32
3.0 1.93
values larger than 0.6, of the statistical parameters characterizing sites and bonds. The mean values obtained by the original SC method for both kinds of elements were always bigger than the requested ones. These deviations were more serious as the overlap increased.17 Table 1 shows the point. Overcoming the Problem At that time we were unable to offer any satisfactory explanation about the origin of the above mentioned distortions. It was not until calculations were made by an “improved self-consistent method” (ISC), i.e., after consistent cyclic sequences were performed inside the network following paths as indicated in Figure 1b, that our results dramatically improved, being statistically correct even at very high (=0.95) degrees of overlap. Table 2, depicting the results for a 2D-square lattice of C ) 4 used to visualize the fluid morphology during the course of several capillary processes,8 as well as Table 3, for a 2D-triangular lattice of C ) 6 employed to study the morphology of dense products of aggregation during the sol-gel transition,13 exemplified the excellent agreement that can be reached between calculated and expected values for statistical Gaussian parameters of twofold distributions. Given the limited size of the explored networks (3 × 104 elements), these results are still not strictly representative. However, even if distribution means turn out to be satisfactory, standard deviations σ
are sometimes appreciably lower than the requested values. This is natural and due to the intense correlations between neighboring elements at so extreme overlaps, as a resolute segregation effect renders entities too close in size, within the restricted domain of the limited network of study. This improvement renders possible the description and understanding of the morphology of networks representing diverse heterogeneous media, as well as to assess the behavior of a variety of capillary and surface processes, among others. The Improved Self-Consistent (ISC) Method We are now in position to offer a quite satisfactory explanation for the errors associated with our previous original method: The precise sequence of size assignment to the elements of the network has been described previously, Figure 1a. Now, consider in detail the situation of assigning a size RS to a γ2 site (note that all the sites bearing a subscript 2 or greater, i.e., most of the sites of the network, are in the same situation). In the original method such an assignment is performed using eq 13, which requires that the values of two βγ bonds have been previously assigned. The problem consists precisely in that these bonds have sizes that are not at all in accordance with the crucial fact that they coincide in the same site, γ2. Each of them, in the original method, have been calculated by virtue of a
Modeling of Porous Media and Surface Structures
trajectory of the type Rβγ2, so that their size correlation is merely by means of site R (third order of correlation), given indeed that their correlation is in fact of first order, since they coincide in the same common site γ2. This erroneous procedure for setting βγ bonds in our original method leads to the following: bond sizes are often abnormally different (one is too big and the other too small), contrary to the justified expectation that they should be very close in size, as they reunite at γ2. The effect of this misunderstanding on eq 13 is as follows: given that RS cannot be smaller that the abnormally big βγ bond, the size of γ2 deviates upward. It is interesting to remark that the abnormally small βγ bond has absolutely no influence over this estimation. Then, such error is transmitted and amplified throughout the network. As bond sizes within the network are determined from abnormally big internal sites, the bond-size distribution is also affected. The correct development of Monte Carlo methods for constructing a (i) self-consistent (obeying the construction principle), (ii) verisimilar (displaying the maximum randomness), and (iii) certified (having their elements in correspondence with some previously fixed statistical parameters) porous network, must unavoidably respect the network character of these media: a network possesses closed paths or cycles and, in the case of the models of our interest, the number of cycles practically equals that of sites. Now, in order to certify the correct correlation of all the elements of a given cycle, it is necessary to undertake a precise sequence of element size assignment. For instance, for the first cycle of the network, R β γ2 β′ R shown in Figure 1b, a correct sequence could be in the following order RS(R), RB(Rβ), RS(β), RB(βγ2), RS(γ2), RB(γ2β′), RS(β′), and afterward, as bond β′R must finally arrive to the original site R, this sequence should continue with RB(β′R), RS(R). Nevertheless, the RS(R) size has been already set. Then, a necessary condition for obtaining a consistent cycle is finding, at the end of such cycle, a value identical to that one previously assigned. This is in accordance with a more general principle of consistent networks: each element possesses a unique size. Then, by assigning sizes to elements in a given cycle, values of experiments leading to inconsistent results are rejected, while those allowing cycles to close consistently are accepted. Discussion The present refinement can be inscribed within a series of three major improvements in reference to the original twofold description of networks: (i) the first one was
Langmuir, Vol. 13, No. 5, 1997 1331
related to interactions of bonds meeting at a common site, this effect particularly concerning porous networks,18 (ii) the second one dealt with the possibility of having sites of different connectivities, a fact frequently found in many real structures, the solution of such problem drastically enhancing the generality of our treatment,19 and (iii) now, recognizing the true multicyclical nature of networks, and performing assignations accordingly in closed sequences, the statistical imposed values will be recovered naturally, without resource to any doubtful constraint that could bias the results. Two complementary aspects of correlated networks, the hierarchical spreading of topological correlations, Figure 1a, and the existence of cyclic sequences, Figure 1b, become consistent in our treatment but only if cycles as small as possible are pursued; e.g., a sequence such as R β γ1 δ2 3 δ2′ γ1′ β′ R would be wrong if only because the two δ sites are considered as linked just by 3, while in fact they are also linked by γ2, this one a forgotten entity that will suffer from a critical mistake that has been remarked all along this work. Multiple ways to solve the problem to choose the elements conforming a cycle are perhaps possible. We may have presented probably the crudest one. The most refined could be related to the establishment of the probability density to find the whole cycle in a multivariate sample space. Nevertheless, the one of these that is optimal, from the point of view of speed of computation or any other pragmatical criterion, is totally irrelevant for the moment. Conclusions Twofold description, as stated previously, could only apply to trees, not to networks. But now, the general principles of Monte Carlo simulations for correlated networks have been outlined and evaluated. Acknowledgment. This work was supported and made possible by the National Council of Science and Technology of Me´xico (CONACyT) under Project No. 5387E, as well as under the joint Project “Cata´lisis, Fisicoquı´mica de Superficies a Interfases Gas-So´lido” by CONICET (Argentina) and CONACyT (Me´xico). LA950812M (18) Mayagoitia, V.; Rojas, F.; Kornhauser, I.; Zgrablich, G.; Faccio, R. J.; Gilot, B.; Guiglion, C. Refinementis of the Twofold Description of Porous Media. Langmuir 1996, 12, 211. (19) Mayagoitia, V.; Rojas, F.; Kornhauser, I.; Ancona, E.; Zgrablich, G.; Faccio, R. J. Twofold Description of Topological Disordered Surfaces. Langmuir 1996, 12, 207.
