Transient Diffusion and Conduction in Heterogeneous Media: Beyond

Jun 15, 1997 - conduction in regular network models of heterogeneous media, including ... cally-random networks, such as the Voronoi network, for whic...
0 downloads 0 Views 285KB Size
Download PDF
Ind. Eng. Chem. Res. 1997, 36, 3043-3052

3043

Transient Diffusion and Conduction in Heterogeneous Media: Beyond the Classical Effective-Medium Approximation Muhammad Sahimi* and Theodore T. Tsotsis Department of Chemical Engineering, University of Southern California, Los Angeles, California 90089-1211

In this paper steady-state and transient diffusion and conduction in heterogeneous media, such as composite solids, are studied. The heterogeneities are of percolation type in the sense that a randomly-selected fraction of the system does not allow transport to take place. We present several new extensions of the classical effective-medium approximation (EMA) for diffusion and conduction in regular network models of heterogeneous media, including an EMA for topologically-random networks, such as the Voronoi network, for which we introduce a novel stochastic Green function, and a cluster (higher order) EMA for transient diffusion and other dynamical properties of heterogeneous media. In all cases the predictions are compared with the results of Monte Carlo simulations. Introduction The problem of transport in heterogeneous media is relevant to a wide variety of phenomena in natural and industrial processes. A partial list of such phenomena includes diffusion in porous media and through biological tissues, flow, dispersion, and displacement processes in rock (for recent reviews, see Sahimi, 1993b, 1995b), transport, mechanical, and rheological properties of heterogeneous materials, such as polymers, glasses, and powders, conduction in composite solids and in amorphous semiconductors, and many others. In the study of transport in such systems, one has to distinguish between macroscopically-homogeneous and -heterogeneous media. In the former class are those media that may have microscopic disorder but on a large enough scale are homogeneous, and therefore their effective transport properties, such as diffusivity and conductivity, are independent of their size or linear dimension. Transport in such systems is described by the classical continuum equations. For example, diffusion and conduction are described by

∂P ) Te∇2P ∂t

(1)

where P is either the mean concentration of the solute or the temperature of the system, at time t, and Te is either the effective diffusivity De or the conductivity Ke. In the second class of heterogeneous media are those that are characterized by large-scale spatial variations of their structural properties. The effective transport properties of such media depend on the length scale of observation or measurements. Such large-scale variations are usually accompanied by long-range correlations. In this paper these types of heterogeneities are not considered; the interested reader is referred to Sahimi (1995a) where this problem was studied. The first step in studying transport in a heterogeneous medium is to develop a model of the system. In this paper, the medium is represented by a network of interconnected bonds and sites. Two- and threedimensional networks are used whose coordination number Z the number of bonds connected to the same site, is either the same everywhere in the network or varies according to a distribution; in the latter case one has a topologically-disordered network. In principle, the bonds represent microscopic (very small) portions of the system (for example, pores of a porous medium) to which effective properties, e.g., conductance or diffusivity, are S0888-5885(96)00602-1 CCC: $14.00

Figure 1. Two-dimensional Swiss Cheese model of disordered continua. Circles represent the inclusions, and the lines indicate the equivalent bonds that represent the transport channels between the inclusions.

assigned. These properties can be selected from appropriate distributions which, in principle, can be measured. In some cases a discrete representation of a heterogeneous medium is not appropriate and/or possible, so that one may have to represent the system by a continuum. For example, granular media composed of random packings of particles, powders, and packed-bed reactors are unconsolidated porous media. It may not be immediately obvious how one represents such media as discrete networks. One can, for example, represent such media by a continumm in which inclusions of a given shape (for example, spherical or ellipsoidal), representing the particles, are inserted in an otherwise homogeneous medium, either at random or in a prescribed fashion. If the inclusions are circular or spherical, then we obtain the so-called Swiss Cheese model (see Figure 1), which has been used extensively in the past. However, if the spheres are not overlapping and if transport takes place through the channels between the spheres, then the problem can be mapped exactly onto the same transport process on the edges of a Voronoi tesselation (Kerstein, 1983), a two-dimensional example of which is shown in Figure 2. On the other hand, if in the Voronoi tesselations one constructs a dual network by connecting the centers of the neighboring © 1997 American Chemical Society

3044 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997

Figure 2. Two-dimensional Voronoi tessellation (light lines) and its dual network (heavy lines) which is a Voronoi network.

