Calculations of the effective diffusion coefcient for

inhomogeneous media
J.R. Kalnin
a
, E.A. Kotomin
a,b,
*
, J. Maier
b
a
Institute of Solid State Physics, 8 Kengaraga str., Riga LV-1063, Latvia
b
Max-Planck-Institut fur Festkorperforschung, Heisenbergstr. 1, D-70569 Stuttgart, Germany
Received 18 November 2000; accepted 31 May 2001
Abstract
We suggest a modication of the well-known Maxwell-Garnett equation for the mobility (effective diffusion coefcient) in
two-phase media (a matrix with inclusions) which permits the description of a wide range of experimental situations. The novel
approach correctly treats the partial trapping of a diffusing particle by an inclusion as well as consequences of an energy barrier
for the particle penetration into an inclusion. Computer simulations show that the presented mean-eld theory reproduces
surprisingly well results for square inclusions without concentration limitation. For inclusions with other shapes (e.g. spherical)
the theory works well up to concentrations at which mobile particles become trapped in `pockets' between inclusions. q 2002
Elsevier Science Ltd. All rights reserved.
1. Introduction
Calculation of transport properties of inhomogeneous
materials has a long history, starting with the pioneering
papers by Maxwell-Garnett [1]. The composite material is
usually modeled as a combination of a host phase (matrix)
characterized by the diffusion coefcient of a probe particle,
D
2
, and spherical inclusions (the second phase) character-
ized by the particle diffusion coefcient D
1
, radius r
0
and
volume fraction F. These two diffusion coefcients are
expressed through the hop length l and the average waiting
time between the two successive hops t:
D
i

l
2
i
2dt
i
; i 1; 2; 1
where d denotes the space dimension (1, 2, or 3).
A very similar problem arises in the description of other
transport coefcients (electrical and thermal conductivity,
dielectric constant, magnetic permeability, elastic moduli,
etc.) in two-phase systems [210]. Examples of systems for
which it is desired to predict such properties are porous
media, polymer blends, foams, and ceramicmetal mixture.
In this paper, we consider cases for which the generally
accepted relation for the effective diffusion fails and
requires generalization.
2. Standard Maxwell-Garnett approach
Let us reproduce briey a typical derivation of what is
generally known as the Maxwell-Garnett (MG) formula.
Experimentally the matrix with inclusions is characterized
by an effective diffusion coefcient, D
eff
, which is a function
of D
1
, D
2
, and F. To determine it, one can use the electro-
static analog. We consider a macroscopically homogeneous
material with the diffusion coefcient D
eff
. Following the
original derivation of MG equation, we imagine that the
particle concentration c has an average gradient ~ g (similarly
to a homogeneous electric eld) along some axis.
Then we insert into the material a spherical inclusion of
radius r
0
surrounded by a spherical shell of a host material
(matrix) with the radius r
1
and assume that the inclusion
does not change the concentration eld outside, i.e. at
r $r
1
. (The radii r
1
and r
0
are dened in such a way that
the inclusion's volume fraction:
F
r
3
0
r
3
1
Journal of Physics and Chemistry of Solids 63 (2002) 449456
0022-3697/02/$ - see front matter q 2002 Elsevier Science Ltd. All rights reserved.
PII: S0022-3697(01)00159-7
www.elsevier.com/locate/jpcs
* Corresponding author at Max-Planck-Institute, Heisenbergstr 1,
70569 Stuttgart, Germany. Tel.: 149-711-6891771; fax: 149-711-
6891722.
E-mail address: [email protected]
(E.A. Kotomin).
is satised) (Fig. 1a). The effective diffusion coefcient D
eff
can be determined from the following steady-state equation:
Dc
r;q
r; q 0; 2
in the coordinates r and q, where q is an angle between ~r
and the external gradient ~ g. The appropriate solution of Eq.
(2) reads:
c
1
r; q Arcosq; 0 , r # r
0
; 3
c
2
r; q Br 1
E
r
2
_ _
cosq; r
0
, r # r
1
4
c
eff
r; q 2grcosq; r
1
, r 5
where c
i
r; q is a local particle concentration in inclusions
(i 1) or in a host material (i 2).
Equations for the unknown constants A, B, E and g arise
from the boundary conditions for the particle concentrations
and uxes:
c
1
r
0
; q c
2
r
0
; q 6
D
1
2c
1
r; q
2r
u
rr
0
D
2
2c
2
r; q
2r
u
rr
0
7
c
2
r
1
; q c
eff
r
1
; q 8
D
2
2c
2
r; q
2r
u
rr
1
D
eff
2c
eff
r; q
2r
u
rr
1
9
From Eqs. (3)(9) we obtain a set of equations:
r
3
0
A 2r
3
0
B 2E 0 10
D
1
r
3
0
A 2D
2
r
3
0
B 12D
2
E 0 11
r
3
1
B 1E 1r
3
1
g 0 12
D
2
r
3
1
B 1D
eff
r
3
1
g 22D
2
E 0 13
Using this set of equations, one gets the MG equation
sought for:
D
eff
D
2
1 1
3D
1
2D
2
F
D
1
12D
2
2D
1
2D
2
F
_ _
14
In general, for an arbitrary space dimension (d 1, 2 and
3) instead of Eq. (14) one nds:
D
eff
D
2
1 1
dD
1
2D
2
F
D
1
1d 21D
2
2D
1
2D
2
F
_ _
15
This result holds not only for a periodic set of the same-
size spherical inclusions but also for a random inclusion
distribution of different r
i
0
if the condition:
F
r
i
0

d
r
i
1

d
(d is space dimension) remains to be fullled (Fig. 1b).
However, the question, at which volume fractions F inclu-
sions begin to `compete' and Eq. (15) is no longer valid,
remains open and could be solved by a comparison with an
analytical theory which takes many-particle effects into
account [8] and/or by means of direct computer simulations.
3. Disadvantages of MG equation
Consider now several cases when the generally accepted
MG Eq. (15) gives incorrect results.
1. Let us begin with a situation when the inclusion is totally
impenetrable, i.e. a diffusing particle is reected at r
0
as
could be described by the condition:
2c
2
r; q
dr
u
rr
0
0
(consider also the limiting case
D
1
D
2
! 0
in Eq. (7)). In this situation Eq. (15) gives:
D
eff
D
2
1 2
2F
1 1F
_ _
; d 2 16
J.R. Kalnin et al. / Journal of Physics and Chemistry of Solids 63 (2002) 449456 450
Fig. 1. (a) Schematic presentation of the matrix with inclusions
in terms of a core-shell model; (b) the case of different-size
inclusions.
D
eff
D
2
1 2
3F
2 1F
_ _
: d 3 17
Such relations are well-known in the reaction-rate theory
[12]. The same result may be obtained from Eqs. (39),
putting there c
1
r; q 0 (as well as A 0, see discus-
sion in Ref. [11]).
However, Eqs. (16) and (17) give an incorrect concentra-
tion dependence (see discussion in Ref. [7]). The correct
F-dependence, as we show below, is:
D
eff

