A Monte Carlo technique that simulates tracer diffusion in multiphase materials of arbitrary complexity has been developed. Effective diffusivities are calculated for structures consisting of either overlapping or nonoverlapping inclusions with diffusivity Dc, distributed in a continuous phase with diffusivity D0\u3eDc. Two-dimensional simulations for various values of D0/Dc generate normalized diffusivities that correspond closely to their three-dimensional counterparts; they nearly collapse to a common curve when a simple scaling relation is applied

Buettner, Helen M.

Muzzio, Fernando J.

Reyes, Sebastian C.

Riley, Mark R

English

UNL | Libraries

University of Nebraska - Lincoln
DigitalCommons@University of Nebraska - Lincoln
Biological Systems Engineering: Papers and
Publications Biological Systems Engineering
1994
Monte Carlo calculation of effective diffusivities in
two- and three-dimensional heterogeneous
materials of variable structure
Mark R. Riley
University of Nebraska-Lincoln, mriley3@unl.edu
Fernando J. Muzzio
Rutgers University - New Brunswick/Piscataway
Helen M. Buettner
Rutgers University - New Brunswick/Piscataway
Sebastian C. Reyes
Corporate Research Laboratories, Exxon Research and Engineering Company, Route 22 East, Annandale, New Jersey
Follow this and additional works at: http://digitalcommons.unl.edu/biosysengfacpub
Part of the Bioresource and Agricultural Engineering Commons, and the Civil and
Environmental Engineering Commons
This Article is brought to you for free and open access by the Biological Systems Engineering at DigitalCommons@University of Nebraska - Lincoln. It
has been accepted for inclusion in Biological Systems Engineering: Papers and Publications by an authorized administrator of
DigitalCommons@University of Nebraska - Lincoln.
Riley, Mark R.; Muzzio, Fernando J.; Buettner, Helen M.; and Reyes, Sebastian C., "Monte Carlo calculation of effective diffusivities in
two- and three-dimensional heterogeneous materials of variable structure" (1994). Biological Systems Engineering: Papers and
Publications. 360.
http://digitalcommons.unl.edu/biosysengfacpub/360
PHYSICAL REVIEW E VOLUME 49, NUMBER 4 APRIL 1994
Monte Carlo calculation of efFective difFusivities in two- and three-dimensional
heterogeneous materials of variable structure
Mark R. Riley, ' Fernando J. Muzzio, '* Helen M. Buettner, ' and Sebastian C. Reyes
Department of Chemical and Biochemical Engineering, Rutgers Uniuersity, P.O. Box 909, Piscataway, New Jersey 08855 090-9
Corporate Research Laboratories, Exxon Research and Engineering Company, Route 22 East, Annandale, Pew Jersey 08801
{Received 30 August 1993;revised manuscript received 17 December 1993)
A Monte Carlo technique that simulates tracer diffusion in multiphase materials of arbitrary complex-
ity has been developed. Effective diffusivities are calculated for structures consisting of either overlap-
ping or nonoverlapping inclusions with diffusivity D, distributed in a continuous phase with diffusivity
Do&D, . Two-dimensional simulations for various values of Do/D, generate normalized diffusivities
that correspond closely to their three-dimensional counterparts; they nearly collapse to a common curve
when a simple scaling relation is applied.
PACS number{s): 51.10.+y, 66.30.Jt, 05.60.+w, 05.40.+j
Multiphase materials occur in a wide variety of sys-
tems, including polymer blends, ceramic precursors, soil,
heterogeneous catalysts, and many types of biological tis-
sue. These systems are often composed of one or more
dispersed phases distributed in a single continuous phase.
In general, the effective diffusivity of molecules through
such materials depends on the volume fraction and solute
mobility in each phase as well as on morphological details
of the individual phases. Here, we present results of a
Monte Carlo method that can be used to evaluate
effective diffusivities of molecules in two-dimensional
(2D) and three-dimensional (3D) multiphase materials
with virtually any type of structure.
A substantial amount of experimental diffusivity data
is available in the literature, but extrapolation of these
measurements to other systems and conditions is general-
ly dif5cult. Several approaches have previously been used
by other researchers to predict diffusion in multiphase
materials. Theoretical methods have led to a number of
models, but the applicability of these models is frequently
limited by restrictions on the geometric structure and
volume fraction of the inclusions. For example,
Maxwell's equation [l] assumes that the discrete phase
occupies a small fraction of the total volume and consists
of nonoverlapping, spherical inclusions randomly distri-
buted in the continuous phase. In the context of
diffusivities, Maxwell s equation can be written as
D, 2Do+D, —2$(Do D,)—
Do 2DO+D, +$(Do D,)
where D, is the macroscopic effective diffusivity, Do is
the molecular diffusivity of the continuous phase, and D,
and P are the molecular diffusivity and the volume frac-
tion of the discrete phase, respectively.
Statistical approaches have been used to generate
bounds for D, (see, for example, [2]). However, the gap
'To whom correspondence should be addressed.
between bounds can be quite large, particularly when D,
differs significantly from Do. A convenient approach is
provided by Monte Carlo techniques, which have previ-
ously been used to study transport properties in a variety
of macrostructures composed of inclusions that were ei-
ther nonconductive [3] or more conductive than the con-
tinuous phase [4,5]. This paper focuses on a different
case where the inclusions are permeable but their trans-
port rates are lower than that of the continuous phase.
Such a case has not been studied in detail before although
it is relevant to many applications.
The first step in this Monte Carlo technique is to gen-
erate a computer representation of the macrostructure.
The methodology is illustrated in this paper by generat-
ing two types of structure composed of either (i) overlap-
pings inclusions or (ii) nonoverlapping inclusions distri-
buted throughout a continuous phase. Circular in-
clusions are placed in a square region to simulate a 2D
structure; spherical inclusions are placed in a cubic re-
gion to generate a 3D system. The placement region is
spatially periodic, so that if an inclusion intersects the
edge of the region, it partially reappears on the opposite
side to account for the overhanging area. The structure
with overlapping inclusions is produced by randomly po-
sitioning inclusions one at a time in the computational
domain. If two (or more) inclusions overlap a given area,
then this area is counted only once toward the total frac-
tional coverage. To produce the nonoverlapping struc-
ture, inclusions are again placed one at a time until the
desired inclusion fraction is attained. However, if the po-
sition of a new inclusion overlaps that of a previously
placed inclusion, the new inclusion is rejected and anoth-
er inclusion is placed in a different location. The process
is repeated until the new inclusion is placed in an entirely
open region. A small region of a representative 3D struc-
ture of nonoverlapping inclusions is depicted in Fig. l.
As noted by other researchers [4—6], this placement
mechanism results in a nonequihbrium structure. How-
ever, as will be described shortly, the placement metho-
dology has a relatively minor effect on the diffusive phe-
nomena investigated here.
1063-651X/94/49{4)/3500{4}/$06.00 49 3500 1994 The American Physical Society
49 BRIEF REPORTS 3501
y
K+.
pp. '~ ~ ~
'F
~yg~i
:::: III"'I' "",:;
y'
FIG. 1. Simulated 3D structure composed of nonoverlapping
spheres randomly placed in a continuous domain. For this
structure, the inclusion volume fraction is /=0. 30.
To determine the diffusivity of molecules in these
structures, tracer particles undergo a simulated random
walk through the structural representation. Molecular
diffusion in the continuous phase is described as a series
of randomly directed steps with length A, =A,Dine, , where e
is a random number with uniform probability in the in-
terval [0,1] and A,o is the mean step length in the continu-
ous phase. Particles inside the inclusions take steps of
length A, = —yk~lns, where y =D, /Do is a constant less
than 1, refiecting the decrease in diffusivity inside the in-
clusions (D, ) relative to that in the continuous phase
(Do). The effective diffusivity is calculated by monitoring
the distance traveled by the tracers in a certain amount of
time. The straight-line displacement X between the ini-
tial and final positions of each tracer is used to calculate
the mean-squared displacement (X ) for a large number
of tracers. For sufficiently many steps, the effective
diffusivity is then calculated using the relationship
(x')
2d+A, ' (2)
where d represents the dimensionality of the system and
gA, , the total length of the walk, is proportional to the to-
tal time along the tracer's trajectory.
Because tracers in each phase are assigned different
mean step lengths to account for different diffusivities in
each phase, it is necessary to know at all times in which
phase the tracer is located. To accomplish this, the simu-
lation domain is discretized into a logical matrix array of
one million nodes (1000X1000 for 2D simulations and
100X100X100 for 3D simulations). The lattice is used
only to discretize the spatial structure; the centers of the
spheres are not regularly positioned on the lattice. Nodes
within one inclusion radius from the center of an in-
clusion are assigned a logical value .TRUE, and all other
nodes are assigned FALSE.. To determine the phase in
which the tracer is located, the algorithm simply recalls
the value of the node nearest the tracer; a returned value
of TRUE. means the particle is inside an inclusion,
while a value of .FALSE. means the particle is in the
continuous phase (a related approach was used by Kim
and Torquato [5] to compute diffusivities in systems of
partially overlapping highly conductive spheres). Al-
though the structure chosen here for illustration is rather
simple, the above discretization procedure is quite flexible
and can be used to investigate essentially any type of
structure because the procedure does not require in-
clusions with a regular geometry. Inclusions with vari-
able morphology (spheroids, cylinders, irregularly shaped
inclusions) or spatial arrangement (distinct, overlapping,
or clustered inclusions) can also be readily described;
clusters of coalesced spheres have been successfully simu-
lated by this method [7]. Potentially, this technique
could also be used to simulate diffusion through real
structures by using digitized photographs of materials to
generate the computational domain for simulations. This
method has the additional advantage that 3D simulations
require a similar amount of computer time as 2D simula-
tions. This behavior is uncommon; as it is well known,
computational time usually increases significantly with
the dimensionality of the system.
Effective diffusivity results for freely overlapping circu-
lar inclusions in 2D are presented in Fig. 2(a). Results
are displayed as normalized diffusivities (D, /Do) for in-
clusion fractions in the range 0 & P & 1 and diffusivity ra-
tios (Do/D, ) ranging from 2 to 5. Results for uniformly
distributed nonoverlapping inclusions in 2D and 3D sys-
tems are presented in Fig. 2(b) and 2(c), respectively, for
inclusion fractions in the range 0&/&0. 50 (2D), and
0& P &0.45 (3D) for diffusivity ratios from 2 to 5. These
ranges of Do/D, and P span experimental values report-
ed for diffusivity measurements in systems of biological
cells encapsulated in gelatinous supports [8]. Each point
in Fig. 2 represents the average of ten structural realiza-
tions; each run was performed using the same number of
inclusions placed in different positions. Each set of sym-
bols corresponds to a given diffusivity ratio for varying
inclusion fractions. Because De/D, ) 1, effective
diffusivities decrease with increasing inclusion fraction.
The diffusivity also decreases as the inclusion diffusivity
decreases; i.e., larger diffusivity ratios correspond to
lower effective diffusivities D, .
All the curves in Fig. 2(a) appear to have the same
shape, except that they approach different lower limits as$~1. This observation can be quantified by using a sim-
ple scaling relationship:
[1 D, (P)/Do]—
(1 D, /Do)— (3)
This scaling removes the dependence on the diffusivity
ratio by stretching the curves so that the bounds for /=0
and P = 1 coincide for all sets of results [i.e.,
D (0)=0, D*(1)=1for all diffusivity ratios]. Rescaling
the curves in Fig. 2(a) according to Eq. (3) makes the re-
sults nearly collapse onto a single curve [Fig. 3(a)]. For
the overlapping inclusion structure, D * is closely approx-
imated for a11 sets of results by the following polynomial
fit (solid curve in Fig. 3):
D =1.7271$—0.8177$ +0.09075$
BRIEF REPORTS
This scaling technique also seems to work for nonover-
lapping inclusions. As shown in Fig. 3(b), once the
curves from Fig. 2(b) are rescaled according to Eq. (3), re-
sults corresponding to diferent diff'usivity ratios again
nearly collapse onto a single curve. Moreover, the same
polynomial fit (4) as in Fig. 3(a) is used to represent re-
scaled diffusivities for nonoverlapping structures [Fig.
3(b)]. For the conditions investigated, the material struc-
ture does not a8'ect D, as the results for overlapping and
nonoverlapping inclusions in 2D are nearly identical for a
1.0:
0.8
0.6
0.4
0.6
Q 0.4
0.2
0.0
0.0 0.2 0.4 0.6
(Inclusion fraction)
00-J
0.0 0.4 0.6
(Inclusion fraction)
1.0
1.0 1.0
0.8 0.8
0.6 ) 0.6
Z
0.4
Ch
0.2 ().2
0.0
0.0
I
0.2 {}4 0.6
I
0.8 1,0
0.0:.
0.0 0.4 0.6 0.8
1.0
(Inclusion fraction) (Inclusion fraction)
0.6 0.6
8
0.4
0.2
0.0
0.0 0.4 0.6
(Inclusion Iraction)
FIG. 2. Normalized effective diffusivities in 2D and 3D sys-
tems. Each point represents the average of ten simulation runs.
Curves are polynomial least-squares fits to the results. (a) 2D
diffusivities for overlapping inclusions; (b) 2D diffusivities for
nonoverlapping inclusions; (c) 3D diffusivities for nonoverlap-
ping inclusions.
0.0 "
1.00.80.20.0 0.4 0.6
(Inclusion fraction)
FIG. 3. Scaled diffusivities D* of normalized diffusivities
displayed in Fig. 2. The solid line corresponds to a third-order
least-squares polynomial fit to the results. (a) 2D scaled
diffusivities for overlapping inclusions; (b) 2D scaled diffusivities
for nonoverlapping inclusions; (b) 3D scaled diffusivities for
nonoverlapping inclusions.
49 BRIEF REPORTS 3503
D, /Do = 1 —(1 D, /Do )( l.7—271$—0.8177ttt
+0.09075$ ) . (5)
This empirical relation can be used to predict effective
diffusivities as a function of the inclusion fraction and the
molecular diffusivities in the individual phases both for
overlapping and nonoverlapping structures in 2D and for
nonoverlapping structures in 3D. The predictions of Eq.
(5) correspond to the three curves in Fig. 4. The simula-
tion results have also been extensively compared with ex-
given Do/D, and tI).
When the scaling Eq. (3) is applied to the 3D diffusivity
predictions [Fig. 3(c)], results for Do/D, =3,4, 5 again
nearly collapse onto one curve. The poorest agreement is
observed for Do/D, =2. For Do/D, ~1, D* ap-
proaches a straight line D *=/ [this limit follows directly
from Maxwell's model, Eq. (1)]. The case of Do/D, =2
can be regarded as a transition state between the
straight-line limit and the fully developed curves for
Do/D, =3,4, 5 (and possibly beyond); rescaled results for
Do/D, =3,4, 5 are essentially identical to one another.
Surprisingly, the 3D results are well represented again by
exactly the same polynomial Pt as the 2D results. Equa-
tion (4), therefore, closely predicts all the effective
diffusivity values with less than 2 —3 % error for
Do/D, =3, 4, or 5 for overlapping and nonoverlapping
structures both in 2D and in 3D.
The similarity of the normalized diffusivities is further
con6rmed by comparing the u'nscaled simulation results
in 2D and 3D both for overlapping and nonoverlapping
inclusions (Fig. 4). For Do/D, =3 and 5, results from 2D
and 3D simulations are almost identical. Results for
Do/D, =4 (not shown for the sake of clarity of display)
are also in excellent agreement with one another. Even
for Do/D, =2, results from 2D and 3D simulations are
quite similar both for overlapping and nonoverlapping
cases. This agreement is entirely empirical and we do not
presume that the correspondence is universal. This ob-
servation, however, has important practical implications,
since it suggests that for some systems, 2D simulations
can be used to predict diffusivities in 3D systems.
Finally, we note that Eqs. (3) and (4) can be combined
to yield
0.8
0.6
0.4
0.2
0.0
0.0
I
0.2 0.4 0.6
(Inclusion fraction)
0.8 1.0
FIG. 4. Comparison of the effective diffusivity in 2D and 3D.
Symbols are as in Figs. 2 and 3.
ACKNOWLEDGMENTS
This work was supported by grants from the Exxon
Foundation and the Merck Foundation to F.J.M. and by
a grant from the New Jersey Center for Biomaterials and
Medical Devices to H.M.B. M.R.R. gratefully acknowl-
edges Exxon Research and Engineering Company for
financial support.
perimental data [7] and efficiently predict the trends in
D, as a function of P and Do/D, .
In summary, a Monte Carlo algorithm has been
developed to calculate effective diffusivities of tracer par-
ticles in multiphase materials with virtually any type of
structure. Results were illustrated for materials consist-
ing of overlapping and nonoverlapping spheres randomly
distributed in a continuous phase of higher diffusivity;
other structures are considered elsewhere [7]. The pre-
dicted normalized diffusivities are largely independent of
the dimensionality of the system; both for 2D and 3D
simulations, diffusivities nearly collapse to the same curve
when scaled using Eq. (3). An empirical relation
developed for the prediction of diffusivities in 2D and 3D
systems works reasonably well for all the cases considered
here.
[1]J. C. Maxwell, A Treatise on Electricity and Magnetism,
3rd ed. (Dover, New York, 1954), Vol. 1.
[2] Z. Hashin and S. Shtrikman, J. Appl. Phys. 33, 3125
(1962).
[3]J. W. Evans, M. H. Abbasi, and A. Sarin, J. Chem. Phys.
72, 2967 (1980); K. A. Akanni and J. W. Evans, Chem.
Eng. Sci. 42, 1945 (1987);S. C. Reyes, E. Iglesia, and Y. C.
Chiew, Mat. Res. Soc. Symp. Proc. 195, 553 (1990); S. C.
Reyes and E. Iglesia, J. Catal. 129, 457 (1991).
[4] S. Y. Sheu, S. Kumar, and R. I. Cukier, Phys. Rev. B 42,
1431 (1990).
[5] I. C. Kim and S. Torquato, J.Appl. Phys. 68, 3892 (1990).
[6] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A.
N. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
[7] M. R. Riley, F. J. Muzzio, H. M. Buettner, and S. C.
Reyes (unpublished).
[8] B. A. Westrin and A. Axelsson, Biotechnol. Bioeng. 38,
439 (1991).



Monte Carlo calculation of effective diffusivities in two- and three-dimensional heterogeneous materials of variable structure

DigitalCommons@University of Nebraska

University of Nebraska - LincolnDigitalCommons@University of Nebraska - LincolnBiological Systems Engineering: Papers andPublications Biological Systems Engineering1994Monte Carlo calculation of effective diffusivities intwo- and three-dimensional heterogeneousmaterials of variable structureMark R. RileyUniversity of Nebraska-Lincoln, mriley3@unl.eduFernando J. MuzzioRutgers University - New Brunswick/PiscatawayHelen M. BuettnerRutgers University - New Brunswick/PiscatawaySebastian C. ReyesCorporate Research Laboratories, Exxon Research and Engineering Company, Route 22 East, Annandale, New JerseyFollow this and additional works at: http://digitalcommons.unl.edu/biosysengfacpubPart of the Bioresource and Agricultural Engineering Commons, and the Civil andEnvironmental Engineering CommonsThis Article is brought to you for free and open access by the Biological Systems Engineering at DigitalCommons@University of Nebraska - Lincoln. Ithas been accepted for inclusion in Biological Systems Engineering: Papers and Publications by an authorized administrator ofDigitalCommons@University of Nebraska - Lincoln.Riley, Mark R.; Muzzio, Fernando J.; Buettner, Helen M.; and Reyes, Sebastian C., "Monte Carlo calculation of effective diffusivities intwo- and three-dimensional heterogeneous materials of variable structure" (1994). Biological Systems Engineering: Papers andPublications. 360.http://digitalcommons.unl.edu/biosysengfacpub/360PHYSICAL REVIEW E VOLUME 49, NUMBER 4 APRIL 1994Monte Carlo calculation of efFective difFusivities in two- and three-dimensionalheterogeneous materials of variable structureMark R. Riley, ' Fernando J. Muzzio, '* Helen M. Buettner, ' and Sebastian C. ReyesDepartment of Chemical and Biochemical Engineering, Rutgers Uniuersity, P.O. Box 909, Piscataway, New Jersey 08855 090-9Corporate Research Laboratories, Exxon Research and Engineering Company, Route 22 East, Annandale, Pew Jersey 08801{Received 30 August 1993;revised manuscript received 17 December 1993)A Monte Carlo technique that simulates tracer diffusion in multiphase materials of arbitrary complex-ity has been developed. Effective diffusivities are calculated for structures consisting of either overlap-ping or nonoverlapping inclusions with diffusivity D, distributed in a continuous phase with diffusivityDo&D, . Two-dimensional simulations for various values of Do/D, generate normalized diffusivitiesthat correspond closely to their three-dimensional counterparts; they nearly collapse to a common curvewhen a simple scaling relation is applied.PACS number{s): 51.10.+y, 66.30.Jt, 05.60.+w, 05.40.+jMultiphase materials occur in a wide variety of sys-tems, including polymer blends, ceramic precursors, soil,heterogeneous catalysts, and many types of biological tis-sue. These systems are often composed of one or moredispersed phases distributed in a single continuous phase.In general, the effective diffusivity of molecules throughsuch materials depends on the volume fraction and solutemobility in each phase as well as on morphological detailsof the individual phases. Here, we present results of aMonte Carlo method that can be used to evaluateeffective diffusivities of molecules in two-dimensional(2D) and three-dimensional (3D) multiphase materialswith virtually any type of structure.A substantial amount of experimental diffusivity datais available in the literature, but extrapolation of thesemeasurements to other systems and conditions is general-ly dif5cult. Several approaches have previously been usedby other researchers to predict diffusion in multiphasematerials. Theoretical methods have led to a number ofmodels, but the applicability of these models is frequentlylimited by restrictions on the geometric structure andvolume fraction of the inclusions. For example,Maxwell's equation [l] assumes that the discrete phaseoccupies a small fraction of the total volume and consistsof nonoverlapping, spherical inclusions randomly distri-buted in the continuous phase. In the context ofdiffusivities, Maxwell s equation can be written asD, 2Do+D, —2$(Do D,)—Do 2DO+D, +$(Do D,)where D, is the macroscopic effective diffusivity, Do isthe molecular diffusivity of the continuous phase, and D,and P are the molecular diffusivity and the volume frac-tion of the discrete phase, respectively.Statistical approaches have been used to generatebounds for D, (see, for example, [2]). However, the gap'To whom correspondence should be addressed.between bounds can be quite large, particularly when D,differs significantly from Do. A convenient approach isprovided by Monte Carlo techniques, which have previ-ously been used to study transport properties in a varietyof macrostructures composed of inclusions that were ei-ther nonconductive [3] or more conductive than the con-tinuous phase [4,5]. This paper focuses on a differentcase where the inclusions are permeable but their trans-port rates are lower than that of the continuous phase.Such a case has not been studied in detail before althoughit is relevant to many applications.The first step in this Monte Carlo technique is to gen-erate a computer representation of the macrostructure.The methodology is illustrated in this paper by generat-ing two types of structure composed of either (i) overlap-pings inclusions or (ii) nonoverlapping inclusions distri-buted throughout a continuous phase. Circular in-clusions are placed in a square region to simulate a 2Dstructure; spherical inclusions are placed in a cubic re-gion to generate a 3D system. The placement region isspatially periodic, so that if an inclusion intersects theedge of the region, it partially reappears on the oppositeside to account for the overhanging area. The structurewith overlapping inclusions is produced by randomly po-sitioning inclusions one at a time in the computationaldomain. If two (or more) inclusions overlap a given area,then this area is counted only once toward the total frac-tional coverage. To produce the nonoverlapping struc-ture, inclusions are again placed one at a time until thedesired inclusion fraction is attained. However, if the po-sition of a new inclusion overlaps that of a previouslyplaced inclusion, the new inclusion is rejected and anoth-er inclusion is placed in a different location. The processis repeated until the new inclusion is placed in an entirelyopen region. A small region of a representative 3D struc-ture of nonoverlapping inclusions is depicted in Fig. l.As noted by other researchers [4—6], this placementmechanism results in a nonequihbrium structure. How-ever, as will be described shortly, the placement metho-dology has a relatively minor effect on the diffusive phe-nomena investigated here.1063-651X/94/49{4)/3500{4}/$06.00 49 3500 1994 The American Physical Society49 BRIEF REPORTS 3501yK+.pp. '~ ~ ~'F~yg~i:::: III"'I' "",:;y'FIG. 1. Simulated 3D structure composed of nonoverlappingspheres randomly placed in a continuous domain. For thisstructure, the inclusion volume fraction is /=0. 30.To determine the diffusivity of molecules in thesestructures, tracer particles undergo a simulated randomwalk through the structural representation. Moleculardiffusion in the continuous phase is described as a seriesof randomly directed steps with length A, =A,Dine, , where eis a random number with uniform probability in the in-terval [0,1] and A,o is the mean step length in the continu-ous phase. Particles inside the inclusions take steps oflength A, = —yk~lns, where y =D, /Do is a constant lessthan 1, refiecting the decrease in diffusivity inside the in-clusions (D, ) relative to that in the continuous phase(Do). The effective diffusivity is calculated by monitoringthe distance traveled by the tracers in a certain amount oftime. The straight-line displacement X between the ini-tial and final positions of each tracer is used to calculatethe mean-squared displacement (X ) for a large numberof tracers. For sufficiently many steps, the effectivediffusivity is then calculated using the relationship(x')2d+A, ' (2)where d represents the dimensionality of the system andgA, , the total length of the walk, is proportional to the to-tal time along the tracer's trajectory.Because tracers in each phase are assigned differentmean step lengths to account for different diffusivities ineach phase, it is necessary to know at all times in whichphase the tracer is located. To accomplish this, the simu-lation domain is discretized into a logical matrix array ofone million nodes (1000X1000 for 2D simulations and100X100X100 for 3D simulations). The lattice is usedonly to discretize the spatial structure; the centers of thespheres are not regularly positioned on the lattice. Nodeswithin one inclusion radius from the center of an in-clusion are assigned a logical value .TRUE, and all othernodes are assigned FALSE.. To determine the phase inwhich the tracer is located, the algorithm simply recallsthe value of the node nearest the tracer; a returned valueof TRUE. means the particle is inside an inclusion,while a value of .FALSE. means the particle is in thecontinuous phase (a related approach was used by Kimand Torquato [5] to compute diffusivities in systems ofpartially overlapping highly conductive spheres). Al-though the structure chosen here for illustration is rathersimple, the above discretization procedure is quite flexibleand can be used to investigate essentially any type ofstructure because the procedure does not require in-clusions with a regular geometry. Inclusions with vari-able morphology (spheroids, cylinders, irregularly shapedinclusions) or spatial arrangement (distinct, overlapping,or clustered inclusions) can also be readily described;clusters of coalesced spheres have been successfully simu-lated by this method [7]. Potentially, this techniquecould also be used to simulate diffusion through realstructures by using digitized photographs of materials togenerate the computational domain for simulations. Thismethod has the additional advantage that 3D simulationsrequire a similar amount of computer time as 2D simula-tions. This behavior is uncommon; as it is well known,computational time usually increases significantly withthe dimensionality of the system.Effective diffusivity results for freely overlapping circu-lar inclusions in 2D are presented in Fig. 2(a). Resultsare displayed as normalized diffusivities (D, /Do) for in-clusion fractions in the range 0 & P & 1 and diffusivity ra-tios (Do/D, ) ranging from 2 to 5. Results for uniformlydistributed nonoverlapping inclusions in 2D and 3D sys-tems are presented in Fig. 2(b) and 2(c), respectively, forinclusion fractions in the range 0&/&0. 50 (2D), and0& P &0.45 (3D) for diffusivity ratios from 2 to 5. Theseranges of Do/D, and P span experimental values report-ed for diffusivity measurements in systems of biologicalcells encapsulated in gelatinous supports [8]. Each pointin Fig. 2 represents the average of ten structural realiza-tions; each run was performed using the same number ofinclusions placed in different positions. Each set of sym-bols corresponds to a given diffusivity ratio for varyinginclusion fractions. Because De/D, ) 1, effectivediffusivities decrease with increasing inclusion fraction.The diffusivity also decreases as the inclusion diffusivitydecreases; i.e., larger diffusivity ratios correspond tolower effective diffusivities D, .All the curves in Fig. 2(a) appear to have the sameshape, except that they approach different lower limits as$~1. This observation can be quantified by using a sim-ple scaling relationship:[1 D, (P)/Do]—(1 D, /Do)— (3)This scaling removes the dependence on the diffusivityratio by stretching the curves so that the bounds for /=0and P = 1 coincide for all sets of results [i.e.,D (0)=0, D*(1)=1for all diffusivity ratios]. Rescalingthe curves in Fig. 2(a) according to Eq. (3) makes the re-sults nearly collapse onto a single curve [Fig. 3(a)]. Forthe overlapping inclusion structure, D * is closely approx-imated for a11 sets of results by the following polynomialfit (solid curve in Fig. 3):D =1.7271$—0.8177$ +0.09075$BRIEF REPORTSThis scaling technique also seems to work for nonover-lapping inclusions. As shown in Fig. 3(b), once thecurves from Fig. 2(b) are rescaled according to Eq. (3), re-sults corresponding to diferent diff'usivity ratios againnearly collapse onto a single curve. Moreover, the samepolynomial fit (4) as in Fig. 3(a) is used to represent re-scaled diffusivities for nonoverlapping structures [Fig.3(b)]. For the conditions investigated, the material struc-ture does not a8'ect D, as the results for overlapping andnonoverlapping inclusions in 2D are nearly identical for a1.0:0.80.60.40.6Q 0.40.20.00.0 0.2 0.4 0.6(Inclusion fraction)00-J0.0 0.4 0.6(Inclusion fraction)1.01.0 1.00.8 0.80.6 ) 0.6Z0.4Ch0.2 ().20.00.0I0.2 {}4 0.6I0.8 1,00.0:.0.0 0.4 0.6 0.81.0(Inclusion fraction) (Inclusion fraction)0.6 0.680.40.20.00.0 0.4 0.6(Inclusion Iraction)FIG. 2. Normalized effective diffusivities in 2D and 3D sys-tems. Each point represents the average of ten simulation runs.Curves are polynomial least-squares fits to the results. (a) 2Ddiffusivities for overlapping inclusions; (b) 2D diffusivities fornonoverlapping inclusions; (c) 3D diffusivities for nonoverlap-ping inclusions.0.0 "1.00.80.20.0 0.4 0.6(Inclusion fraction)FIG. 3. Scaled diffusivities D* of normalized diffusivitiesdisplayed in Fig. 2. The solid line corresponds to a third-orderleast-squares polynomial fit to the results. (a) 2D scaleddiffusivities for overlapping inclusions; (b) 2D scaled diffusivitiesfor nonoverlapping inclusions; (b) 3D scaled diffusivities fornonoverlapping inclusions.49 BRIEF REPORTS 3503D, /Do = 1 —(1 D, /Do )( l.7—271$—0.8177ttt+0.09075$ ) . (5)This empirical relation can be used to predict effectivediffusivities as a function of the inclusion fraction and themolecular diffusivities in the individual phases both foroverlapping and nonoverlapping structures in 2D and fornonoverlapping structures in 3D. The predictions of Eq.(5) correspond to the three curves in Fig. 4. The simula-tion results have also been extensively compared with ex-given Do/D, and tI).When the scaling Eq. (3) is applied to the 3D diffusivitypredictions [Fig. 3(c)], results for Do/D, =3,4, 5 againnearly collapse onto one curve. The poorest agreement isobserved for Do/D, =2. For Do/D, ~1, D* ap-proaches a straight line D *=/ [this limit follows directlyfrom Maxwell's model, Eq. (1)]. The case of Do/D, =2can be regarded as a transition state between thestraight-line limit and the fully developed curves forDo/D, =3,4, 5 (and possibly beyond); rescaled results forDo/D, =3,4, 5 are essentially identical to one another.Surprisingly, the 3D results are well represented again byexactly the same polynomial Pt as the 2D results. Equa-tion (4), therefore, closely predicts all the effectivediffusivity values with less than 2 —3 % error forDo/D, =3, 4, or 5 for overlapping and nonoverlappingstructures both in 2D and in 3D.The similarity of the normalized diffusivities is furthercon6rmed by comparing the u'nscaled simulation resultsin 2D and 3D both for overlapping and nonoverlappinginclusions (Fig. 4). For Do/D, =3 and 5, results from 2Dand 3D simulations are almost identical. Results forDo/D, =4 (not shown for the sake of clarity of display)are also in excellent agreement with one another. Evenfor Do/D, =2, results from 2D and 3D simulations arequite similar both for overlapping and nonoverlappingcases. This agreement is entirely empirical and we do notpresume that the correspondence is universal. This ob-servation, however, has important practical implications,since it suggests that for some systems, 2D simulationscan be used to predict diffusivities in 3D systems.Finally, we note that Eqs. (3) and (4) can be combinedto yield0.80.60.40.20.00.0I0.2 0.4 0.6(Inclusion fraction)0.8 1.0FIG. 4. Comparison of the effective diffusivity in 2D and 3D.Symbols are as in Figs. 2 and 3.ACKNOWLEDGMENTSThis work was supported by grants from the ExxonFoundation and the Merck Foundation to F.J.M. and bya grant from the New Jersey Center for Biomaterials andMedical Devices to H.M.B. M.R.R. gratefully acknowl-edges Exxon Research and Engineering Company forfinancial support.perimental data [7] and efficiently predict the trends inD, as a function of P and Do/D, .In summary, a Monte Carlo algorithm has beendeveloped to calculate effective diffusivities of tracer par-ticles in multiphase materials with virtually any type ofstructure. Results were illustrated for materials consist-ing of overlapping and nonoverlapping spheres randomlydistributed in a continuous phase of higher diffusivity;other structures are considered elsewhere [7]. The pre-dicted normalized diffusivities are largely independent ofthe dimensionality of the system; both for 2D and 3Dsimulations, diffusivities nearly collapse to the same curvewhen scaled using Eq. (3). An empirical relationdeveloped for the prediction of diffusivities in 2D and 3Dsystems works reasonably well for all the cases consideredhere.[1]J. C. Maxwell, A Treatise on Electricity and Magnetism,3rd ed. (Dover, New York, 1954), Vol. 1.[2] Z. Hashin and S. Shtrikman, J. Appl. Phys. 33, 3125(1962).[3]J. W. Evans, M. H. Abbasi, and A. Sarin, J. Chem. Phys.72, 2967 (1980); K. A. Akanni and J. W. Evans, Chem.Eng. Sci. 42, 1945 (1987);S. C. Reyes, E. Iglesia, and Y. C.Chiew, Mat. Res. Soc. Symp. Proc. 195, 553 (1990); S. C.Reyes and E. Iglesia, J. Catal. 129, 457 (1991).[4] S. Y. Sheu, S. Kumar, and R. I. Cukier, Phys. Rev. B 42,1431 (1990).[5] I. C. Kim and S. Torquato, J.Appl. Phys. 68, 3892 (1990).[6] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A.N. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).[7] M. R. Riley, F. J. Muzzio, H. M. Buettner, and S. C.Reyes (unpublished).[8] B. A. Westrin and A. Axelsson, Biotechnol. Bioeng. 38,439 (1991).

https://digitalcommons.unl.edu/cgi/viewcontent.cgi?article=1362&amp;context=biosysengfacpub

Monte Carlo calculation of effective diffusivities in two- and three-dimensional heterogeneous materials of variable structure

Abstract

Similar works

Full text

Available Versions

UNL | Libraries

DigitalCommons@University of Nebraska