polyhedra (polygons) to each other, one obtains a Voronoi network (see Figure 2). Jerauld et al. (1984a,b) have already studied steady-state conduction in such networks. They did not, however, consider unsteady diffusion or conduction in such random structures, a much more complex problem; a goal of this paper is to study this problem. If the inclusions do not have the same size, then van der Marck (1996) has proposed a general method for mapping such a medium onto an equivalent network. Thus, in most cases one can map a continuum onto an equivalent network. One of the simplest prototypes of a disordered medium is a percolation system. In such a system, a fraction of the bonds (or inclusions in the continuum models) do not allow any transport to take place at all. For example, in electrical conduction through a composite solid, a certain volume fraction of the solid may be insulating and not allowing any conduction to take place. In porous catalysts a fraction of the pores may be inaccessible to diffusion of the reactants because such pores are too small to allow the molecules to diffuse through (see, e.g., Deen, 1987; Sahimi, 1992). Thus, the interconnectivity of those regions of the system through which transport does take place and the existence of a sample-spanning cluster of such regions play important roles in the overall behavior of transport processes in such heterogeneous media. Even in a heterogeneous medium which is well-connected, percolation can play an important role if the distribution of the heterogeneities is very broad. For example, if the distribution of the pore permeabilities of a porous medium is very broad, then a finite fraction of the pores has very small permeabilities, and hence their contribution to the overall permeability of the medium is very small and may be neglected. One may set the permeabilities of such pores to be zero and thus reduce the problem to a percolation process. This idea, usually called the critical path analysis, was first developed by Ambegaokar et al. (1971) for determining the hopping conductivity of semiconductors. It was further developed by Katz and Thompson (1986, 1987) for determining the permeability and electrical conductivity of microscopically-disordered but macroscopically-homogeneous porous media saturated by a Newtonian fluid, by Sahimi (1993a) for estimating the effective permeability of a porous medium saturated by a polymer, and by Sahimi and Mukhopadhyay (1996) for determining the effective permeability of a macroscopically-heterogeneous porous medium. Percolation theory has found fruitful applications in many branches of science and technology (see

Sahimi, 1994, for a review). This paper also considers percolation disorder and its effect on the effective transport properties of a heterogeneous medium. Over the past two decades, diffusion and conduction in heterogeneous media have been studied by a variety of methods. Chief among them has been the effectivemedium approximation (EMA), a method by which one replaces a heterogeneous medium by a homogeneous one whose properties mimic those of the heterogeneous system. The EMAs developed so far have been restricted to regular networks. One prime goal of this paper is to extend it to topologically-random networks that are of considerable practical interests. It will be shown that this extension involves the introduction of a novel stochastic Green function which is absent in the treatment of transport in regular networks. Since it is well-known that the predictions of the EMA near the percolation threshold are not very accurate, new cluster EMAs are also introduced for transient diffusion and conduction in heterogeneous media. The plan of this paper is as follows. In the next section the basic concepts, ideas, and results of percolation theory that are used in the rest of the paper will be outlined. Next, diffusion and conduction in heterogeneous media are considered, and an exact but implicit solution of the problem is given. The development of the EMA is the subject of the next section of the paper. The EMA is then extended to topologically-random networks. In the last section a cluster EMA is developed and applied to transient diffusion and conduction in heterogeneous networks. Percolation Concepts Consider a 2D or 3D network of coordination number Z, in which a randomly-selected fraction p of the bonds is open to transport (for example, their conductance is not zero). For convenience they are called the conducting bonds. If p < pc, where pc, is the bond percolation threshold of the network, no sample-spanning cluster of the conducting bonds exists and macroscopic transport does not exist, whereas for p > pc a samplespanning cluster of the conducting bonds is formed and macroscopic transport does take place. The value of pc depends on Z and the dimensionality d of the system. For example, for the square network (Z ) 4), pc ) 1/2, and for the simple-cubic network (Z ) 6), pc = 0.249; in general pc = d/[Z(d - 1)]. One can define a correlation length ξp(p), which is the average distance between two sites belonging to the same cluster. The physical significance of ξp is that for any length scale L . ξp, the system is macroscopically homogeneous. However, for l , L , ξp, where l is the length of a bond, the system is macroscopically heterogeneous, so that transport processes at such length scales cannot be described by the classical continuum equations such as (1). Another important quantity is the accessible fraction of the conducting bonds XA(p), which is the fraction of the conducting bonds that are in the sample-spanning cluster (some of the conducting bonds are not connected to the sample-spanning cluster and are isolated). The sample-spanning cluster can be divided into two parts: the dead-end bonds that carry no flow (of current, heat, or a fluid) and the backbone which is the multiplyconnected part of the cluster that carries all of the flow. Consequently, the backbone fraction XB(p) is defined as the fraction of the conducting bonds that are in the backbone. It should be obvious that XB < XA and that the isolated fraction XI(p) of the conducting bonds is given by XI ) XA - XB.

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3045

The above quantities are geometrical characteristics of the system. The effective transport properties of the system can also be defined in a straightforward way. Suppose that the conductance of the open bonds have been assigned either at random or according to correlated distribution. If a unit potential gradient (for example, temperature, concentration, or voltage gradient) is applied to two opposite faces of the system and if p > pc, then macroscopic conduction or diffusion will take place. The overall effective conductivity Ke(p) of the system obviously depends on p, with Ke(p ) pc) ) 0. In a similar manner, an effective diffusivity De is defined for the system which, as a result of the Einstein relation, is proportional to Ke. One of the most important features of percolating systems is their universal properties near pc. As pc is approached, the correlation length diverges according to the power law

ξp ∼ (p - pc)-ν

(2)

Moreover, near pc one has

Figure 3. Two-dimensional sample-spanning percolation cluster at pc. Table 1. Currently-Accepted Values of the Critical Exponents and the Fractal Dimensions in d Dimensionsa

XA ∼ (p - pc)β ∼ ξp-β/ν

(3)

d

β

ν

βB

Df

DB

µ

2 3

5/36 0.41

4/3 0.88

0.48 0.99

91/48 2.52

1.64 1.87

1.3 2.0

XB ∼ (p - pc)βB ∼ ξp-βB/ν

(4)

Ke ∼ (p - pc)µ ∼ ξ-µ/ν p

(5)

length is short), one must replace ξp in eqs 5 and 6 by L, implying that

De ∼ (p - pc)µ-β ∼ ξ-θ p

(6)

Ke ∼ L-µ/ν

(8)

where θ ) (µ - β)/ν. If there are at most short-range correlations in the system, then the critical exponents ν, β, and βB are universal and do not depend on the microscopic details of the system; they depend only on the dimensionality d of the system. The transport exponent µ is also universal if (Kogut and Straley, 1979; Sahimi et al., 1983a)

De ∼ L-θ

(9)

f-1 )

∫0∞

f(x) dx < ∞ x

(7)

where f(x) is the conductance distribution. Note that although the above scaling laws are supposed to be valid only in the critical region close to pc, in many cases they are actually valid over a much broader region. For any length scale L , ξp, the sample-spanning cluster is a fractal and self-similar object. This means that the cluster looks similar at various length scales; see Figure 3. If Nc is the number of the conducting bonds in the sample-spanning cluster, the fractal dimension Df of the cluster is defined by Nc ∼ LDf and it can be shown that Df ) d - β/ν. Figure 3 shows the sample-spanning cluster at pc. For length scales L , ξp, the backbone is also a fractal and self-similar object whose fractal dimension DB is given by DB ) d - βB/ν and clearly DB < Df < d. Note that at p ) pc, the correlation length ξp is divergent and therefore the sample-spanning cluster and its backbone are both fractal objects at any length scale. The most important implication of fractality of the sample-spanning cluster is that all of its properties at length scales L , ξp are scale-dependent. Since for L , ξp the only relevant length scale of the system is L, if the conductances are distributed randomly (or if the conductance correlation

a

Rational numbers represent exact values.

It will be shown below that such scale-dependent transport coefficients imply diffusion and conduction processes that cannot be described by eq 1. Table 1 provides a list of the current values of the critical exponents and the fractal dimensions. For more discussion of various properties of percolating systems, see Stauffer and Aharony (1992), and for various applications, consult Sahimi (1994). Diffusion and Conduction in Disordered Networks: Exact Formulation Since in a heterogeneous medium the transport coefficients are spatially-varying quantities, eq 1 should be rewritten as ∂P/∂t ) ∇‚(T∇P), where T denotes the spatially-varying transport coefficient. Discretizing this equation by a finite-element or finite-difference method yields