D
2
1 2F
1 2
3F
2 1F
_ _
; d 3 18
D
eff

D
2
1 2F
1 2
2F
1 1F
_ _
: d 2 19
The reason for this incorrectness lies in the use of relation
(Eq. 8). In fact, the concentration of diffusing particles in
the matrix cannot be equal to that in the effective medium
because in the latter all particles are stirred over a whole
system's volume and thus their averaged concentration
should be less by the factor of 1 2F. This indicates that
Eq. (8) should be corrected as:
c
2
r
1
; q k
1
c
eff
r
1
; q; 20
i.e. in reality there is a jump in concentration on the core
(inclusion)-shell (matrix) boundary, r r
1
. The question
is, how to get the coefcient k
1
? We propose to obtain it,
as a much better approximation, from the expression for
the average equilibrium particle concentration of the
system (total particle number divided by total volume):
c
eff
c
1
F1c
2
1 2F: 21
In the particular case of impenetrable inclusions (c
1
0)
we obtain from Eq. (21):
c
eff
c
2
1 2F; k
1

1
1 2F
: 22
2. Another restriction of the use Eq. (15) arises from the
Maxwell's boundary condition c
1
r
0
; q c
2
r
0
; q, Eq.
(6). In fact it can be shown that Eq. (15) is valid in the
case of different diffusion coefcients in the matrix and
the inclusion, D
1
D
2
, only if the particle velocities in
the matrix and inclusions coincide,
l
1
t
1

