2964
Ind. Eng. Chem. Res. 1997, 36, 2964-2969
A Rigorous Matrix Approach to Site Percolation for Rectangular Two-Dimensional Grids Jean W. Beeckman*,† Specialty Chemicals Department, Washington Research Center, W. R. Grace & Co.sConn., 7500 Grace Drive, Columbia, Maryland 21044
This paper presents a rigorous approach to calculate detailed site percolation probabilities in a rectangular n × m grid. The probability of finding a site in a particular grid point is allowed to be a function of the grid point location. It is shown that the exact site percolation probabilities can be obtained via a connectivity matrix. Notions of symbolic algebra further allow the site percolation probabilities to be expressed as polynomials, and their coefficients can be deduced exactly. It is shown that when the probability of finding a site in a grid point is independent of the grid point location, then long-range site percolation diminishes exponentially with the length of the grid, while the coefficient in the exponent is an explicit function of an eigenvalue of the connectivity matrix. 1. Introduction Percolation has been a subject of study since it was first introduced in the area of polymer gelation (Flory, 1941; Stockmayer, 1943). A mathematical treatment of the problem was first introduced by Broadbent and Hammersley (1957). An excellent textbook on percolation was written by Stauffer and Aharony (1992). Applications of the concepts of percolation theory to the description of the performance of porous catalysts are many. Amongst those, one finds Beeckman and Froment (1982) on pore blockage in the presence of diffusional limitations in Bethe networks, Marin et al. (1986) on catalyst deactivation by coking during butene dehydrogenation, Beyne and Froment (1993) on zeolite catalyst deactivation by coking, Theodorou and Wei (1983) used Monte Carlo simulations for reaction and diffusion in pore networks, Sahimi and Tsotsis (1985) on deactivation in cubic networks, and Reyes and Scriven (1991) on coke formation with diffusional limitations in a cubic lattice. In the present paper, the problem of site percolation is rigorously dealt with for a simple two-dimensional grid. The n × m grid is characterized by m rows stacked on top of each other, and each row contains n grid points. Each grid point contains a site with probability R or is empty with probability 1 - R, while R is allowed to be a function of the grid point location within the grid. Figure 1 shows a possible configuration for a 10 × 10 grid. The question at hand is to find the probability that sites on the top row are connected through some path with sites on the bottom row (it is assumed that two sites are connected if they are nearest neighbors in the grid). Brute Force Approach. In order to find the number of successful cases for an n × m grid that lead to connectivity between the top and bottom rows, clearly one approach that could be followed would be to actually generate every single arrangement of occupied and unoccupied sites in that grid. Every arrangement could then be tested for connectivity between the top row and the bottom row. This task is enormous even for a 5 × * I dedicate this paper to my mentor, Prof. Gilbert F. Froment, for his 65th birthday and for his continual and inspiring encouragement over the years. † Current address: Mobil Technology Company, P.O. Box 480, Paulsboro, NJ 08066-0480. S0888-5885(96)00663-X CCC: $14.00
Figure 1. Example configuration in a 10 × 10 grid.
5 grid since 25×5 ) 3.36 × 107 different arrangements have to be checked. For a 10 × 10 grid, it becomes impossible since it leads to 210×10 ) 1.26 × 1030 different arrangements. If every arrangement would take 10-4 s to generate and check, it would take on the order of 4.0 × 1018 years to check all of them, a task clearly out of reach. 2. Definitions A rectangular n × m grid is defined as having m rows stacked on top of each other, each row having n grid points. The grid point location in the grid is characterized by a pair (i, j) where i is the location in a row (i ) 1, ..., n) and j is the row number (j ) 1, ..., m). The bottom row is defined by j ) 1, and the top row is defined by j ) m. Every grid point is per definition occupied by a site with a certain probability R(i, j) or empty with a probability 1 - R(i, j). From the above, for any particular grid generated by using a random number generator for instance, certain grid points will be empty and others will have a site located on it according to the local probability R(i, j). The status (occupied or empty) of a group of grid points in a single row is called an assembly, and its occupancy with sites can be represented by a string of zeros and ones. Such a group is denoted by zn(s1,s2,...,sn) with every si (i ) 1, 2, ..., n) being equal to 1 or 0 depending on whether one has a site in location i in the row or not. The total number of possible assemblies in a row is clearly 2n, and this group is denoted by Zn. Figure 2 shows, for instance, a possible assembly for row m with n ) 12. Two sites located in grid points (i, j) and (k, l) are called paired if only one of the indices k or l is different © 1997 American Chemical Society
Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 2965 Table 1. Number of Row Configurations (N(n))
Figure 2. Schematic diagram of an assembly and a row configuration for n ) 12.
from i and, respectively, j by not more than 1. Two sites located in grid points (i, j) and (k, l) are called connected if there is at least one contiguous path that leads from grid point (i, j) to grid point (k, l) via paired sites. Clearly, paired sites are per definition connected. Looking at the top row of an n × m grid, then certain sites may be paired, certain sites may be connected but not paired, and some can be paired and connected. For instance, Figure 2 shows that the sites located in (1,m) and (2,m) are paired and connected, while sites in grid points (4,m) and (6,m) are connected but not paired. Further, sites in grid points (8,m) and (9,m) are paired and connected and are also connected to the sites located in (1,m) and (2,m). We will give any group (one or more) of sites on a single row that are connected an identical numerical value. This value will be numerically equal to 1 if the connection extends through some path all the way to sites (one or more) located on the bottom row. If there is no path extending to the sites on the bottom row, then a unique numerical value greater than 1 will be assigned to that group. According to their connections, the sites on row m can now all be assigned a specific numerical value, and this entire string of numbers is called a configuration. From Figure 2, it can be observed for the assembly as shown on row m and together with the connections shown that a unique configuration is obtained. For the case illustrated, the sites located in grid points (1,m), (2,m), (8,m), and (9,m) are connected but no connection exists to the bottom row; sites in grid points (4,m) and (6,m) are connected but again with no connection to the bottom row; sites in grid points (11,m) and (12,m) are paired and connected to the bottom row; and grid points (3,m), (5,m), (7,m), and (10,m) are empty. The probability of the occurrence of a certain configuration on the top row will be denoted by pn,m(r1,r2,...,rn) and is also the specific site percolation probability for that configuration in the n × m grid. The sum P(n,m) of all pn,m(r1,r2,...,rn) where each configuration has at least one ri (i ) 1, ..., n) equal to 1 is the full probability for bottom row to top row connectivity. P(n,m) is the site percolation probability, specific for the n × m grid. Table 1 shows the number of configurations N(n) possible that have at least one connection to the bottom row, as a function of n. The successive ratios N(n)/N(n1) show that N(n) increases roughly by a factor of 2.75 as n increases by unity. Table 2 shows as an example the configurations that have at least one connection to the bottom row for n ) 2 and 3. Percolation along a String, n ) 1. For n equal to 1, one gets a string of sites that cannot be interrupted in order to have long-range percolation. At the end of the string, we have three configurations possible of which only one is of interest at this point, namely, p1,m(1).
a
na
N(n)
N(n)/N(n-1)
1 2 3 4 5 6 7 8 9 10
1 3 9 25 69 189 518 1422 3915 10813
3.000 3.000 2.778 2.760 2.739 2.741 2.745 2.753 2.762
Number of row grid points.
Table 2. Row Configurations for Grids with Two- and Three-Row Grid Points n ) 2: n ) 3:
(1,1), (0,1), (1,0) (1,1,1), (0,1,1), (1,0,1) (0,0,1), (1,1,0), (0,1,0) (1,0,0), (2,0,1), (1,0,2)
Clearly, one obtains
p1,m(1) ) R(m)p1,m-1(1) and by iteration with p1,1(1) ) R(1) m
p1,m(1) )
R(i) ∏ i)1
Since there are no other configurations leading to longrange percolation, we clearly also have
P(1,m) ) p1,m(1) Percolation on a Ladder, n ) 2. For n equal to 2, we obtain a ladder-type grid. At the top of the ladder, we have three configurations that contribute to longrange percolation, namely, p2,m(1,1), p2,m(0,1), and p2,m(1,0). We will illustrate here the case for R as a function of m only:
p2,m(1,1) ) R(m)2p2,m-1(1,1) + R(m)2p2,m-1(0,1) + R(m)2p2,m-1(1,0) p2,m(0,1) ) R(m)(1 - R(m))p2,m-1(1,1) + R(m)(1 - R(m))p2,m-1(0,1) p2,m(1,0) ) R(m)(1 - R(m))p2,m-1(1,1) + R(m)(1 - R(m))p2,m-1(1,0) or in matrix notation
|p2,m| ) |A(m)| × |p2,m-1| with
| |
p2,m(1,1) |p2,m| ) p2,m(0,1) p2,m(1,0) and
|A(m)| )
|
R(m)2 R(m)2 R(m)2 R(m)(1 - R(m)) R(m)(1 - R(m)) 0 R(m)(1 - R(m)) 0 R(m)(1 - R(m))
|
Matrix A is called the connectivity matrix. By introducing an imaginary zeroth row |p2,0| defined by
2966 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997
||
1 |p2,0| ) 0 0 one obtains by iteration m
|A(i)|) × |p2,0| ∏ i)1
|p2,m| ) (
The column matrix |p2,m| now represents the probability for each configuration on the top row to be connected to at least one site on the bottom row irrespective of the location of that site on the bottom row. The probability for long-range percolation is then obtained from
Figure 3. Schematic diagram of row-to-row connectivity for n ) 10. Table 3. Numerical Values for the Coefficients c(2,m,i) i
P(2,m) ) p2,m(1,1) + p2,m(0,1) + p2,m(1,0)
m
1
2
2
Special Case: r ) Constant. (a) Clearly, P(2,m) will be a polynomial of the following form:
1 2 3 4 5
1 2
i)2m
P(2,m) )
∏
c(2,m,i)Ri(1 - R)2m-i
i)m
where c(2,m,i) is a coefficient that is numerically equal to the total number of configurations that lead to longrange percolation for a ladder grid of length m with exactly i occupied sites. Using notions of symbolic algebra, the coefficients c(2,m,i) can readily be obtained from the above linear equations, and their values are shown in Table 3. The following recursive relationship exists amongst the coefficients:
c(2,m,i) ) c(2,m-1,i-1) + c(2,m-1,i-2) + c(2,m-2,i-3) Since c(2,m,m) ) 2 and c(2,m,2m) ) 1, the above relation allows us to calculate the coefficients effortlessly for any value of m. (b) Calculating the successive values of p2,m(1,1), p2,m(0,1), and p2,m(1,0) for m ) 1, 2, 3, etc., and as a function of R reveals that
limmf∞
p2,m(1,1) p2,m(0,1) ) limmf∞ ) p2,m-1(1,1) p2,m-1(0,1) p2,m(1,0) limmf∞ ) λ(R) p2,m-1(1,0)
Also, numerically, these ratios become constant and equal for relatively small values of m like 5 and upwards. The function λ(R) can easily be obtained from the above by putting
|p2,m| ) λ(R)|p2,m-1| and substituting in the matrix equation. Clearly λ(R) becomes an eigenvalue of matrix A and is obtained as one of the roots of the resulting quadratic.
R λ(R) ) (1 + x1 + 4R(1 - R)) 2 For relatively large values of l and k with k > l, one then obtains
P(2,k) ) P(2,l)e-(k-l)ln(1/λ(R))
3
4
4 2
1 8 2
5
6
7
8
9
10
6 12 2
1 18 16
8 38
1 32
10
1
which clearly shows the exponential decay of the longrange percolation probability for a ladder grid. 3. General Formulation Row-to-Row Connectivity Algorithm. Knowing the configuration for row m of an n × m grid, we will now illustrate how to obtain the configuration of an additional row (m + 1) with a certain assembly as illustrated in Figure 3. In passing the connectivity from row m onto row m + 1, we find that (a) the sites in grid points (3,m+1) and (10,m+1) become connected, (b) the sites in grid points (5,m+1) and (8,m+1) become connected, and (c) the site in grid point (1,m+1) is connected to the bottom row. Updating the connectivity in row m + 1 alone now, we find that (a) the sites in grid points (1,m+1), (2,m+1), and (3,m+1) become connected to the bottom row and (b) since the sites in grid points (3,m+1) and (10,m+1) are connected from the above, also the site in grid point (10,m+1) is connected to the bottom row. Feeding the updated connectivity of row m + 1 back to row m, we find that (a) the site in grid point (8,m) now becomes connected to the bottom row successively via the sites in grid points (10,m+1), (9,m+1), and (8,m+1) and (b) the site in grid point (5,m) now also becomes connected to the bottom row (where it previously was not) via its connection with the site in grid point (8,m). Finally, passing the updated connectivity of row m back onto row m + 1, we find that the site in grid point (6,m+1) is connected to the bottom row via its connection with the sites in grid points (5,m+1) and (5,m). Further connectivity exchange between rows m and m + 1 no longer changes either configuration, which means that the final configuration of row m + 1 has been achieved. It can be observed that all sites in row m + 1 are now actually connected to the bottom row, while before, without row m + 1, only the site in grid point (1,m) was. Also, it is clear that the full connectivity in a row (and not only the connectivity to the bottom row) is required to correctly describe the percolation behavior. Figure 4 describes in general terms the connection algorithm illustrated above. Percolation Probabilities. Consider a rectangular grid which has m rows stacked on top of each other and
Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 2967
where B is a diagonal N(n) × N(n) matrix of the form
|
b1,1 0 b2,2 0 |B| ) · · · 0 0
0 0
· · · · · ·
0 0
0
· · ·
bN(n),N(n)
with the diagonal elements given by
|
n
bj,j ) Figure 4. Connection algorithm.
where each row has n grid points. The top row has a certain number of possible configurations, each one of them describing the connectivity of the sites in the grid points (a) along that row, (b) with the network below it, and (c) connected to the bottom row. The probability of occurrence of each configuration is written as pn,m(r1, r2, ..., rn-1, rn), where pn,m stands for the probability that in a row of n grid points the site in the first grid point has a value r1, the second has a value r2, etc. If, for instance, r2 ) 0, it means that the particular grid point is empty, r2 ) 1 means that the grid point is occupied by a site which is connected to the bottom row, and r2 > 1 means that the grid point is occupied by a site and is potentially connected to other grid points on that row (if other r values exist that are equal to r2) via the network below or along the row itself but not connected to the bottom row. Knowing the configuration pn,m(r1, r2, ..., rn-1, rn) and locating a particular assembly zn(s1, s2, ..., sn-1, sn) for row m + 1 on top contains the required information to determine the new configuration pn,m+1(t1, t2, ..., tn-1, tn) from the row-to-row connectivity algorithm described before. Clearly what is required of zn(s1, s2, ..., sn-1, sn) for any pn,m+1(t1, t2, ..., tn-1, tn) is that sk ) 1 if tk > 0 and sk ) 0 if tk ) 0. The contribution of pn,m(r1, r2, ..., rn-1, rn) to pn,m+1(t1, t2, ..., tn-1, tn) can now be written as
pn,m+1(t1, t2, ..., tn-1, tn) ) pn,m(r1, r2, ..., rn-1, rn) × n
zn
δ(r 98 t)
R(k,m+1) ∏ k)1
sk
1-sk
(1 - R(k,m+1))
zn
The function δ(r 98 t)is equal to 1 if pn,m(r1, r2, ..., rn-1, rn) together with assembly zn(s1, s2, ..., sn-1, sn) leads to pn,m+1(t1, t2, ..., tn-1, tn) and 0 if it does not. The complete probability for pn,m+1(t1, t2, ..., tn-1, tn) can now be written:
pn,m+1(t1, t2, ..., tn-1, tn) )
∑
n
R(k,m+1) ∏ k)1
sk
1-sk
(1 - R(k,m+1))
Writing all pn,m(r1, r2, ..., rn-1, rn) as a column matrix, we get
|pn,m+1(i1, i2, ..., in) ) |A| × |pn,m(i1, i2, ..., in)| Matrix A is called the connectivity matrix and is a product matrix that can be written as follows:
|A| ) |B| × |Ω|
k
where the sk are fully determined by the configuration of row element j in the column matrix |pn,m+1(i1, i2, ..., in)|. The matrix |Ω| is also an N(n) × N(n) matrix with zn
elements δ(r 98 t) which are all equal to 0 or 1 as determined by zn
δ(r 98 t) ) 1 if pn,m(r1, r2, ..., rn-1, rn) yields pn,m+1(t1, t2, ..., tn-1, tn) with assembly zn(s1, s2, ..., sn-1, sn) and zn
δ(r 98 t) ) 0 if it does not. The matrix |Ω| is called the transfer matrix. The total probability that the bottom row is connected to the top row is then given by
P(n,m) )
∑
pn,m(r1, r2, ..., rn-1, rn)
all configs
By introducing the imaginary row |pn,0| of which all elements are put equal to 0 except for element pn,0(1, 1, ..., 1) which is put equal to unity, one obtains in matrix notation
P(n,m) )
{|B| × |Ω|}m × |pn,0| ∑ all elements
P(n,m) expresses the general site percolation probability for having at least one site on the top row connected to at least one site on the bottom row. As an example, consider n ) 3 and R(i,j) ) R ) constant; as given in Chart 1. Then,
|p3,m| ) |B| × |Ω| × |p3,m-1| and
∑
|p3,m|
all elements
all configs
δ(r 98 t)
k
P(3,m) )
pn,m(r1, r2, ..., rn-1, rn) × zn
R(k,m+1)s (1 - R(k,m+1))1-s ∏ k)1
|
The polynomial expression for P(3,m) can quickly be generated for any m by initializing p3,m as follows:
1 0 0 0 |p3,0| ) 0 0 0 0 0
2968 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 Chart 1
|
R3 0 0 0 |B| ) 0 0 0 0 0
|| |
p3,m(1,1,1) p3,m(0,1,1) p3,m(1,0,1) p3,m(0,0,1) |p3,m| ) p3,m(1,1,0) p3,m(0,1,0) p3,m(1,0,0) p3,m(2,0,1) p3,m(1,0,2) 0 R2(1 - R) 0 0 · · ·
0 0 R2(1 - R) 0 0 · · ·
1 1 1 1 |Ω| ) 1 1 1 0 0
· · · 0 0 R(1 -R)2 0 0 · · ·
1 1 0 1 1 1 0 1 0
1 1 1 1 1 0 1 0 0
· · · 0 0 R2(1 - R) 0 0 · · ·
1 1 0 1 0 0 0 1 0
1 1 0 0 1 1 1 0 1
· · · 0 0 R(1 - R)2 0 0 · · ·
1 1 0 0 1 1 0 0 0
1 0 0 0 1 0 1 0 1
· · · 0 0 R(1-R)2 0 0
1 1 0 1 0 0 0 1 0
1 1 0 0 1 0 1 0 1
|
· · · 0 0 R2(1 - R) 0
0 0 0 0 0 0 0 0 R2(1 - R)
|
Table 4. Exact Number S of m × m Grid Configurations That Have Bottom-to-Top Row Connectivitya
a
grid
S
2×2 3×3 4×4 5×5 6×6 7×7 8×8 9×9 10 × 10
7 197 22,193 10,056,959 18,287,614,751 133,267,613,878,665 3,888,492,110,032,895,921 454,016,084,146,597,129,480,639 212,041,997,127,527,144,580,159,558,415
Total possible configurations is 2m×m.
Figure 5. Site percolation polynomials.
For instance,
P(3,3) ) 3R3(1 - R)6 + 22R4(1 - R)5 + 59R5(1 - R)4 + 67R6(1 - R)3 + 36R7(1 - R)2 + 9R8(1 - R) + R9 4. Discussion Figure 5 shows a number of total site percolation probabilities as a function of R (R is assumed not to be a function of the grid point location) for square grids with n ) m and m ranging from 1 to 10. The crossover point on the x axis of successive P(m,m) for increasing values of m approaches, of course, the percolation threshold, Rc. Rc is defined as the value of R where the second derivative of P(m,m) with respect to R is 0 for m infinitely large. Table 4 shows the number of grid configurations with bottom row to top row connectivity as a function of m. This number of successful grid configurations is equal to 2m×m × P(m,m) (with P(m,m) evaluated at R ) 1/2). Figure 6 shows P(5,m) for R independent of the grid point location with m ranging from 1 to 10. As m increases, P(5,m) decreases over the entire R range (except for R equal to unity obviously) since every extra row drops the probability of having the bottom row and top row connected. Figure 7 shows λ(R)/R again for R independent of the grid point location with n ranging from 1 to 9 (and with λ defined in the same way as for
Figure 6. Site percolation polynomial P(5,m).
the case of the ladder grid, n ) 2). It can easily be shown that |d(λ(R)/R)/dR|R)1 ) -1 and is therefore independent of n. This is not surprising since for values of R ever closer to unity, one primarily hass all grid points occupied by sites and long-range connectivity is basically assured. Figure 8 shows the ratio of P(9,m)/P(9,m-1) for R independent of the grid point location and m ranging from 2 to 20. It can be observed that the ratio rather quickly becomes constant for relatively small values of m. Figure 9 shows P(m,m) and a predictor P(m-1,m1)2/P(m-2,m-2). It can be observed, except for m ) 3 at R values above 0.7, that the predictor does an excellent job at forecasting the next higher level P(m,m).
Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 2969
Nomenclature Arabic Symbols A ) connectivity matrix B ) diagonal matrix with diagonal elements b; see text m ) number of rows in the rectangular grid n ) number of columns in the rectangular grid N ) total number of configurations p ) probability of occurrence of a configuration; see text P ) total site percolation probability r ) configuration designator for a grid point s ) assembly designator for a grid point t ) configuration designator for a grid point z ) assembly; see text
Figure 7. λ(R)/R vs R.
Greek Symbols R(i, j) ) probability to have a site in grid point location (i, j) δ ) transfer function between configurations λ ) eigenvalue of matrix A Ω ) transfer matrix with elements δ; see text
Literature Cited
Figure 8. Ratio P(9,m)/P(9,m-1) vs m.
.
Figure 9. P(m,m) and P(m-1,m-1)2/P(m-2,m-2) as predictor for P(m,m).
5. Conclusion This paper elaborates on an approach to calculate exact site percolation probabilities for a general n × m grid. The approach allows for the probability of finding a site in a grid point to be a function of the location of the grid point. The general site percolation polynomial that represents the probability of having the bottom and top row connected is presented in a closed form.
Beeckman, J. W.; Froment, G. F. Deactivation of Catalysts by Coke Formation in the Presence of Internal Diffusional Limitations. Ind. Eng. Chem. 1982, 21, 243-250. Beyne, A. O. E.; Froment, G. F. A Percolation Approach for the Modeling of Deactivation of Zeolite Catalysts by Coke Formation: Diffusional Limitations and Finite Rate of Coke Growth. Chem. Eng. Sci. 1993, 48 (3), 503-511. Broadbent, S. R.; Hammersley, J. M. Percolation Processes. I. Crystal and Mazes. Proc. Camb. Phil. Soc. 1957, 53, 629. Flory, P. J. J. Am. Chem. Soc. 1941, 63, 3083, 3091, 3906. Marin, G. B.; Beeckman, J. W.; Froment, G. F. Rigorous Kinetic Models for Catalyst Deactivation by Coke Deposition: Application to Butene Dehydrogenation J. Catal. 1986, 97, 416-426. Reyes, S. C.; Scriven, L. E. Analysis of Zeolite Catalyst Deactivation during Catalytic Cracking Reactions. Ind. Eng. Chem. Res. 1991, 30, 71. Sahimi, M.; Tsotsis, T. T. A Percolation Method of Catalyst Deactivation by Site Coverage and Pore Blockage. J. Catal. 1985, 96, 552. Stauffer, D.; Aharony, A. Introduction to Percolation Theory, 2nd ed.; Tayler & Francis: London, Washington, DC, 1992. Stockmayer, W. H. J. Chem. Phys. 1943, 11, 45. Theodorou, D.; Wei, J. Diffusion and Reaction in Blocked and High Occupancy Zeolite Catalyst. J. Catal. 1983, 83, 205.
Received for review October 17, 1996 Revised manuscript received January 8, 1997 Accepted January 20, 1997X IE960663I
Acknowledgment I thank Cristian Libanati for helpful discussions and W. R. Grace & Co. for letting me publish this work.
X Abstract published in Advance ACS Abstracts, June 15, 1997.