∂Pi ) ∂t

∑j Wij[Pj(t) - Pi(t)]

(10)

where Wij is related to the transport coefficients Ti and Tj at points i and j of the system. For example, for a cubic grid Wij ) (Ti + Tj)/(2δ2), where δ is the distance between i and j. Equation 10 can also be viewed as a probabilistic description of diffusion and conduction in a random network. For example, for diffusion Pi(t) is the probability that the diffusing particle will be found at site i at time t, and Wij is the transition rate between sites i and j (i.e., the probability of diffusing from i to j), taken to be nonzero in this paper only when sites i

3046 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997

and j are nearest neighbors. This probabilistic interpretation indicates how diffusion and conduction in heterogeneous media can be studied by random walk processes (see below). The exact solution of eq 10 is not available for any two- or three-dimensional model of heterogeneous media. Aside from numerical simulations, a powerful method of treating transient diffusion and conduction in heterogeneous media is by perturbation expansion. The general technique for developing a perturbation expansion for such problems has been given by Sahimi et al. (1983a). The method as it applies to diffusion and conduction is summarized below. Taking the Laplace transform of eq 10 yields

λP ˜ i(λ) - δi0 )

∑j Wij[P˜ j(λ) - P˜ i(λ)]

(11)

where λ is the Laplace transform variable conjugate to t and the initial condition Pi(0) ) δi0 has been used. An effective-medium network is now introduced in which all the transition rates are equal to W ˜ e(λ), with site occupation probabilities being P ˜ ei (λ), so that

λP ˜ ei (λ) - δi0 )

∑j W˜ e(λ)[P˜ ej (λ) - P˜ ei (λ)]

(12)

where the effective transition rate W ˜ e is yet to be determined. If one subtracts eq 12 from 11, after some rearrangements, one obtains

˜ i(λ) (Zi + )[P

P ˜ ei (λ)]

-

) ∑j [P˜ j(λ) ˜ i(λ) - P ˜ j(λ)] - ∑ ∆ij[P j P ˜ ej (λ)]

(13)

where Zi is the coordination number of site i,  ) λ/W ˜ e, and ∆ij ) (Wij - W ˜ e)/W ˜ e. A site-site Green function Gij is now introduced by the equation

(Zi + )Gik -

∑j Gjk ) -δik

(14)

∑j ∑k Gij∆jk[P˜ j(λ) - P˜ k(λ)]

(15)

1

d

∞ exp[-x(Z + )/2] ∏ Im (x) dx ∫ 0 2 i)1

i

(16)

(17)

Effective-Medium Approximation The effective transition rate W ˜ e is calculated by the following self-consistency condition: Solve eq 17 for an arbitrary cluster of bonds with transition rates ∆β. Then select any bond R of the cluster and require that

˜ eR 〈Q ˜ R〉 ) Q

(18)

where 〈 〉 denotes an averaging over all possible transition rates of all the bonds in the cluster. Equation 18, when rewritten as 〈Q ˜R - Q ˜ eR〉 ) 0, implies that the fluctuations in the flux in bond R from its value in the effective network must vanish on the average. If one substitutes the effective transition rate so obtained in eq 12 and inverts the Laplace transform, one obtains

∂Pi(t)

∑j ∫0 We(t - τ)[Pj(τ) - Pi(τ)] dτ t

(19)

Equation 19 indicates that matching a heterogeneous medium to a uniform one induces memory effects. The Laplace transform of mean square displacement 〈R ˜ 2(λ)〉 of a diffusing particle is given by

〈R ˜ 2(λ)〉 )

Gij ) G(m1, m2,...,md) )

∑R ∆β γRβQ˜ β(λ)

Similar to eq 15, 17 is also exact and represents the working equation in this paper. It expresses the flux Q ˜ R in the bond R in the heterogeneous medium as the sum of the flux Q ˜ eR in the same bond but in the effective medium, plus the fluctuations in it which are the result of the heterogeneity of the network.

∂t

Equation 15 is an exact but implicit solution of 10 which applies to any network, irrespective of its dimensionality or topological structure. However, for random structures, such as the Voronoi network, further progress requires a statistical treatment of the Green function appropriately coupled to the topological disorder of the network. The statistical treatment of Gij is discussed below. For regular networks, for which the coordination number is constant, the derivation of Gij is relatively simple. For example, for a simple-cubic network in d dimensions (Z ) 2d) if one labels by d integers (m1, m2, ..., md) the relative position of two sites i and j on the network, separated by md bond lengths along the principal axis k, then the Green function Gij is given by

-

Q ˜ R(λ) ) Q ˜ eR(λ) +

)

Gij is the response of the system at j if a unit current is injected into the network at i. In terms of Gij eq 13 is rewritten as

˜ ei (λ) + P ˜ i(λ) ) P

Here Imi(x) denotes the modified Bessel function of order mi. The Green function for many other networks is given by Sahimi et al. (1983a). Equation 15 can also be rewritten in terms of the “flux” Q ˜ ij(λ) ) P ˜ i(λ) - P ˜ j(λ). We denote bonds with Greek letters, assign directions to them, and let γRβ ) (Gil + Gjk) - (Gjl + Gik) be a bondbond Green function, where i and l (j and k) are the network sites with tails (heads) of arrows on bonds R and β, respectively. Similar to Gij, the bond-bond Green function γRβ is the response of the system in bond β if a unit current is injected into bond R of the network. In terms of γRβ one can rewrite eq 15 as

(λZ )W˜ (λ) 2

e

(20)

In general, 〈R2(t)〉 and De are related through the following equation

〈R ˜ 2(t)〉 ) 2dDet

(21)