l
2
t
2
:
In a general case it is necessary to introduce the measur-
able coefcient k connecting c
1
and c
2
:
c
1
r
0
; q kc
2
r
0
; q: 23
In order to improve the MG equation, in the above-
presented standard derivation of the effective diffusion
coefcient, we can use Eqs. (20) and (23) instead of the
original Maxwell's Eqs. (6) and (8). When doing so, instead
of standard Eq. (15) we arrive at:
D
eff
D
2
k
1
1 1
dD
1
k 2D
2
F
kD
1
1d 21D
2
2kD
1
2D
2
F
_ _
24
In the 1-D case this equation coincides with the exact
solution [13].
The two coefcients k and k
1
are related through the equi-
librium concentrations in inclusions and the matrix and
volume fraction F according to Eqs. (20), (21), (23):
k
c
1
c
2
25
and:
c
2
k
1
c
eff
: 26
From Eq. (21) we obtain the volume dependence of the
coefcient k
1
:
k
1

1
1 2F1
c
1
c
2
F
: 27
Using Eq. (27), one obtains the following relation instead
of Eq. (24):
D
eff

D
2
1 2F1
c
1
c
2
F
1 1
d D
1
c
1
c
2
2D
2
_ _
F
d 21D
2
1
c
1
c
2
D
1
2
c
1
c
2
D
1
2D
2
_ _
F

_
_
_

_
_
_
:
28
Similarly to the MG theory, this equation reproduces
correctly both limiting cases, as F strives for zero and
unity. Eq. (28) is a basic result of our theory. It should be
recalled that c
1
and c
2
are average concentrations of diffus-
ing particles in the two phasesthe inclusions and the
matrix.
It is convenient to express the ratio c
1
/c
2
entering Eq. (28)
through the kinetic parameters of inclusions and the matrix.
In equilibrium, the steady-state situation uxes of particles
to and from inclusions are equal:
c
1
l
1
t
1
c
2
l
2
t
2
: 29
Remember that the diffusion coefcients D
1
and D
2
are
dened by Eq. (1).
Fig. 2 shows several important situations for the potential
energy proles of the diffusing particle, modeling its partial
trapping by an inclusion (potential energy well) and the
(partial) reection from it due to the energy barrier, respec-
tively. To describe these situations, let us introduce the
penetration probabilities p
1
from the inclusion to the matrix
and p
2
from the matrix to inclusions, respectively. Thus, in
the general case one gets:
c
1
c
2

l
2
p
2
t
1
l
1
p
1
t
2
: 30
J.R. Kalnin et al. / Journal of Physics and Chemistry of Solids 63 (2002) 449456 451
In the case of a potential barrier the penetration probabil-
ity (per unit time) is dened entirely by the activation energy
E
12
[14]:
p
2
t
2
z exp
2E
12
kT
_ _
: 31
The same is true for the particle hop from the inclusion
with the probability p
1
. Eqs. (28) and (30) allow one to
describe many diffusion-controlled processes in composite
media with trapping and release of mobile particles.
Now let us compare our results with previous theories. In
the 1-D case Eq. (28) reproduces the exact result derived for
a periodical inclusion distribution in the Kronig-Penny
model with particle reections (Eq. (4) in Ref. [13]) which
reads in our notations as:
D
eff