and therefore eqs 20 and 21 enable one to calculate the effective diffusivity (and thus conductivity) of the medium. Before proceeding further, we must point out that although the quantity We so determined is independent of spatial position, it must contain some time dependence. To see this, we need only consider the most naive approximate analysis of eq 10, in which all of the transition rates Wij are replaced by a single, constant (independent of time) rate Wc. The resulting equation is a discretization of eq 1 for a uniform system and therefore can predict only macroscopic diffusive behavior. Consider now a distribution of transition rates

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3047

which is percolation-like, i.e., only a fraction p < 1 of all bonds of the network are open (have Wij > 0). We know that if p < pc, no connected paths of open bonds span the network, and thus macroscopic diffusion is precluded, whereas the naive replacement of all the transition rates with a constant Wc still predicts macroscopic diffusion, and thus such an approximation fails. Only when the analysis is performed in the Laplace transform space with a time-dependent W ˜ e(λ) are these difficulties automatically circumvented. In practice, however, one cannot solve eq 17 for an arbitrary cluster of bonds, and thus suitable approximate schemes have to be developed. In an EMA approach, one assigns to all but a finite cluster of bonds in the network the effective transition rate W ˜ e(λ) (so that ∆β * 0 only for a finite set of bonds) and proceeds as above, now averaging over the transition rates of the bonds in the cluster to determine W ˜ e(λ). In the simplest case, only a single bond R has transition rate W differing from W ˜ e(λ); this is called single-bond EMA. In this case eqs 17 and 18 yield

〈1 - γ1 ∆ 〉 ) 1 RR

(22)

R

If f(W) is the transition rate distribution, one obtains

f(W)

∫0∞1 - γRR∆R dW ) 1

(23)

We must point out that, if instead of the fluxes Q ˜ R we use the potentials P ˜ i and require that 〈P ˜ i〉 ) P ˜ ei , we would obtain the same equation as 23. Both methods imply that W ˜ e(λ) should be calculated by requiring that e the average excess flux 〈Q ˜R - Q ˜ R 〉 or excess potential e 〈P ˜i - P ˜ i 〉 vanishes. Since percolation disorder gives rise to fractal media, consider a simple case in which

f(W) ) (1 - p)δ+(W) + ph(W)

(24)

i.e., a fraction 1 - p of the bonds do not conduct (do not allow diffusion or conduction to take place) while the transition rates (or, equivalently, the diffusivity or conductivity) of the rest of the bonds are distributed according to h(W), where h(W) is any normalized probability density function. In this case one obtains

p

h(W)

∫0∞ 1 - pc + pcG + pc(1 - G)(W/W˜ e) dW ) p - pc + pcG (25) 1 - pc + pcG

where G() ) -Gii(). Here pc ) 2/Z is the EMA prediction of the bond percolation threshold of a regular network of coordination number Z. In particular, eq 25 predicts that pc ) 1/2 for the square network, an exact result, and pc ) 1/3 for the simple-cubic network, which is about 34% larger than the numerical estimate pc = 0.249 mentioned above. Since at p ) pc the samplespanning cluster of the network is a fractal object, let us briefly consider this case. Setting p ) pc, eq 25 predicts that for a three-dimensional regular network (Sahimi et al., 1983a)

W ˜ e(λ) ∼

[

pcG(0) (1 - pc)h-1

]

1/2

λ1/2

(26)

where h-1 ) ∫∞0 dW h(W)/W is the first inverse moment of h(W). Therefore, if one uses eq 20, one immediately obtains the following equation for large times

〈R2(t)〉 ∼

[

]

pcG(0) 12 t1/2 1/2 (1 - p )h π c -1

(27)

which indicates that the mean square displacement of the diffusants grows slower than it does linearly with t. This unusual phenomenon is called anomalous diffusion (Gefen et al., 1983) or fractal diffusion (Sahimi et al., 1983a). A random walk fractal dimension Dw is now introduced which is defined by

〈R2(t)〉 ∼ t2/Dw

(28)

so that normal or Fickian diffusion corresponds to Dw ) 2. Then, it is obvious that EMA predicts Dw ) 4 for diffusion in percolating networks at pc (or above pc for L , ξp). The accuracy of the result Dw ) 4 can be tested by relating Dw to the critical exponents of percolation defined above. Suppose that diffusion takes place only on the sample-spanning cluster. From eq 6 one has De 2 1/2 , ∼ ξ-θ p , and since fractal diffusion occurs only if 〈R 〉 2 1/2 2 1/2 ξp, one can replace ξp with 〈R 〉 (since for 〈R 〉 , ξp the only relevant length scale of the system is 〈R2〉1/2) and write De ∼ 〈R2〉-θ/2. On the other hand, eq 21 implies that De ∼ d〈R2〉/dt ∼ 〈R2〉-θ/2, which after integration yields 〈R2〉 ∼ t2/(2+θ), and thus

Dw ) 2 + θ ) 2 +

µ-β ν

(29)

This result was first derived by Gefen et al. (1983). Using the numerical values of the exponents given in Table 1, one finds that Dw(d ) 2) ) 2.87 and Dw(d ) 3) = 3.8, indicating fractal diffusion in both two and three dimensions. The three-dimensional value of Dw is only about 5% smaller than the EMA prediction Dw ) 4, and thus this prediction of EMA is reasonably accurate. The crossover between fractal and Fickian diffusion takes w place at a crossover time tco ∼ ξD p , such that for t , tco diffusion is fractal, whereas for t , tco diffusion is Fickian. This implies that for percolation disorder in -3.34, whereas the three dimensions, tco ∼ ξ3.8 p ∼ (p - pc) EMA discussed above predicts that tco ∼ (p - pc)-2. However, both the scaling theory and the EMA predict that as pc is approached, the time scale tco for reaching Fickian diffusion diverges, which is due to the fact that as pc is approached the correlation length ξp also diverges, and thus the heterogeneous medium is a fractal object on a rapidly increasing length scale. This has important implications for interpreting experimental data for the diffusivity, since one has to ensure that a Fickian regime has been reached before assigning constant values to the diffusivity. Equation 27, which predicts a universal value Dw ) 4, is valid as long as h-1, or more generally f-1, is finite. However, Dw will not be universal if h-1 (or f-1) is divergent. For example, if h(W) ) (1 - R)W-R, where 0 < R < 1, EMA predicts that

Dw )

2-R 2(1 - R)

(30)

indicating that, in this case, Dw depends on the details of the distribution function. It should be pointed out that such power-law distributions are not uncommon.

3048 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997

the network with W ) We as the sum of two contributions, a current I0 injected at A and extracted at a very large distance in all directions and an equal current injected at infinity and extracted at B. In both cases, the current flowing through each of the Z equivalent bonds at the point of current entrance is I0/Z, so that a total current of 2I0/Z flows through AB. Thus WAB ) ZWe/2, W′AB ) (Z/2 - 1)We, and therefore

We - W V0 ) Ve W + (Z/2 - 1)We Figure 4. Embedding a bond and its true conductance W in an effective-medium network, in which all the bonds have the same conductance We (left), and its equivalent circle (right). I0 is the injected current, and W′AB is the equivalent conductance of the rest of the network.

In fact, the conductance distribution of the conducting channels in the Swiss Cheese model discussed above has this form (Feng et al., 1987), and therefore such distributions are of practical importance. Effective-Medium Approximations for Topologically-Disordered Networks The predictions of the EMA for steady-state and transient diffusion and conduction in regular networks have been discussed extensively in the past. However, little work has been done for topologically-random structures such as the Voronoi networks. This section considers this problem. Consider first the steady-state limit, λ f 0. The key quantity is γRR which in a sense reflects the topology of the network. Since γRR ) (-2/ Z) + (2G()/Z), in the limit λ f 0 one obtains γRR ) -2/Z, and eq 23 becomes

W -W

e f(W) dW ) 0 ∫0∞ W + (Z/2 - 1)We

(31)

which was first derived by Kirkpatrick (1973). To extend eq 31 to topologically-random networks, we first rederive (31) by an argument from Kirkpatrick (1973). Consider a regular network of coordination number Z, e.g., a simple-cubic network, in which the bond conductances are distributed according to a distribution. Now construct the effective-medium equivalent of the network by replacing the conductances of all the bonds with We. Applying an external potential gradient to this network induces a uniform potential Ve (the analog of Q ˜ eR) along those bonds that are oriented along the external field. Take one of such bonds and replace its conductance We by its true value W in the original heterogeneous network (see Figure 4). This induces an extra current I0 ) Ve(We - W) in the bond, which in turn induces an extra or excess potential V0 along the bond. As discussed above, in the EMA approach one calculates V0 and insists that its average must vanish

∫0∞V0 f(W) dW ) 0

(32)

from which We is determined. To calculate V0 one only needs to calculate W′AB, the conductance of the network between A and B since V0 ) I0/(W + W′AB). But because WAB ) W′AB + We, all one needs to do is calculate WAB, the conductance between A and B in the uniform effective network, which can be determined by a symmetry argument: Express the current distribution in

(33)

Substituting eq 33 into 32 yields 31. Consider now a topologically-disordered network in which each bond has a conductance W, distributed according to f(W). Replace all the conductances by their effective-medium value We. Consider now one bond such that the coordination numbers of its two ends are Z1 and Z2, distributed according to some probability density function, and replace We by the true conductance W of the bond in the original network. If the distribution of the angles between the bonds that are connected to the same site is essentially uniform, as in the case of the Voronoi network, then injecting a current I0 at one end of a bond distributes it essentially as I0/Z1 in each bond connected to that site. Therefore, following Kirkpatrick’s argument, the total current through bond AB is now I0/Z1+I0/Z2 and hence

We - W V0 ) Ve W + [Z1Z2/(Z1 + Z2) - 1]We

(34)

Therefore, if P(Z1, Z2) is the joint probability distribution of Z1 and Z2, then the steady-state single-bond EMA for a topologically-random network, obtained by substituting (34) into (32), is given by

W -W

∫∫∫ W + [Z Z /(Ze 1 2

1 + Z2) - 1]We f(W) P(Z1,Z2) dW dZ1 dZ2 ) 0 (35)

Note that one has to take the average of V0 with respect to both f(W) and P(Z1,Z2). Equation 35 implies that, since Z1 and Z2 are random variables, for topologicallyrandom networks the stochastic bond-bond Green function γRR is given by

γRR ) -

Z1 + Z2 Z1Z2

(36)

which reduces to γRR ) -2/Z given above when Z1 ) Z2 ) Z. Note that Vrettos et al. (1989), who also attempted to derive an EMA for conduction in Voronoi networks, obtained an equation identical to (31), which in our opinion is incorrect. Equations 35 and 36 were tested by numerical simulations. Two-dimensional L × L Voronoi networks were used, where L refers to the square root of the number of Poisson points in the network; L ) 100 was used with 10 different realizations of the network. The bond-bond Green function γRR was then calculated. Since the local environment and topology around each bond of the network are different, one obtains a distribution of the values of γRR, depending on which bond of the network is selected as the reference bond R. The resulting frequency distribution of γRR is shown in Figure 5. It is seen that one can obtain a value of γRR

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3049

Figure 5. Frequency distribution of the bond-bond Green function γRR for the two-dimensional Voronoi network.