F
D
1
k
1
1 2F
D
2
_ _
21
z
1
1 2F1Fk
: 32
where k c
1
/c
2
.
In turn, in the 3-D case our Eq. (28) for D
eff
coincides with
that obtained in Ref. [7] for a particular case of periodic
distribution of inclusions using mathematically very com-
plicated formalism of irreversible thermodynamics with
chemical potentials. The analytical results presented in
Ref. [15] also demonstrate the presence of a distinctive
co-factor:
1 2F1
c
1
c
2
F
_ _
21
entering the D
eff
.
In the case of a complete particle reection from the
inclusions; D
1
0;
c
1
c
2
0, Eq. (28) transforms into Eqs.
(18) and (19) as quoted above. For a small volume fraction
of the inclusions, F !1, one arrives at:
D
eff
D
2
1 2
F
2
_ _
; d 3 33
D
eff
D
2
1 2F: d 2: 34
Eq. (33) was received earlier in Ref. [16] whereas
Eq. (34) was derived in Ref. [11] using the effective medium
approximation. The expression (Eq. 33) was also derived
calculating the effective self-diffusion constant of the
mobile species in solution [17].
Let us consider now several limiting cases for D
eff
. In the
case of impenetrable inclusions (complete reection of
particles, c
1
/c
2
! 0), Eq. (28) is simplied:
D
eff

D
2
d 21
d 21 1F
: 35
As F ! 1, one gets in 2-D and 3-D cases D
eff
D
2
/2 and
2D
2
/3, respectively, while the MG equation yields a zero
effective diffusion coefcient. This result means that in
these limiting cases the particle diffusion becomes in fact
one- or two-dimensional as shown in Fig. 3 for square and
cubic inclusions.
To stress the role of a potential barrier in D
eff
, let us
consider the limiting case of a strong trapping, c
2
/c
1
! 0.
From Eq. (28) one obtains:
D
eff

D
2
kF
z
1 1d 2F
1 2F
; 36
which demonstrates that D
eff
is independent of the diffusion
in inclusions. The same results in the limiting case D
2
/
D
1
! 0.
J.R. Kalnin et al. / Journal of Physics and Chemistry of Solids 63 (2002) 449456 452
Fig. 2. Different cases of energy barriers between the matrix and
inclusions, l
i
, E
i
are a hop length and an activation energy for diffu-
sion in the two phases, i 1, 2: (a) an inclusion with the diffusion
coefcient in the inclusions smaller than in the matrix; there is no
energy barrier between them; (b) an energy barrier E
a
for the pene-
tration into inclusion, p
2
#1; (c) partial trapping of particles inside
inclusions, p
1
#1. The detrapping energy is E
t
.
4. Random walk simulations
In order to compare Eq. (28) with computer simulations,
we have modeled periodical arrays of spherical and square
inclusions of the same size in 2-D varying the kinetic para-
meters l
1
, t
1
, and l
2
, t
2
for particle diffusion in the matrix
and inclusions in a very wide range of magnitudes thus
simulating very different situations mentioned above. For
a periodical distribution of inclusions in the matrix we moni-
tored a particle diffusion and calculated D
eff
by the standard
formula:
D
eff

, r
2
.
2dt
; 37
where t is diffusion time and the mean-square particle
displacement during its random walks on the lattice:
, r
2
.

N
i1
r
2
i
N
38
was averaged over more than (typically) N 10
5
10
6
runs.
The waiting time t was chosen to be sufcient to satisfy the
standard condition: ,r
2
. is much larger than the squared
distance between adjacent inclusions L
2
. For this purpose we
used the rst-passage algorithm [8]. We modeled cases of
both impenetrable and penetrable inclusions. The hop length
l was always chosen to be much smaller than both the short-
est distance between boundaries of the two nearest inclu-
sions and inclusion radius. Results of computer simulations
are discussed below.
5. Simple exactly solvable model
There is one particular case in which D
eff
could be deter-
mined exactly for the two-phase inhomogeneous media in
all dimensions: l
1
l
2
, t
1
t
2
, D
1
D
2
. That is, the waiting
times in the matrix and inclusions differ but hopping lengths
are equal (Fig. 2a). In this case after N walks we get from
Eq. (37):
D
eff

, r
2
.
2dN
1
t
1
1N
2
t
2

39
Here N
1
and N
2
are numbers of particle walks in the
phases 1 and 2, respectively. For sufciently large N (diffus-
ing particle visits inclusions many times) one obtains,
obviously:
N
1
FN; N
2
1 2FN 40
Substituting Eq. (40) and ,r
2
. 2dNl
2
into Eq. (39),
we receive results well-known for conductivity in inhomo-
geneous media [5,6]:
1
D
eff