anywhere from -0.26 to -0.43, depending on where the reference bond R is located. Most values are, however, concentrated in the relatively narrow range between -0.3 and -0.37. If the average 〈γRR〉 is calculated from its frequency distribution, one obtains 〈γRR〉 ) -0.322. On the other hand, if one uses eq 36 and averages γRR over the distributions of Z1 and Z2, one obtains γRR = -0.33, a difference of about 2% from its value obtained from the frequency distribution of Figure 5, which is presumably due to the finite size of the networks used in the simulations. One may also be tempted to use the naive approximation 〈γRR〉 ) -2/〈Z〉, where 〈Z〉 is the average coordination number of the Voronoi network. Since 〈Z〉 = 6.29 (Jerauld et al., 1984a), one obtains 〈γRR〉 = -0.317, surprisingly close to the numerical estimates. Note that, although in general the coordination numbers Z1 and Z2 of the two sites of any bond are correlated with each other, as Jerauld et al. (1984a) demonstrated, the correlations are weak and thus the two coordination numbers may be assumed to be independent of each other. The numerical accuracy of eq 35 was also tested by considering steady-state conduction in two-dimensional Voronoi networks. A simple case was considered in which each bond of the network had a unit conductance with probability p or zero conductance (i.e., insulating) with probability 1 - p. p was varied, and the effective conductivity Ke ()We) was calculated by imposing on the network a unit potential gradient in one direction and using a periodic boundary condition in the other direction. Ke was calculated by performing two averagings. One was over the spatial distribution of the conducting bonds in a given Voronoi network, while the second averaging was over all the Voronoi networks that one can obtain with a given linear size L (since the initial Poisson points that are used to construct the network are randomly distributed, the topology of the Voronoi networks vary from realization to realization). The results were obtained with 10 different realizations for each of the averaging procedures. Figure 6 presents the results. Equation 35 was also used to calculate Ke ) We using three different approaches. In one the correlations between Z1 and Z2 were taken into account, by the numerical construction of the correlation function for the two coordination numbers of any bond. In the second approach Z1 and Z2 were assumed to be independent of each other, so that P(Z1,Z2) ) ζ(Z1)ζ(Z2) ) ζ(Z)2, where ζ(Z) is the distribution of the local coordination numbers of the

Figure 6. Comparison of the EMA predictions for the effective conductivity Ke of the two-dimensional Voronoi network and its variation with the fraction of the conducting bonds p. The results are for the case in which the two coordination numbers of the reference bond are independently distributed (solid curve) or are correlated with each other (large broken lines). The third curve repesents the case in which an average coordination number has been used in the EMA.

network. In the third approach eq 31 was used with Z = 6.29, the average coordination number of the Voronoi network. As Figure 6 indicates, there are some differences between the three sets of predictions and also with the Monte Carlo results, although the differences are certainly not too large. Also studied were transient diffusion and conduction in the Voronoi networks. The problem was solved by Monte Carlo calculations and also by using the EMA, eq 25. For simplicity, the calculations were performed in the Laplace transform space, and thus W ˜ e(λ) (and De) was calculated for several values of  ) λ/W ˜ e(λ). When the EMA was used, eq 25 was averaged not only over f(W) but also over the distribution of γRR shown in Figure 5. The results are shown in Figure 7, indicating close agreement between the simulation and EMA results, especially for values of  that are not too large (i.e., at relatively long or intermediate times). It should be pointed out that the apparent agreement between the EMA predictions and the Monte Carlo results in two dimensions does not extend all the way to the percolation threshold pc. In fact, very close to pc the difference between the EMA predictions and the simulation results are quite large. The critical exponent µ of the conductivity defined by eq 5 is about 1.3 and 2 (see Table 1) in two and three dimensions, whereas the EMA predicts µ ) 1 in both cases. Near pc the fluctuations in the numerical results are large, and therefore it is difficult to make a meaningful comparison between the EMA predictions and the numerical simulation results. In three dimensions, the EMA predictions for pc are largely in error (see above), and therefore over a relatively large range of pc a comparison between the EMA predictions and the simulation results is meaningless. Cluster EMA for Transient Diffusion and Conduction As discussed above, the predictions of the single-bond EMA are not very accurate near the percolation threshold. The most important reason for the relative failure of the EMA is that, near pc, the currents and potentials

3050 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997

Figure 8. Two- and three-dimensional clusters used in the cluster EMAs. R is the reference bond.

Figure 7. Comparison of the EMA predictions (lines) with the Monte Carlo simulation results (circles) for the effective transition rate W ˜ e(λ) of the two-dimensional Voronoi network. p is the fraction of conducting bonds, and  ) λ/W ˜ e(λ).

in the conducting bonds are highly correlated, so that the state of any bond can affect that of many other bonds in its vicinity. However, in the single-bond EMA such correlations are totally absent, as one embeds only one bond with a random conductance in a uniform (effectivemedium) network and calculates the resulting potential or flux fluctuations in it. Over the past two decades considerable efforts have been devoted to improving the performance of the EMA near pc (see, for example, Sahimi, 1988, and references therein). Most of such efforts have focused on steady-state diffusion and conduction. An obvious method of improving the performance of the single-bond EMA is to use higher order EMAs. As eq 17 indicates, the current Q ˜ R is written in terms of the fluctuations in the conductivities ∆R of the bond R itself and of the other bonds β that surround R. Thus, instead of embedding in the effective-medium only a single bond R with a random conductance, one can embed a cluster of such bonds in the effective-medium network, calculate Q ˜ R for a reference bond R within the cluster, and insist that the average 〈Q ˜ R〉 must be equal to Q ˜ eR, when the averaging is taken with respect to the distribution of the conductivities of all the bonds in the cluster. In practice, one has to apply eq 17 to every bond of the cluster and calculate Q ˜ R numerically. Since near the percolation threshold the currents and potentials in the conducting bonds are correlated, a cluster of bonds includes to some extent the effect of such correlations, and thus one may obtain more accurate predictions. For the steady-state problems this method has already been used by several authors (Blackman, 1976; Turban, 1978; Ahmed and Blackman, 1979; Sheng, 1980; Nagatani, 1981) with various degrees of success. For the transient problem, Haus et al. (1983) employed this method. To obtain rapid convergence with increasing cluster size to the infinite system results, it is important to choose a suitably symmetrical cluster of bonds whose transition rates are to be allowed to be distributed. The cluster EMA must also provide more accurate predictions than the single-bond EMA.

Since the single-bond EMA predicts the exact bond percolation threshold of the square network pc ) 1/2, any cluster EMA must also do so. In three dimensions, since the single-bond EMA predicts pc = 1/3 for the simplecubic network, any cluster EMA must provide a more accurate prediction of pc, closer to the true value pc = 0.249. This is not the case in the work of Haus et al. (1983). These authors used a two-bond cluster EMA which predicted pc = 0.51 and pc = 0.344 for the square and simple-cubic networks, respectively, both of which are less accurate than the predictions of the single-bond EMA. The origin of this discrepancy can be traced back to the fact that a two-bond cluster does not have the symmetry of the network itself. In this paper the smallest clusters that preserve the point symmetry of the networks were used. Since the computations are very invloved and require the calculations of the Green functions Gij and γRβ for many sites and pairs of bonds, the square and simple-cubic networks were used for which eq 16 provides the necessary expression for the Green functions. If one were to use the Voronoi network, one would have to treat all the required Green functions as stochastic variables and determine them numerically, which would necessitate very intensive computations. In two dimensions the cluster consists of nine bonds, whereas in three dimensions it contains fifteen bonds; see Figure 8. The twodimensional cluster reproduces the exact value pc ) 1/2, while the three-dimensional cluster predicts pc = 0.31, an improvement of about 10% over the single-bond EMA. Since the averaging must be taken with respect to all the bonds in the cluster, one must consider all the possible configurations of the cluster as the bond conductivities are varied. Thus, if there are b bonds in the cluster, one has 2b possible configurations (since the conductivity of a bond takes on only 2 values, either 0 or a finite value) which means that for the threedimensional cluster there are 32 768 possible configurations. Use of symmetry can greatly reduce this number, so that one has only 168 distinct configurations in the two-dimensional cluster (instead of 512) and only 2844 in the three-dimensional case. Instead of calculating the effective conductivity or diffusivity, let us study another important related property of heterogeneous media and materials, namely N(ω), the density of states at the frequency ω. N(ω) dω is the number of vibrational modes of a material with a frequency between ω and ω + dω. This is an important quantity for calculating the specific heat and thermal conductivity of heterogeneous materials, such as polymers and other composite solids, and it is also important in the study of proteins and other biomaterials (Elber and Karplus, 1986). The density of states at frequency ω is given by

N(ω) )

1 Im〈P ˜ 0(-ω2)〉 2π

(37)

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3051

Figure 10. Density of states as predicted by the single-bond EMA (solid curves) and cluster EMA (dashed curves) in three dimensions. The results are for p ) 0.316 (A), p ) 0.334 (B), p ) 0.35 (C), and p ) 0.5 (D). Figure 9. Density of states as predicted by the single-bond EMA (solid curves) and cluster EMA (dashed curves) in two dimensions. The results are for p ) 0.5005 (A), p ) 0.505 (B), p ) 0.55 (C), p ) 0.65 (D), and p ) 0.95 (E).

where Im denotes the imaginary part of the complex quantity. Since P0(λ) ) -G(λ)/We(λ), one can calculate 〈P0(λ)〉 by an EMA, make the replacement λ ) -ω2, and thus calculate N(ω). Thus the solution to transient diffusion problem yields directly the density of states. Suppose that the heterogeneities of the system are of percolation type. Then for p > pc, L . ξp, and ω < ωco, where ωco is a critical or crossover frequency, one has N(ω) ∼ ωd-1, where d is the Euclidean dimensionality of the system. However, for ω > ωco one has a different type of scaling for the density of states. In particular, at p ) pc, or for p > pc at length scales L , ξp (i.e., when the sample-spanning cluster is a fractal object) one has (Alexander and Orbach, 1982)

N(ω) ∼ ωDs-1

(38)

where Ds ) 2Df /Dw ) 2(νd - β)/(2ν + µ - β) is called the spectral dimension. Numerically, Ds = 4/3 in both two and three dimensions, and thus N(ω) ∼ ω1/3. It is clear that ωco, which separates the two regimes of the density of states, must be related to the crossover time tco that separates Fickian diffusion from fractal diffuw sion, and in fact 1/ω2co ∼ tco. Since tco ∼ ξD p , one has ωco -Dw/2 νD /2 w ∼ ξp ∼ (p - pc) . Thus, ωco ∼ (p - pc)1.26 in two dimensions and ωco ∼ (p - pc)1.67 in three dimensions, if we use the numerical values of ν and Dw. On the other hand, the EMA predicts that ωco ∼ (p - pc) in both two and three dimension since, as was shown above, the EMA predicts that tco ∼ (p - pc)2. Thus, eq 38 provides another testing ground for the EMA and its refinements. Hence, N(ω) was calculated in both two and three dimensions for several values of p, the fraction of the conducting bonds, in order to assess the accuracy of the single-bond and cluster EMA both near and away from pc. Figures 9 and 10 present the results obtained with both the single-bond and the cluster EMAs. In two dimensions (Figure 9) the difference between the predictions of the single-bond and cluster EMAs becomes smaller as p approaches pc, since both methods predict

the exact value pc ) 1/2, and at p ) pc ) 1/2 the two sets of predictions are virtually identical. The difference between the two EMAs is more pronounced in three dimensions (Figure 10). Since the single-bond EMA predicts pc ) 1/3 for the simple-cubic network, it cannot provide any predictions for p < 1/3, whereas the cluster EMA can. As Figures 9 and 10 indicate, both the singlebond and cluster EMAs predict that for ω > ωco the density of states N(ω) approaches a constant value. That is, both the single-bond and cluster EMAs predict that for ω > ωco one has Ds ) 1 in both two and three dimensions, in qualitative agreement with the prediction (see above) that in both two and three dimensions the spectral dimension has essentially the same value, though the numerical value of the EMA prediction, Ds ) 1, is not in agreement with the true value Ds = 4/3. Therefore, a cluster EMA provides more accurate predictions of the properties of interest, but the convergence of the predictions to the true values is relatively slow. For steady-state diffusion and conduction a combination of the EMA and the real-space renormalization group method (Sahimi et al., 1983b, 1984; Sahimi, 1988; see also Zhang and Seaton, 1992) provides a highly accurate method for estimating the efective conductivity and diffusivity of heterogeneous media even near pc. However, for the corresponding transient problem such a method has not yet been developed. This is due to the fact that it is not yet clear how to develop a real-space renormalization group transformation for transient diffusion and conduction, because it is not obvious how to rescale the time in a renormalization transformation. Thus the problem of accurately predicting transient diffusion and conduction near pc remains unsolved. Summary and Conclusions This paper extended the classical effective-medium approximation to topologically-random networks. It was shown that this entails the introduction of a novel stochastic Green function, which is absent in the EMA for regular networks. Cluster EMAs for transient diffusion and density of states of heterogeneous materials were also presented. In all cases, the predictions were compared with the Monte Carlo simulation results. In two dimensions, the EMA predictions were in close

3052 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997