1
D
1
F1
1
D
2
1 2F 41
Note that this equation is often considered to be valid only
for 1-D but as we have demonstrated, in fact it could be used
in any space dimension. (Compare this equation with
Eq. (32) where energy barriers for particle penetration to/
from inclusions are incorporated.)
The same result also follows immediately from our
general Eq. (28), taking into account that at l
1
l
2
:
c
1
c
2

t
1
t
2
42
which leads to
D
eff

D
2
1 2F1
t
1
t
2
F
: 43
J.R. Kalnin et al. / Journal of Physics and Chemistry of Solids 63 (2002) 449456 453
Fig. 3. Transformation of the 2-D (3-D) diffusion into 1-D and 2-D
diffusion as square (cubic) inclusion fractional volume approaches
unit ((a) and (b), respectively).
Eq. (43) coincides with Eq. (41), we just have to
replace: t
1
=t
2
by D
2
=D
1
:
The effective diffusion coefcient could be easily related
to the fractions of time which mobile particle spends in
matrix (t
1
) and inclusions (t
1
):
D
eff
D
1
t
1
t
1
1t
2
1D
2
t
2
t
1
1t
2
44
From Eqs. (28) and (44) one obtains:
t
2
t
1
1 2
1
A
; A

k1 2d 2F
1 2F1kF k
D
1
D
2
1 2F 1d 21 1F
_ _
45
6. Results of computer simulations
It should be mentioned that at high concentrations of the
inclusions their shape becomes important (Fig. 4)the
effect never discussed earlier in the literature. If inclusions
are circular or spherical and touch each other at one point,
D
eff
! 0 as F ! F
lim
(which could be easily calculated as
F
lim
p/4 and 3/p for circular inclusions in the square
lattice and hexagonal lattice, respectively). The reason is
that a particle spends most of its time in a `pocket' formed
by nearest inclusions (Fig. 5).
Let us discuss here results of modeling inclusions of
different shapes as shown in Fig. 4. (Very preliminary
results were presented in Ref. [18]). First of all, our
computer simulations clearly demonstrated the correct-
ness of our theory in the case of large concentrations of
the square inclusions (Fig. 6) in a whole range of
concentrations and failure of the MG theory. Our simu-
lations agree within 1% with exact but complicated
result [19]. This gure presents also the results of cal-
culations for the effective diffusion coefcient at com-
pletely reecting circular inclusions in the 2-D square
and hexagonal lattices (Fig. 4a,c), D
1
=D
2
0:
The immediate conclusion can be drawn that the com-
puter simulations also coincide, with a precision of 1%,
with our theory, Eqs. (28) and (35), for the case of spherical
inclusions, up to inclusion volume fractions as large as
F 0.6 and 0.8, for square and hexagonal lattices, respec-
tively. The discrepancy at larger volume fraction is due to
J.R. Kalnin et al. / Journal of Physics and Chemistry of Solids 63 (2002) 449456 454
Fig. 4. Different types of lattices and inclusions: spherical and
square shape inclusions in the square lattice (a, b), spherical inclu-
sions in the hexagonal lattice (c).
Fig. 5. A pocket effect when particles become trapped between
inclusions which results in zero effective diffusion coefcient.
the just explained `pocket' effect (Fig. 5) neglected in the
analytical theory. What should be stressed here is that
the MG Eq. (14) (dotted line) gives a rather incorrect F-
dependence even at small F since it neglects differences
in the particle concentrations in inclusions and in matrix.
This becomes very transparent in the case of impenetrable
inclusions.
7. Discussion and conclusions
In this paper a modied Maxwell-Garnett Eq. (28) was
derived which gives surprisingly good description of the
mobility in the 2-D heterogeneous media with square inclu-
sions. (We believe that this is also true for 3-D case.) The
main advantage of our mean-eld theory is its simplicity,
transparency and validity for any space dimensionin
contrast to previous mathematically very complicated
approaches based on irreversible thermodynamics with
chemical potentials [7], cluster expansion [15], etc. Validity
of our universal relation for the effective diffusion coef-
cient is proved by computer simulations. These latter
demonstrated for the rst time the effect of inclusion
shape at high inclusion concentrations.
Note that the concentration jump on the inclusion bound-
ary, Eq. (25), commonly known as partition coefcient, is
not new and was observed in all previous theories
[15,20,21]. The more so, this partition coefcient was
found also experimentally, e.g. for metallic alloys Al-Cu,
Al-Ge, Al-In, Al-Sn [20]. The concentration jump is not
surprising which becomes obvious in the limiting case of
impenetrable inclusions (c
2
.0, c
1
0).
The expression suggested for the effective diffusion
coefcient permits treatment of solid-state inhomogeneous
systems (composites, ceramics) with very different properties
of inclusions and the host matrix, including a partial reec-
tion of diffusing particles from inclusions and a trapping
inside inclusions. The more so, it could be useful in quite
different elds like mathematical modeling of the release of
antimicrobial agent from packaging material to a food
product [21], transport through membranes, or self-diffusion
of small molecules in colloidal systems containing poly-
mers, proteins, micells, etc. [7].
Acknowledgements
This study has been partly supported by the European
Center of Excellence in Advanced Materials Research and
Technology in Riga, Latvia (contract no. ICA-I-CT-2000-
7007).
References
[1] J.C. Maxwell-Garnett, Philos. Trans. R. Soc. London 203
(1904) 385.
[2] G.W. Milton, Pure and Applied Mathematics XLIII (1990) 63
and references therein.
[3] N. Oshima, S. Nomura, J. Composite Materials 19 (1985) 287.
[4] D. Bergman, Phys. Rep. 43 (1978) 377.
[5] R.M. Christensen, Mechanics of Composite Materials, Wiley-
Interscience, New York, 1979.
[6] R.M. Christensen, K.H. Lo, Journal of Mechanics and Physics
of Solids 27 (1979) 315.
[7] B. Jonsson, H. Wennerstrom, P.G. Nilsson, P. Line, Colloidal
and Polymer Science 264 (1986) 77.
[8] S. Torquato, J. Stat. Phys. 65 (1991) 1173.
[9] P.W.M. Jacobs, I.E. Hooton, J. Phys. Chem. Sol. 51 (1990)
1207.
[10] C.R.A. Catlow, Phil. Mag. A 64 (1991) 1011.
[11] Yu.R. Kalnin, Proc. Latvian Acad. Sci. Physics 6 (1985) 13.
J.R. Kalnin et al. / Journal of Physics and Chemistry of Solids 63 (2002) 449456 455
Fig. 6. Comparison of generally accepted Maxwell-Garnett theory (dotted line, Eq. (14)), and our Eq. (28) (solid line) with 2-D computer
simulations for periodically distributed reecting inclusions versus their dimensionless concentration; square symbols are for the spherical
inclusions on a square lattice; triangles for the spherical inclusions on a hexagonal lattice; and full circles are for square inclusions on a square
lattice.
[12] B.U. Felderhof, J.M. Deutch, J. Chem. Phys. 64 (1976) 4551.
[13] Yu.H. Kalnin, P. Zapol, Radiation Effects and Defects in
Solids 137 (1995) 295.
[14] J.R. Manning, Diffusion Kinetics for Atoms in Crystals, Prin-
ceton, Toronto, 1968.
[15] G.M. Ronis, Phys. Rev. A 36 (1987) 1908.
[16] R.I. Cukier, J. Phys. Chem. Sol. 87 (1983) 582.
[17] G.M. Bell, Trans. Farad. Soc. 60 (1961) 1752.
[18] Yu.R. Kalnin, E.A. Kotomin, J. Phys: Math Gen. 31 (1998)
7227.
[19] J. Vucans, Proceedings of IFIP WG 7.2 on Modelling and
Optimisation of Distributed Parameter Systems with Applica-
tions to Engineering, Chapman and Hall, London, 1995, p. 321.
[20] P.M. Smith, M.J. Aziz, Acta Metallurgica and Materialia 42
(1994) 3515.
[21] J.H. Han, Foodtechnology 54 (2000) 56.
J.R. Kalnin et al. / Journal of Physics and Chemistry of Solids 63 (2002) 449456 456