agreement with the Monte Carlo results, except very close to the percolation threshold. In three dimensions, the cluster EMA provided more accurate predictions than those of the single-bond EMA, although the convergence of the predictions to the Monte Carlo results appears to be slow. Clearly, more accurate analytical approximations for transient diffusion and conduction in heterogeneous media are yet to be developed. Acknowledgment This paper is dedicated to Professor Gilbert F. Froment on the occasion of his 65th birthday. Partial support of this work by the Petroleum Research Fund, administered by the American Chemical Society, is gratefully acknowledged. Literature Cited Ahmed, G.; Blackman, J. A. On Theories of Transport in Disordered Media. J. Phys. C 1979, 12, 837-853. Alexander, S.; Orbach, R. Density of States on Fractals: Fractons. J. Phys. Lett. 1982, 43, L625-631. Ambegaokar, V.; Halperin, B. I.; Langer, S. Hopping Conductivity in Disordered Systems. Phys. Rev. B 1971, 4, 2612-2620. Blackman, J. A. A Theory of Conductivity in Disordered Resistor Networks. J. Phys. C 1976, 9, 2049-2071. Deen, W. N. Hindered Transport of Large Molecules in LiquidFilled Pores. AIChE J. 1987, 33, 1409-1425. Elber, R.; Karplus, M.; Low-Frequency Modes in Proteins: Use of Effective-Medium Approximation to Interpret the Fractal Dimension Observed in Electron-Spin Relaxation Measurements. Phys. Rev. Lett. 1986, 56, 394-397. Feng, S.; Halperin, B. I.; Sen, P. N. Transport Properties of Continuum Systems near the Percolation Threshold. Phys. Rev. B 1987, 35, 197-214. Gefen, Y.; Aharony, A.; Alexander, S. Anomalous Diffusion on Percolating Clusters. Phys. Rev. Lett. 1983, 51, 77-81. Haus, J. W.; Kehr, K. W.; Kitahara, K. Transport in a Disordered Medium: Analysis and Monte Carlo Simulation. Z. Phys. B 1983, 50, 161-169. Jerauld, G. R.; Hatfield, J. C.; Scriven, L. E.; Davis, H. T. Percolation and Conduction on Voronoi and Triangular Networks: A Case Study in Topological Disorder. J. Phys. C 1984a, 17, 1519-1529. Jerauld, G. R.; Scriven, L. E.; Davis, H. T. Percolation and Conduction on the 3D Voronoi and Regular Networks: A Second Case Study in Topological Disorder. J. Phys. C 1984b, 17, 34293439. Katz, A. J.; Thompson, A. H. Quantitative Prediction of Permeability in Porous Rock. Phys. Rev. B 1986, 34, 8179-8181. Katz, A. J.; Thompson, A. H. Prediction of Rock Electrical Conductivity from Mercury Injection Measurements. J. Geophys. Res. 1987, B92, 599-607. Kerstein, A. R. Equivalence of the Void Percolation Problem for Overlapping Spheres and a Network Problem. J. Phys. C 1983, 16, 3071-3075.

Kirkpatrick, S. Percolation and Conduction. Rev. Mod. Phys. 1973, 45, 574-588. Kogut, P. M.; Straley, J. P. Distribution-Induced Non-universality of the Percolation Conductivity Exponents. J. Phys. C 1979, 12, 2151-2159. Nagatani, T. A Two-Bond Theory of Conductivity in BondDisordered Resistor Networks. J. Phys. C 1981, 14, 3383-3391. Sahimi, M. On the Determination of Transport Properties of Disordered Systems. Chem. Eng. Commun. 1988, 64, 179-195. Sahimi, M. Transport of Macromolecules in Porous Media. J. Chem. Phys. 1992, 96, 4718-4728. Sahimi, M. Nonlinear Transport Processes in Disordered Media. AIChE J. 1993a, 39, 369-386. Sahimi, M. Flow Phenomena in Rocks: From Continumm Models to Fractals, Percolation, Cellular Automata, and Simulated Annealing. Rev. Mod. Phys. 1993b, 65, 1393-1534. Sahimi, M. Applications of Percolation Theory; Taylor and Francis: London, 1994. Sahimi, M. Effect of Long-Range Correlations on Transport Phenomena in Disordered Media. AIChE J. 1995a, 41, 229240. Sahimi, M. Flow and Transport in Porous Media and Fractured Rock; VCH: Weinheim, Germany, 1995b. Sahimi, M.; Mukhopadhyay, S. Scaling Properties of a Percolation Model with Long-Range Correlations. Phys. Rev. E 1996, 38703880. Sahimi, M.; Scriven, L. E.; Davis, H. T. On the Improvement of the Effective-Medium Approximation to the Percolation Conductivity Problem. J. Phys. C 1984, 17, 1941-1948. Sahimi, M.; Hughes, B. D.; Scriven, L. E.; Davis, H. T. Stochastic Transport in Disordered Systems. J. Chem. Phys. 1983a, 78, 6849-6870. Sahimi, M.; Hughes, B. D., Scriven, L. E.; Davis, H. T. Real-Space Renormalization and Effective-Medium Approximation to the Percolation Conduction Problem. Phys. Rev. B 1983b, 28, 307311. Sheng, P. Pair-Cluster Theory for the Dielectric Constant of Composite Media. Phys. Rev. B 1980, 22, 6364-6368. Stauffer, D.; Aharony, A. Introduction to Percolation Theory, 2nd ed.; Taylor and Francis: London, 1992. Turban, L. On the Effective-Medium Approximation for BondPercolation Conductivity. J. Phys. C 1978, 11, 449-459. van der Marck, S. C. Network Approach to Void Percolation in a Pack of Unequal Spheres. Phys. Rev. Lett. 1996, 77, 1785-1788. Vrettos, N. A.; Imakoma, H.; Okazaki, M. An Effective Medium Treatment of the Transport Properties of a Voronoi Tesselated Network. J. Appl. Phys. 1989, 66, 2873-2878. Zhang, L.; Seaton, N. A. Prediction of the Effective Diffusivity in Pore Networks near the Percolation Threshold. AIChE J. 1992, 38, 1816-1824.

Received for review October 1, 1996 Revised manuscript received February 3, 1997 Accepted February 7, 1997X IE960602K

X Abstract published in Advance ACS Abstracts, June 15, 1997.